Wenjun Zha1,2†Changyan Li1†Yan Wu1†
Junxiao Chen1†Sanhe Li1Minshan Sun3Bian Wu1
Shaojie Shi1Kai Liu1Huashan Xu1Peide Li1Kai Liu1Guocai Yang1Zhijun Chen1Deze Xu1,2*Lei Zhou1,2*
Aiqing You1,2*- 1Key Laboratory of Crop Molecular Breeding, Ministry of Agriculture and Rural Affairs, Hubei Key Laboratory of Food Crop Germplasm and Genetic Improvement, Food Crops Institute, Hubei Academy of Agricultural Sciences, Wuhan, China
- 2Hubei Hongshan Laboratory, Wuhan, China
- 3Henan Assist Research Biotechnology Co., Ltd., Zhengzhou, China
The brown planthopper (BPH) (Nilaparvata lugens) sucks rice sap causing leaves to turn yellow and wither, often leading to reduced or zero yields. Rice co-evolved to resist damage by BPH. However, the molecular mechanisms, including the cells and tissues, involved in the resistance are still rarely reported. Single-cell sequencing technology allows us to analyze different cell types involved in BPH resistance. Here, using single-cell sequencing technology, we compared the response offered by the leaf sheaths of the susceptible (TN1) and resistant (YHY15) rice varieties to BPH (48 hours after infestation). We found that the 14,699 and 16,237 cells (identified via transcriptomics) in TN1 and YHY15 could be annotated using cell-specific marker genes into nine cell-type clusters. The two rice varieties showed significant differences in cell types (such as mestome sheath cells, guard cells, mesophyll cells, xylem cells, bulliform cells, and phloem cells) in the rice resistance mechanism to BPH. Further analysis revealed that although mesophyll, xylem, and phloem cells are involved in the BPH resistance response, the molecular mechanism used by each cell type is different. Mesophyll cell may regulate the expression of genes related to vanillin, capsaicin, and ROS production, phloem cell may regulate the cell wall extension related genes, and xylem cell may be involved in BPH resistance response by controlling the expression of chitin and pectin related genes. Thus, rice resistance to BPH is a complicated process involving multiple insect resistance factors. The results presented here will significantly promote the investigation of the molecular mechanisms underlying the resistance of rice to insects and accelerate the breeding of insect-resistant rice varieties.
Introduction
Rice (Oryza sativa L.) is one of the most important cereal crops and has been domesticated for approximately 7,500 years (Zong et al., 2007). Stored and unharvested rice can be attacked by >800 species of insect pests (Ghaffar et al., 2011). The brown planthopper (BPH), Nilaparvata lugens (Stål), is one of the most economically important insects which can cause huge destruction of rice plants (Xue et al., 2014). BPH can damage rice growth by spreading plant viruses (such as rice grassy and rice-ragged stunt viruses) (Cabauatan et al., 2009) and sucking plant sap.
BPH has co-evolved to adapt strongly to its host, rice (Sezer and Butlin, 1998; Zheng et al., 2021). The widely used insecticides could effectively protect from the BPH and other pests-caused damages. However, overuse of insecticides not only promotes BPH obtaining the adaptability and resistance to insecticides but also causes serious environmental pollution (Senthil-Nathan et al., 2009). However, a complex defense system against BPH also exists in rice. Therefore, developing rice varieties resistant to insects using their insect-resistance genes could be an ideal complement and alternative to existing insect-control measures. Since the first BPH-resistant rice variety was discovered, >40 genes associated with BPH resistance (Akanksha et al., 2019; Li et al., 2019), such as bph1-40, have been identified. Some of these BPH-resistant genes, such as bph1-4, have been successfully used in breeding BPH-resistant rice varieties (Athwal et al., 1971; Laksminarayana and Khush, 1977). These BPH-resistance genes not only promote rice BPH-resistance but also decrease BPH reproduction and prolong the period of BPH development (Du et al., 2009; Senthil-Nathan et al., 2009; Nguyen et al., 2019). The interaction between BPH and rice is a complex and dynamic process. Currently, it is understood that BPH sucks phloem sap by inserting the stylet bundle with an accompanying salivary sheath into the plant (Spiller, 1990). On its way to the phloem, the mouth stylet bundle pierces through various cells, such as the epidermis, mesophyll, phloem, etc. Thus, multiple resistant mechanisms may interrupt the BPH feeding process at several cellular locations. Moreover, the different cells encountered by the stylet may mount different resistant functions depending on the rice variety. Presently, the identity of the cells involved in insect resistance in rice and the underlying molecular mechanisms remain unknown.
Single-cell RNA sequencing (scRNA-seq) provides a method to examine the expression of all genes in each cell at the transcriptional level (Brennecke et al., 2013; Islam et al., 2014). Recently, scRNA-seq was applied in many fields to explore the heterogeneity of cells (Luecken and Theis, 2019). In addition, studies on heterogeneity in animal cells have been abundantly reported (Butler et al., 2018; Aran et al., 2019; Zhang et al., 2019). However, because of the presence of plant cell walls, the application of scRNA-seq is limited in plants. Currently, several single-cell studies in plants focusing on tissue and cellular functional differentiation–such as differentiation of roots and stems in Arabidopsis (Zhang et al., 2021b), stomatal lineage and developing leaf (Lopez-Anido et al., 2021), ploidy-dependence in Arabidopsis female gametophytes (Song et al., 2020), leaf and root differentiation in rice (Liu Q. et al., 2021; Wang et al., 2021; Zhang et al., 2021a), differentiation of maize ears facilitates (Xu et al., 2021b), poplar xylem formation (Xie et al., 2022), Nicotiana attenuate corolla cells formation (Kang et al., 2022)–have been reported. However, studies on the molecular mechanisms underlying the response of different plant cells to biotic stresses are rarely reported.
Here, we used scRNA-seq to explore the differences in the molecular responses mounted by the various cell types to BPH biotic stress (to BPH infestation) in two rice varieties differing in their resistance to BPH. This study’s results will be a reference for future studies revealing the molecular mechanisms of rice insect resistance and provide a theoretical basis for better breeding of varieties resistant to BPH.
Materials and methods
Rice plants cultivation and infesting rice plants with brown planthopper
For this study, seedlings of the rice varieties TN1 (susceptible to BPH) and YHY15 (moderately resistant to BPH) were cultivated in a climate chamber under a 12-h light/12-h dark cycle at 30 ± 1°C and 70% relative humidity. Then, the 2- to 3-instar hopper nymphs were used to infest the seedlings at the third-leaf stage (13 or 14 days old) at a density of eight insects per seedling under the conditions of 70% relative humidity, 25 ± 1°C, and 16-h light/8-h dark cycle. And we covered the stems of each seedling with breathable plastic tubes to prevent the brown planthoppers from escaping. After 48 hours, the leaf sheaths, including herbivore-exposed local (damaged) parts of TN1 and YHY15, were collected for further study. Three seedlings were sampled for each rice variety.
Rice protoplast isolation
The rice protoplast was isolated using the previous methods (Jabnoune et al., 2015; Wang et al., 2021). Briefly, the finely cut rice seedling leaf sheaths were immediately incubated in an enzymes solution containing 10 mM MES, 0.6 M mannitol, 0.75% macerozyme R-10, and 1.5% cellulase RS (pH 5.7) for 3 h (shaking at 70 rpm) at 28°C. After incubation, the enzyme solution was filtered out. The digested protoplasts were washed with the W5 solution (154 mM NaCl, 125 mM CaCl2, 5 mM KCl, and 2 mM MES at pH 5.7) and collected by centrifugation at 300 × g. The protoplasts were resuspended in a solution containing 0.6 M D-Mannitol, 15 mM MgCl2, and 4 mM MES (pH 5.7). A small amount of the single-cell suspension was added to an equal volume of 0.4% trypan blue dye. The concentration of viable cells was adjusted to the desired concentration (1000 to 2000 cells/µL) by counting the cells using Countess® II Automated Cell Counter.
Single-cell RNA sequencing
To generate the single-cell Gel Bead-In-Emulsions (GEMs) from cellular suspensions, a GemCode Single-cell instrument (10x Genomics, Pleasanton, CA, USA) was used. The Chromium Next GEM Single Cell 3’ Reagent Kits v3.1 (CG000183, RevC, 10x Genomics) was used to prepare the library and for sequencing. To construct the library, the barcoded, reverse-transcribed, and full-length cDNAs were amplified using PCR assay. Next, the PCR products were ligated to the adapters for sequencing, followed by PCR amplification according to the DNA template concentrations. Then a paired-end sequencing was conducted using the Illumina platform (Illumina Inc., San Diego, USA) for sequencing libraries.
Preprocessing and cell clustering analysis
To generate counts quantification, alignment, and FASTQ files from raw Illumina BCL files, the Cell Ranger software (version 3.1.0, 10x Genomics) was used. (Lun et al., 2019). To align the sequencing, reads of rice scRNA-seq samples were aligned to the rice reference genome (Kawahara et al., 2013). The Seurat software (V3.1.1, Satija Lab, New York, USA) was used for the analysis of the downstream genes through importing the gene matrices cell for each sample individually (Butler et al., 2018; Stuart et al., 2019). Using DoubletFinder (v2.0.3), we filtered out the cells with doublet GEMs, no less than 8000 UMIs and no less than 10% mitochondrial genes (McGinnis et al., 2019). After normalizing the data, the harmony algorithm was used to correct the batch effect (Korsunsky et al., 2019). Principal component analysis dimensional reduction was applied as a reference according to the previous report (Chung and Storey, 2015). Cells clustering was analyzed using the Louvain method to maximize modularity (Rotta and Noack, 2011). The combination of cell type annotation with these previously reported cell type marker genes was conducted using the R packages SingleR (Aran et al., 2019) and Cellassign (Zhang et al., 2019).
Differentially expressed genes analysis
The Wilcoxon rank sum test was used to compare difference in expression of every detected gene between the given cluster and other cells (Camp et al., 2017). We identified the significantly upregulated genes using the following criteria–genes had to be at least 1.28-fold overexpressed in a target cluster having >25% of same type of cells, and have a p value < 0.01. Subsequently, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analyses and Gene Ontology (GO) functional annotation were conducted to analyze these DEGs (Boyle et al., 2004; Kanehisa and Goto, 2000; Kanehisa et al., 2008).
Gene set variation analysis and cell cycle analysis
GSVA was performed using a collection of gene sets from Molecular Signatures Database (MSigDB) (Liberzon et al., 2015) to identify pathways and cellular processes enriched in different clusters. GSVA was performed as implemented in the GSVA R package version 1.26 (Hanzelmann et al., 2013) based on the cluster-averaged log-transformed expression matrix. According to the expression of the genes associated with G1/S phase (n = 100), S phase (n = 113), G2/M phase (n = 133), M phase (n = 106), and M/G1 phase (n = 106), the cell cycle score was assigned using the package “Seurat” (Macosko et al., 2015). Cells with the highest score <0.3 were defined as non-cycling cells (Neftel et al., 2019). In addition, the Plant Transcription Factor Database (PlantTFDB) was used in the transcript factors annotation (Jin et al., 2017).
RNA in situ hybridization assay
The situ hybridization assay was conducted based on the previously published protocol (Umeda et al., 1999). Briefly, hydration and paraffin embedding of the fresh rice seedling leaf sheaths were performed after 12 h fixation with FAA solution (3.7% formaldehyde, 5% acetic acid, and 50% ethanol). Subsequently, the paraffin-embedded rice seedling leaf sheaths were sectioned into 10 µm thick sections and treated for 2 h at 62°C (KD-P, Zhejiang Jinhua Kedi Instrumental Equipment Co., Ltd, China) and xylene (twice for 15 min) to remove paraffin and serially rehydrated using different concentrations of ethanol. Then, the leaf sheath sections were hybridized with RNA probes following a 15-minute Proteinase K (20 µg/ml) (G1234, Nanjing Zoonbio Biotechnology, Ltd., China) digestion (at 37°C) and serial dehydration using different concentrations of ethanol. Next, the probes were transcribed in vitro using a Digoxigenin RNA labeling kit (Roche, USA). The transcribed probes were then incubated with the leaf sheath tissue sections. Finally, they were washed and incubated with an anti-digoxigenin-AP (200-052-156, Jackson ImmunoResearch Inc., PA, USA). The RNA hybridization signals were detected at room temperature by staining with a nitro-blue tetrazolium/5-bromo-4-chloro-30-indolyphosphate stock solution (NBT/BCIP solution; Boster Bio, CA, USA). Images were taken in the bright field mode using a microscope (Nikon Eclipse ci, Nikon Instruments Inc., NY, USA).
Results
Characteristics of the constructed single-cell transcriptome library of the leaf sheaths of BPH-resistant rice variety
We isolated protoplasts from rice leaf sheaths after a 48 h infestation with the BPH to generate a single-cell transcriptome of the rice resistant to BPH. The 10x Genomics Chromium and the Illumina sequencing platforms were used to generate scRNA-Seq libraries (Figure 1A). For TN1 and YHY15 samples, we obtained 23,346 and 19,775 reads per cell, respectively. 1,670 expressed genes and 4,600 unique molecular identifiers (UMI) were generated for each cell. In TN1 and YHY15, we also detected 27,740 and 26,710 genes, respectively (Supplementary Tables 1–3). Using t-distributed stochastic neighbor embedding (tSNE) projection, the single-cell transcriptomes were plotted, and largely overlapping distributions between TN1 and YHY15 were observed, suggesting a high reproducibility (Figure 1B). To categorize the single cells, the package “Seurat” was applied. The results showed 16 major clusters of cell transcriptomes in two-dimensional space (Figure 1C). The total cells for TN1 and YHY15 were 14699 and 16237, respectively. Among these 16 clusters, the proportion of cells in most clusters was not significantly different, but the proportion of YHY15 cells was relatively low in clusters 11–14. In contrast, in cluster 15, the cell ratio of YHY15 was higher (Supplementary Table 4). Next, we examined the genes significantly upregulated in the 16 clusters, and the top 5 highly expressed genes were selected in each cluster (Supplementary Table 5); the expression pattern of these genes showed apparent cluster specificity (Figure 1D).
Figure 1 Quality control of the rice seedlings scRNA-seq analysis. (A) Representative schematic graph showing the workflow of scRNA-seq; (B) Representative plot of dimensional reduction of TN1 and YHY15; [Orange- TN1, Blue- YHY15] (C) Representative plot of the 16 major clusters of the TN1 and YHY15 cell transcriptomes; (D) Top 5 genes exhibiting upregulated expression in 16 clusters. [Yellow- high expression, Purple- low expression].
Tissue-specific marker gene analysis of the BPH-inoculated rice leaf sheath transcriptome library revealed nine cell/tissue types
For the assignment of tissue/cell type to clusters, the accumulated transcripts in our single-cell population were analyzed for significant expression of 52 marker genes in leaf sheath tissue/cell types (Supplementary Table 6). Nine major cell-type clusters were observed in both TN1 and YHY15 (Figure 2A) transcriptomes. Through comparing the changes in cell proportions for each cell type, we found that the proportions of mesophyll did not change much, but the proportions of procambium, guard cell, mestome sheath cell, and phloem were quite different (Figure 2B). The expression distribution of some existing marker genes, such as—1) Mesophyll marker genes (LOC_Os07g38960 (CAB7), LOC_Os01g41710 (CAB2R), LOC_Os12g19470 (RBCS), and LOC_Os12g19381 (RBCS)); 2) Procambium marker genes (LOC_Os02g08100 (4CL3), LOC_Os12g04080 (TBT1), and LOC_Os11g42290 (TBT1)); 3) Bulliform marker genes (LOC_Os06g14540 (GLU13)); 4) Phloem marker genes (LOC_Os06g41090 (FTIP1), LOC_Os01g06500 (PP2A1), and LOC_Os03g07480 (SUT1)); 5) Guard cell marker genes (LOC_Os03g41460 (SPARK10) and LOC_Os04g48530 (SLAC1)); and 6) Mestome sheath marker genes (LOC_Os01g68540 (GDI1))—were analyzed in each cell-type cluster. The expression of these genes was consistent with previously reported (Figures 2C–P). Therefore, we focused on detecting expression patterns of three selected genes, viz., LOC_Os02g08100 (4CL3), LOC_Os01g41710 (CAB2R), and LOC_Os12g19381 (RBCS). The RNA in situ hybridization results confirmed the marker gene LOC_Os02g08100 (4CL3) expression in procambium cells at the center of immature young leaf sheath (Figure 2Q). In mesophyll cells, a high RNA abundance of the mesophyll marker genes–LOC_Os12g19381 (RBCS) and LOC_Os01g41710 (CAB2R)–was seen, indicating that the annotation of cell types was reliable (Figures 2R, S).
Figure 2 Single-cell transcriptome atlas for the rice leaf. (A) Representative plot of the combined accumulation of transcript from the tested marker genes (listed in Supplemental Table S1); (B) The percentage of cells in cell-type clusters; (C–P) tSNE plots of marker genes predicting the identities of clusters [The color scale indicates normalized expression level]; (C–F) Mesophyll maker genes, LOC_Os07g38960 (CAB7), LOC_Os01g41710 (CAB2R), LOC_Os12g19470 (RBCS), and LOC_Os12g19381 (RBCS); (G–I) Procambium maker genes, LOC_Os02g08100 (4CL3), LOC_Os12g04080 (TBT1), and LOC_Os11g42290 (TBT1); (J) Bulliform maker gene, LOC_Os06g14540 (GLU13); (K–M) Phloem maker genes, LOC_Os06g41090 (FTIP1), LOC_Os01g06500 (PP2A1), and LOC_Os03g07480 (SUT1); (N, O) Guard cell maker genes, LOC_Os03g41460 (SPARK10), LOC_Os04g48530 (SLAC1); (P) Mestome sheath maker gene, LOC_Os01g68540 (GDI1); (Q–S) Representative images showing the results of in situ hybridization. [Black triangles- the identified cell types and gene ID; the scale bare is shown in images].
Functional enrichment of each cluster
The DEGs upregulated in each cluster were identified to obtain the cell types’ basic information in all clusters; 3,780 upregulated genes were screened (Supplementary Tables 7, 8). In all, 208 to 856 DEGs were identified. Procambium cells, followed by mesophyll and epidermal cells, had the highest number of DEGs; the least number of DEGs were found in fiber cells (Supplementary Table 7). Next, the potential enriched pathways and functions were determined by KEGG and GO analyses. GO analysis revealed that all eight cell-type clusters (except the phloem cluster) were significantly enriched in the “response to stimulus” biological process. The guard cell, mesophyll, and procambium clusters were significantly enriched in the “metabolic process.” Bulliform, epidermal, fiber, guard cell, mesophyll, and procambium clusters were significantly enriched in the “cellular process” (Figure 3A). The molecular function results showed bulliform and mestome sheath were significantly enriched in the “binding” terms. Epidermal and procambium clusters were significantly enriched in “catalytic activity” and “transporter activity” (Figure 3A).
Figure 3 Function enrichment analyses. (A) GO analysis of all clusters; (B) KEGG analysis of all clusters. The level 2 GO terms, such as biological process, molecular function, cellular components, and the top 5 Go pathways and terms, are shown. The number of enriched genes was represented as the size of the circle. The significant enrichments were presented in red, and the insignificant enrichments were presented in blue.
KEGG analysis revealed that the bulliform cluster was significantly enriched in “ribosome,” “protein processing in the endoplasmic reticulum,” and “phagosome” pathways. The epidermal cluster was significantly enriched in “endocytosis,” “glycerolipid-,” and “glycerophospholipid-” metabolism pathways. Fiber cluster was significantly enriched in “ubiquitin-mediated proteolysis,” “ribosome,” “oxidative phosphorylation,” “β-Alanine metabolism,” “fatty acid degradation,” “valine, leucine, and isoleucine degradation,” and “histidine metabolism” pathways. Guard cell cluster was significantly enriched in “plant-pathogen interaction” and “MAPK signaling” pathways. Mesophyll cluster was significantly enriched in “ribosome,” “photosynthesis,” “carbon fixation in photosynthetic organisms,” “glycolysis/gluconeogenesis,” “oxidative phosphorylation,” “carbon metabolism,” “mannose and fructose metabolism,” ‘photosynthesis-antenna proteins,” and “valine, leucine, and isoleucine degradation” pathways. Procambium cluster was significantly enriched in “citrate cycle (TCA cycle),” “phenylalanine metabolism,” “proteasome,” “tryptophan, tyrosine, and phenylalanine biosynthesis,” “α-linolenic acid metabolism,” “biosynthesis of amino acids,” “carbon metabolism,” “biosynthesis of secondary metabolites,” “glutathione metabolism,” “pentose phosphate pathway,” “oxidative phosphorylation,” and “stilbenoid, diarylheptanoid, and gingerol biosynthesis” pathways (Figure 3B).
The bulliform cluster was similar to the mestome sheath, and the epidermal cluster was similar to the procambium in the Go enrichment profiles.
Susceptible and resistant rice varieties mounted different responses to BPH feeding, which also differed based on cell type
The pseudotime trajectory analysis of all cell-type clusters was conducted to evaluate the responses of the rice cells to BPH infestation (Figure 4A). The development and response processes to BPH feeding could be divided into nine differentiation states (states 1–9) (Figure 4B). Pseudotime path clustering of DEGs revealed branching in the gene expression pattern of the differentiation states (Figure 4C, Supplementary Table 9). Susceptible and resistant rice varieties can mount different gene expression responses to BPH infestation. By comparing the differences in each cell type in TN1 and YHY15, we found that the procambium, epidermal, and fiber cells are similar across all differentiation states. Between TN1 and YHY15, significant differences in differentiation states were observed as follows–a) mestome sheath cells-states 1, 4, 6, and 8; b) guard cells-states 4, 6, 8, 9; c) mesophyll cells-state 8; d) xylem cells-states 1, 6, 8, 9; e) bulliform cells-states 8, 9; and f) phloem cells-states 1, 4, 6, 8, 9 (Figure 4D).
Figure 4 Pseudotime trajectory of all cluster’s cells. (A) Pseudotime analysis using Monocle for cell transcriptomes; (B) The state information of differentiation; (C) Gene expression heatmap for all cluster genes; (D) Expression profile in differentiation state of all subgroup cells.
Expression of genes related to vanillin, capsaicin, and ROS production in the mesophyll cluster depends on BPH-susceptibility
Although epidermal cells are the primary barrier of plants when BPH suck the phloem sap, their stylet bundle mainly pierces mesophyll cells; thus, mesophyll cells form a secondary barrier against further BPH feeding. The pseudotime trajectory analysis results of the mesophyll cells showed that the mesophyll cells were mainly divided into five states (Figure 5A). Comparing the gene expression in TN1 and YHY15 showed that the differences were primarily concentrated in branch 5 (Figure 5B). The overall expression level in cluster 5 was increased compared to that in other branches (Figure 5C). KEGG enrichment analysis of cluster 5 showed that these genes were significantly enriched in “phenylalanine, tyrosine, and tryptophan biosynthesis,” “biosynthesis of amino acid,” “MAPK signaling,” and “phenylalanine metabolism” pathways (Figure 5D, Supplementary Table 10). Phenylalanine metabolism is one of the downstream regulatory pathways of phenylalanine, tyrosine, and tryptophan biosynthesis. In this pathway, the phenylalanine ammonia-lyase (PAL) gene expression–viz. (LOC_Os02g41670, LOC_Os02g41680), CYP73A (LOC_Os02g26810, LOC_Os05g25640) and 4CL (LOC_Os08g34790)–were mainly affected. These are part of gene regulatory pathways that affect the 4-Coumaroyl-CoA levels and, consequently, the levels of vanillin or capsaicin. The GO enrichment analysis of cluster 5 genes showed that the main enriched GO terms include “response to stimulus” (GO: 0050896), “response to biotic stimulus” (GO: 0009607), etc. (Figure 5E). Three genes–two PAL homologs, two CYP73A homologs, and a 4CL homolog (Figures 5G–K)–in the phenylalanine and vanillin synthesis pathways were significantly different between the two rice varieties (Figure 5F). Further, cell number and expression of five genes in cluster 5 were quite different between TN1 and YHY15 varieties. Furthermore, the ethylene response pathway genes–RAN1, EBF1/2, and ChiB–under the MAPK signaling pathway had significant differences in expression levels (Figures 5L–P); thus, they may be involved in the trauma response of rice to BPH feeding. In addition, the expression of cluster 5 gene LOC_Os11g33120 (encoding a ROBH gene involved in reactive oxygen species (ROS) production) was significantly higher in TN1 than in YHY15. Thus, BPH feeding may have triggered differences in ROS production (Figure 5Q).
Figure 5 Pseudotime trajectory and function analysis of Mesophyll cells. (A) The state information of Mesophyll cells differentiation; (B) Samples information of Mesophyll cells; (C) Pseudotime trajectory of Mesophyll cells; (D) KEGG enrichment analysis of cluster 5 of Mesophyll cells; (E) GO enrichment analysis of cluster 5 of Mesophyll cells; (F) The schematic diagram of phenylalanine pathway; (G–K) The expression profiles of genes LOC_Os02g41670 (PAL) (G); LOC_Os02g41680 (PAL) (H); LOC_Os02g26810 (CYP73A) (I); LOC_Os05g25640 (CYP73A) (J), and LOC_Os08g34790 (4CL) (K) in TN1 and YHY15 samples; (L) The schematic diagram of ethylene response pathway; (M–Q) The expression profiles of genes LOC_Os06g45500 (RAN1) (M); LOC_Os02g10700 (EBF1/2) (N); LOC_Os06g40360 (EBF1/2) (O); LOC_Os06g51060 (ChiB) (P) and LOC_Os11g33120 (ChiB) (Q) in TN1 and YHY15 samples.
Phloem cell
Since BPH sucks phloem sap, the feedback of phloem cells to stress is inextricably linked to rice BPH resistance. The pseudotime trajectory analysis showed that the phloem cells were mainly divided into three clusters (Figure 6A). The gene expression level was lowest in cluster 1 and up-regulated in clusters 2 and 3 (Figure 6B). Cluster 1 mainly consists of TN1 cells, while clusters 2 and 3 have more YHY15 cells (Figure 6C). The gene expression level of cluster 1 gradually decreased along pseudotime, while clusters 2 and 3 showed an increasing trend (Figure 6D). The main GO enrichment terms for these genes include “response to stress,” “response to abiotic stimulus,” and “response to stimulus” (Figure 6E). Cluster 1 has multiple genes closely related to cell wall extension, such as LOC_Os06g48160 (XTH22, xyloglucan endotransglucosylase/hydrolase protein 22) (Figure 6F), LOC_Os08g40690 (RIXI, xylanase inhibitor protein 1) (Figure 6G). It also has energy production genes LOC_Os11g10480 (ADH1, alcohol dehydrogenase I) (Figure 6H). Multiple ABC transporter G family members in cluster 3–such as LOC_Os09g29660 (ABCG11) (Figure 6I), LOC_Os08g29570 (ABCG44) (Figure 6J), LOC_Os07g33780 (ABCG43) (Figure 6K), LOC_Os01g42410 (ABCG37) (Figure 6L), LOC_Os01g261460 (ESK1, promotes xylan acetylation) (Figure 6M), LOC_Os08g21040 (ASPG1, aspartic protease in guard cell) (Figure 6N)–exhibited significantly elevated expression.
Figure 6 Pseudotime trajectory and function analysis of Phloem cells. (A) Pseudotime trajectory of Phloem cells; (B) The state information of Phloem cells differentiation; (C) Samples information of Phloem cells; (D) Gene expression heatmap for all cluster genes of Phloem cells; (E) GO enrichment analysis of Phloem cells; (F–N) Representative graph showing the trend of the selected DEGs expression along pseudotime trajectory during differentiation for each cell type- (F) LOC_Os06g48160 (XTH22); (G) LOC_Os08g40690 (RIXI); (H) LOC_Os11g10480 (ADH1); (I) LOC_Os09g29660 (ABCG11); (J) LOC_Os08g29570 (ABCG44); (K) LOC_Os07g33780 (ABCG43); (L) LOC_Os01g42410 (ABCG37); (M) LOC_Os01g261460 (ESK1); (N) LOC_Os08g21040 (ASPG1). [One single cell was represented as one point. The entire X-axis was defined as “Pseudotime,” entire Y-axis was defined as “Relative expression”].
Xylem cells
BPH also sucks xylem sap (Seo et al., 2009). The xylem cells were divided into 3 clusters (Figure 7A). The overall expression level was lowest in cluster 1 and higher in clusters 2 and 3 (Figure 7B). Among the three clusters, the main difference between the two varieties was found in cluster 2 (Figures 7C, D). The gene expression level in cluster 1 gradually decreased; it was enriched in genes associated with GO terms “response to stress.” On the other hand, the gene expression level in clusters 2 and 3 showed a gradual increase; they were enriched in genes associated with GO terms, including “transport,” “localization,” “response to stress,” etc. Cluster 3 was enriched with genes associated with GO terms, such as “catalytic activity,” “secondary metabolic process,” “biosynthetic process,” and “response to stress” (Figure 7E). However, enrichment of KEGG pathways for the three clusters showed that cluster 1 was mainly enriched in pathways associated with “glycolysis/gluconeogenesis,” cluster 2 was enriched primarily in pathways such as “phenylalanine, tyrosine, and tryptophan biosynthesis,” including “amino sugar and nucleotide sugar metabolism,” especially chitin related genes, such as LOC_Os01g47070 (CHIT3), LOC_Os02g39330 (Cht6), LOC_Os04g52730 (UEL-2), LOC_Os05g29990 (UXS2), and so on. In addition, some pectin-related genes (such as LOC_Os02g29530 (GAUT8), LOC_Os02g51130 (GAUT9), etc.) showed significant changes (Figures 7F–K). The major KEGG pathways enriched in cluster 3 are “phenylpropanoid biosynthesis,” “phenylalanine metabolism,” and “phenylalanine, tyrosine, and tryptophan biosynthesis.”
Figure 7 Pseudotime trajectory and function analysis of Xylem cells. (A) The state information of Xylem cells differentiation; (B) Pseudotime trajectory of Xylem cells; (C) The expression pattern of Xylem cells in TN1; (D) The expression pattern of Xylem cells in YHY15; (E) Gene expression heatmap and GO enrichment analysis for all cluster genes of Xylem cells; (F–K) Representative graphs showing the profiles of the selected DEGs expression and trend- (F) LOC_Os01g47070 (CHIT3); (G) LOC_Os02g39330 (Cht6); (H) LOC_Os04g562730 (UEL-2); (I) LOC_Os05g29990 (UXS2); (J) LOC_Os02g29530 (GAUT8); (K) LOC_Os02g51130 (GAUT9). [One single cell was represented as one point. The entire X-axis was defined as “Pseudotime,” entire Y-axis was defined as “Relative expression”].
Discussion
BPH is one of the most damaging rice pests, causing substantial economic losses. Breeding BPH-resistant rice varieties has been a successful strategy, but BPH can eventually co-evolve and adapt to the resistance mechanisms of rice plants. BPH has a stylet bundle mouthpart that can pierce and suck the phloem sap of rice. However, when feeding on resistant rice varieties, BPH may fail to reach the phloem and stop feeding due to the presence of repellent substances in any cell type along the way. Therefore, analyzing the BPH-feeding stimulated cell expression patterns of different rice plant tissues will help identify the cells that mediate rice resistance to BPH and reveal the underlying molecular mechanisms. Through enabling transcriptomic analysis at single-cell resolution, the application of scRNA-seq has significantly revolutionized the study of cell and molecular biology. Moreover, it has dramatically enhanced our ability to characterize cell states and gene expression responses to BPH feeding. Here, we first constructed a single-cell atlas of rice leaf sheath response to BPH infestation.
Previous studies have used scRNA-seq technology to investigate tissue and organ development and differentiation in rice, thus, identifying marker genes that served as important references for our research (Xu et al., 2021b; Zhang et al., 2021a; Zong et al., 2022). However, reports on how different cell types respond–either by producing resistance chemicals or activating other molecular processes –to BPH feeding are lacking. Therefore, we selected 48 h post-infestation time for this study based on previous results that show significant changes in gene expression occurring between 24–48 hours of BPH infestation of rice leaves (Liu Y. et al., 2021; Xu et al., 2021a; Xue et al., 2023).
Here, we found that mainly mestome sheath cells, guard cells, mesophyll cells, xylem cells, bulliform cells, and phloem cells showed differential gene expression in resistant (YHY15) and susceptible (TN1) rice varieties after BPH infestation (Figure 4). This suggests that the factors imparting resistance to BPH may originate from multiple sources; different cells may contribute to insect resistance. Thus, the combined effects of numerous resistance factors may confer insect-resistance properties on plants. Previous studies using electrical penetration graphs and honeydew clocks have shown that rice resistance to BPH is determined by differences in sustained phloem ingestion, not by phloem location (Ghaffar et al., 2011). However, in the susceptible TN1 variety, BPH ingests phloem sap continuously without interruption (Ghaffar et al., 2011). Our results further support that rice resistance to BPH does not arise from a single resistance factor in the phloem but from a combination of resistance factors present in multiple locations in the rice plant.
Here, we found that the cell numbers were higher for mesophyll, procambium, and epidermal cells, with a significant difference in the proportion of mesophyll cells (cluster 4 in Figure 1, Supplementary Table 4) and procambium cells (cluster 2 in Figure 1, Supplementary Table 4). In contrast, the proportion of epidermal cells was less different in TN1 and YHY15 rice varieties (clusters 5, 7, and 8 in Figure 1, Supplementary Table 4). In addition, although the proportion of guard, mestome sheath, and bulliform cells differed approximately by 1% between TN1 and YHY15 rice varieties, the number of these cell types was relatively small. Moreover, procambium cells are the precursors of differentiated mature cells. Further, BPH feeding may mainly involve xylem and phloem sap ingestion (Sōgawa, 1982; Seo et al., 2009). Therefore, we focused on mesophyll, xylem, and phloem cells in this study.
Basal resistance is present in both susceptible and resistant rice varieties. The release of green leaf volatiles, which can prevent BPH infestation, is promoted by BPH feeding (Qi et al., 2011). Additionally, MAPKs, ethylene, and salicylic acid (SA) related signaling pathways were also activated by BPH feeding (Du et al., 2009; Hu et al., 2011; Lu et al., 2011). Here, the DEGs of mesophyll cells in TN1 and YHY15 varieties were mainly related to “phenylalanine, tyrosine, and tryptophan biosynthesis,” “MAPK signaling pathway,” and “phenylalanine metabolism.” PALs is a crucial enzyme that mediates the resistance to BPH by regulating the biosynthesis and accumulation of SA and lignin (He et al., 2020). We found that the downstream genes of the PAL pathway (CYP73A and 4CL), which may mainly induce the synthesis of vanillin and other compounds, were significantly differentially expressed in the two varieties. Vanillin-containing plants have vigorous insecticidal and insect-repellent activities (Kim et al., 2012; Songkro et al., 2012; Kletskova et al., 2017). MAPK signaling pathway is an essential part of the ethylene signaling transduction, which also plays a crucial role in plant defense response to insects (Hu et al., 2011; Hettenhausen et al., 2015; Zhou et al., 2019). The copper-transporting ATPase RAN1 is essential for the biogenesis of ethylene receptors (Binder et al., 2010). F-box proteins EBF1/EBF2 can form SCF complex to degrade EIN3 protein and regulate the expression level of downstream gene ChiB which is involved in the rice defense response to BPH (Zhu and Guo, 2008).
In phloem cells, the functions associated with DEGs in the two varieties mainly include cell wall extension, energy production, etc. However, cells of the susceptible variety TN1 were mainly concentrated in cluster 1, primarily associated with cell wall extension function. While clusters 2 and 3 had more cells of the BPH-resistant YHY15 variety, whose function was mainly related to “response to biotic stimulus” (Figure 6). Clusters 2 and 3 were also enriched in “phenylpropanoid biosynthesis pathway genes” such as LOC_Os02g41670 (PAL), LOC_Os05g35290 (PAL), LOC_Os08g34790 (4CL5), LOC_Os02g08100 (4CL3), LOC_Os09g04050 (CCR1), and LOC_Os01g73200 (PRDX6); these genes may promote lignin biosynthesis.
In xylem cells, we found that although DEGs were mainly enriched with genes significantly associated with stress response, their mainly enriched KEGG pathways were practically different, especially cluster 2 (Figure 7), which contained multiple genes related to chitin and pectin metabolism. In plants, chitin oligosaccharides induce various defense responses across multiple plant cells (Kaku et al., 2006). Pectin is a critical component of the cell wall. Therefore, pectin metabolism may be crucial in cell wall integrity and mediate plant defense responses (Wang et al., 2023).
In summary, the differences in immune responses stimulated by BPH infestation–such as MAPK signaling pathway and lignin biosynthesis for preventing callose and cell wall degradations that restrict BPH feeding and disrupt BPH digestion–are caused by the existing differences in TN1 and YHY15 varieties. Compared to the susceptible rice variety, the amplified and accelerated responses of the resistant rice variety may protect the plant from further BPH attack, allowing them to survive. The results of scRNA-seq suggest that multiple resistive factors may work together to make rice plants resistant to BPH, and different cell types may have other molecular mechanisms underlying BPH resistance.
In this study, we successfully—1) mapped a single cell transcriptome atlas of rice leaf sheath affected by BPH infestation; 2) observed the relationship of differentiation among the cell clusters; and 3) identified multiple cell types that may be involved in BPH-defense response mounted by rice. However, further experimental evidence is needed to elucidate the role of multiple cells and genes in BPH resistance. Nevertheless, our investigation provides a basis for mining insect-resistance genes and deciphering the molecular mechanism underlying insect-resistance and importantly guides insect-resistant plant molecular breeding.
Data availability statement
The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive in National Genomics Data Center (https://ngdc.cncb.ac.cn/), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences, under project number: PRJCA014345, accession number CRA009555 that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.
Author contributions
WZ, DX, LZ, and AY conceived and designed the experiments. WZ performed most of the experiments. CL, JC, and YW analyzed the data, authored or reviewed article drafts, and approved the final draft. MS analyzed the data. SL and BW provided help with the Q16 RNA in situ hybridization experiments. SS, KL, and HX prepared figures and/or tables. PL and KL helped to collect the samples. Finally, GY and ZC helped to revise the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National Key Research and Development Program of China (No.2021YFD1401100) and the National Natural Science Foundation of China (No.31501654).
Acknowledgments
We are grateful to Dr. Minshan Sun and Henan Assist Research Biotechnology Co., Ltd. (Zhengzhou, China) for assisting in bioinformatics analysis. And we would like to thank TopEdit for the English language editing of this manuscript. We show our appreciation to Guangzhou Gene Denovo Biotechnology Co., Ltd. for support with single-cell sequencing.
Conflict of interest
Author MS is employed by Henan Assist Research Biotechnology Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2023.1200014/full#supplementary-material
References
Akanksha, S., Lakshmi, V. J., Singh, A. K., Deepthi, Y., Ram, T. (2019). Genetics of novel brown planthopper Nilaparvata lugens (Stål) resistance genes in derived introgression lines from the interspecific cross O. sativa var. swarna ×. O. nivara. J. Genet. 98, 113. doi: 10.1007/s12041-019-1158-2
Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20, 163–172. doi: 10.1038/s41590-018-0276-y
Athwal, D. S., Pathak, M. D., Bacalangco, E. H., Pura, C. D. (1971). Genetics of resistance to brown planthoppers and green leafhoppers in oryza sativa l. Crop Sci. 11, 47–750. doi: 10.2135/cropsci1971.0011183X001100050043x
Binder, B. M., Rodríguez, F. I., Bleecker, A. B. (2010). The copper transporter RAN1 is essential for biogenesis of ethylene receptors in Arabidopsis. J. Biol. Chem. 285, 37263–37270. doi: 10.1074/jbc.M110.170027
Boyle, E. I., Weng, S., Gollub, J., Jin, H., Botstein, D., Cherry, J. M., et al. (2004). GO::TermFinder–open source software for accessing gene ontology information and finding significantly enriched gene ontology terms associated with a list of genes. Bioinformatics 20, 3710–3715. doi: 10.1093/bioinformatics/bth456
Brennecke, P., Anders, S., Kim, J. K., Kolodziejczyk, A. A., Zhang, X., Proserpio, V., et al. (2013). Accounting for technical noise in single-cell RNA-seq experiments. Nat. Methods 10, 1093–1095. doi: 10.1038/nmeth.2645
Butler, A., Hoffman, P., Smibert, P., Papalexi, E., Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36, 411–420. doi: 10.1038/nbt.4096
Cabauatan, P. Q., Cabunagan, R. C., Choi, I. R. (2009). Rice viruses transmitted by the brown planthopper Nilaparvata lugens stål (International Rice Research Institute: Los Baños, Philippines: Planthoppers: new threats to the sustainability to of intensive rice production systems in Asia), 357–368.
Camp, J. G., Sekine, K., Gerber, T., Loeffler-Wirth, H., Binder, H., Gac, M., et al. (2017). Multilineage communication regulates human liver bud development from pluripotency. Nature 546, 533–538. doi: 10.1038/nature22796
Chung, N. C., Storey, J. D. (2015). Statistical significance of variables driving systematic variation in high-dimensional data. Bioinformatics 31, 545–554. doi: 10.1093/bioinformatics/btu674
Du, B., Zhang, W., Liu, B., Hu, J., Wei, Z., Shi, Z., et al. (2009). Identification and characterization of Bph14, a gene conferring resistance to brown planthopper in rice. Proc. Natl. Acad. Sci. U.S.A. 106, 22163–22168. doi: 10.1073/pnas.0912139106
Ghaffar, M. B., Pritchard, J., Ford-Lloyd, B. (2011). Brown planthopper (N. lugens stal) feeding behaviour on rice germplasm as an indicator of resistance. PloS One 6, e22137. doi: 10.1371/journal.pone.0022137
Hanzelmann, S., Castelo, R., Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf. 14, 7. doi: 10.1186/1471-2105-14-7
He, J., Liu, Y., Yuan, D., Duan, M., Liu, Y., Shen, Z., et al. (2020). An R2R3 MYB transcription factor confers brown planthopper resistance by regulating the phenylalanine ammonia-lyase pathway in rice. Proc. Natl. Acad. Sci. 117, 271–277. doi: 10.1073/pnas.1902771116
Hettenhausen, C., Schuman, M. C., Wu, J. (2015). MAPK signaling: a key element in plant defense response to insects. Insect Sci. 22, 157–164. doi: 10.1111/1744-7917.12128
Hu, J., Zhou, J., Peng, X., Xu, H., Liu, C., Du, B., et al. (2011). The Bphi008a gene interacts with the ethylene pathway and transcriptionally regulates MAPK genes in the response of rice to brown planthopper feeding. Plant Physiol. 156, 856–872. doi: 10.1104/pp.111.174334
Islam, S., Zeisel, A., Joost, S., La Manno, G., Zajac, P., Kasper, M., et al. (2014). Quantitative single-cell RNA-seq with unique molecular identifiers. Nat. Methods 11, 163–166. doi: 10.1038/nmeth.2772
Jabnoune, M., Secco, D., Lecampion, C., Robaglia, C., Shu, Q. Y., Poirier, Y. (2015). An efficient procedure for protoplast isolation from mesophyll cells and nuclear fractionation in rice. BIO-PROTOCOL 5, e1412. doi: 10.21769/BioProtoc.1412
Jin, J., Tian, F., Yang, D. C., Meng, Y. Q., Kong, L., Luo, J., et al. (2017). PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res. 45, D1040–D1045. doi: 10.1093/nar/gkw982
Kaku, H., Nishizawa, Y., Ishii-Minami, N., Akimoto-Tomiyama, C., Dohmae, N., Takio, K., et al. (2006). Plant cells recognize chitin fragments for defense signaling through a plasma membrane receptor. Proc. Natl. Acad. Sci. 103, 11086–11091. doi: 10.1073/pnas.0508882103
Kanehisa, M., Araki, M., Goto, S., Hattori, M., Hirakawa, M., Itoh, M., et al. (2008). KEGG for linking genomes to life and the environment. Nucleic Acids Res. 36, D480–D484. doi: 10.1093/nar/gkm882
Kanehisa, M., Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28, 27–30. doi: 10.1093/nar/28.1.27
Kang, M., Choi, Y., Kim, H., Kim, S. G. (2022). Single-cell RNA-sequencing of Nicotiana attenuata corolla cells reveals the biosynthetic pathway of a floral scent. New Phytol. 234, 527–544. doi: 10.1111/nph.17992
Kawahara, Y., de la Bastide, M., Hamilton, J. P., Kanamori, H., Mccombie, W. R., Ouyang, S., et al. (2013). Improvement of the Oryza sativa nipponbare reference genome using next generation sequence and optical map data. Rice (N Y) 6, 4. doi: 10.1186/1939-8433-6-4
Kim, S. I., Yoon, J. S., Baeck, S. J., Lee, S. H., Ahn, Y. J., Kwon, H. W. (2012). Toxicity and synergic repellency of plant essential oil mixtures with vanillin against Aedes aegypti (Diptera: culicidae). J. Med. Entomol 49, 876–885. doi: 10.1603/me11127
Kletskova, A. V., Potkin, V. I., Dikusar, E. A., Zolotar, R. M. (2017). New data on vanillin-based isothiazolic insecticide synergists. Nat. Prod Commun. 12, 105–106. doi: 10.1177/1934578X1701200130
Korsunsky, I., Millard, N., Fan, J., Slowikowski, K., Zhang, F., Wei, K., et al. (2019). Fast, sensitive and accurate integration of single-cell data with harmony. Nat. Methods 16, 1289–1296. doi: 10.1038/s41592-019-0619-0
Laksminarayana, A., Khush, G. S. (1977). New genes for resistance to the brown planthopper in rice. Crop Sci. 17, 96–100. doi: 10.2135/cropsci1977.0011183X001700010028x
Li, Z., Xue, Y., Zhou, H., Li, Y., Usman, B., Jiao, X., et al. (2019). High-resolution mapping and breeding application of a novel brown planthopper resistance gene derived from wild rice (Oryza. rufipogon griff). Rice (N Y) 12, 41. doi: 10.1186/s12284-019-0289-7
Liberzon, A., Birger, C., Thorvaldsdottir, H., Ghandi, M., Mesirov, J. P., Tamayo, P. (2015). The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004
Liu, Q., Liang, Z., Feng, D., Jiang, S., Wang, Y., Du, Z., et al. (2021). Transcriptional landscape of rice roots at the single-cell resolution. Mol. Plant 14, 384–394. doi: 10.1016/j.molp.2020.12.014
Liu, Y., Wang, W., Li, Y., Liu, F., Han, W., Li, J. (2021). Transcriptomic and proteomic responses to brown plant hopper (Nilaparvata lugens) in cultivated and bt-transgenic rice (Oryza sativa) and wild rice (O. rufipogon). J. Proteomics 232, 104051. doi: 10.1016/j.jprot.2020.104051
Lopez-Anido, C. B., Vaten, A., Smoot, N. K., Sharma, N., Guo, V., Gong, Y., et al. (2021). Single-cell resolution of lineage trajectories in the Arabidopsis stomatal lineage and developing leaf. Dev. Cell 56, 1043–1055 e1044. doi: 10.1016/j.devcel.2021.03.014
Lu, J., Ju, H., Zhou, G., Zhu, C., Erb, M., Wang, X., et al. (2011). An EAR-motif-containing ERF transcription factor affects herbivore-induced signaling, defense and resistance in rice. Plant J. 68, 583–596. doi: 10.1111/j.1365-313X.2011.04709.x
Luecken, M. D., Theis, F. J. (2019). Current best practices in single-cell RNA-seq analysis: a tutorial. Mol. Syst. Biol. 15, e8746. doi: 10.15252/msb.20188746
Lun, A. T. L., Riesenfeld, S., Andrews, T., Dao, T. P., Gomes, T., Participants in the 1st Human Cell Atlas, J, et al. (2019). EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data. Genome Biol. 20, 63. doi: 10.1186/s13059-019-1662-y
Macosko, E. Z., Basu, A., Satija, R., Nemesh, J., Shekhar, K., Goldman, M., et al. (2015). Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell 161, 1202–1214. doi: 10.1016/j.cell.2015.05.002
McGinnis, C. S., Murrow, L. M., Gartner, Z. J. (2019). DoubletFinder: doublet detection in single-cell RNA sequencing data using artificial nearest neighbors. Cell Syst. 8, 329–337 e324. doi: 10.1016/j.cels.2019.03.003
Neftel, C., Laffy, J., Filbin, M. G., Hara, T., Shore, M. E., Rahme, G. J., et al. (2019). An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell 178, 835–849 e821. doi: 10.1016/j.cell.2019.06.024
Nguyen, C. D., Verdeprado, H., Zita, D., Sanada-Morimura, S., Matsumura, M., Virk, P. S., et al. (2019). The development and characterization of near-isogenic and pyramided lines carrying resistance genes to brown planthopper with the genetic background of japonica rice (Oryza sativa l.). Plants (Basel) 8, 498. doi: 10.3390/plants8110498
Qi, J., Zhou, G., Yang, L., Erb, M., Lu, Y., Sun, X., et al. (2011). The chloroplast-localized phospholipases d alpha4 and alpha5 regulate herbivore-induced direct and indirect defenses in rice. Plant Physiol. 157, 1987–1999. doi: 10.1104/pp.111.183749
Rotta, R., Noack, A. (2011). Multilevel local search algorithms for modularity clustering. ACM J. Exp. Algorithmics 16. doi: 10.1145/1963190.1970376
Senthil-Nathan, S., Choi, M. Y., Paik, C. H., Seo, H. Y., Kalaivani, K. (2009). Toxicity and physiological effects of neem pesticides applied to rice on the Nilaparvata lugens stal, the brown planthopper. Ecotoxicol Environ. Saf. 72, 1707–1713. doi: 10.1016/j.ecoenv.2009.04.024
Seo, B. Y., Kwon, Y.-H., Jung, J. K., Kim, G.-H. (2009). Electrical penetration graphic waveforms in relation to the actual positions of the stylet tips of Nilaparvata lugens in rice tissue. J. Asia-Pacific Entomology 12, 89–95. doi: 10.1016/j.aspen.2009.02.002
Sezer, M., Butlin, R. K. (1998). The genetic basis of host plant adaptation in the brown planthopper (Nilaparvata lugens). Heredity 80, 499–508. doi: 10.1046/j.1365-2540.1998.00316.x
Sōgawa, K. (1982). The rice brown planthopper: feeding physiology and host plant interactions. Annu. Rev. Entomology 27, 49–73. doi: 10.1146/annurev.en.27.010182.000405
Song, Q., Ando, A., Jiang, N., Ikeda, Y., Chen, Z. J. (2020). Single-cell RNA-seq analysis reveals ploidy-dependent and cell-specific transcriptome changes in Arabidopsis female gametophytes. Genome Biol. 21, 178. doi: 10.1186/s13059-020-02094-0
Songkro, S., Jenboonlap, M., Boonprasertpon, M., Maneenuan, D., Bouking, K., Kaewnopparat, N. (2012). Effects of glucam p-20, vanillin, and fixolide on mosquito repellency of citronella oil lotions. J. Med. Entomol 49, 672–677. doi: 10.1603/me11141
Spiller, N. J. (1990). An ultrastructural study of the stylet pathway of the brown planthopper Nilaparvata lugens. Entomologia Experimentalis Applicata 54, 191–193. doi: 10.1111/j.1570-7458.1990.tb01329.x
Stuart, T., Butler, A., Hoffman, P., Hafemeister, C., Papalexi, E., Mauck, W. M., 3rd, et al. (2019). Comprehensive integration of single-cell data. Cell 177, 1888–1902 e1821. doi: 10.1016/j.cell.2019.05.031
Umeda, M., Umeda-Hara, C., Yamaguchi, M., Hashimoto, J., Uchimiya, H. (1999). Differential expression of genes for cyclin-dependent protein kinases in rice plants. Plant Physiol. 119, 31–40. doi: 10.1104/pp.119.1.31
Wang, Y., Huan, Q., Li, K., Qian, W. (2021). Single-cell transcriptome atlas of the leaf and root of rice seedlings. J. Genet. Genomics 48, 881–898. doi: 10.1016/j.jgg.2021.06.001
Wang, D., Kanyuka, K., Papp-Rupar, M. (2023). Pectin: a critical component in cell-wall-mediated immunity. Trends Plant Sci. 28, 10–13. doi: 10.1016/j.tplants.2022.09.003
Xie, J., Li, M., Zeng, J., Li, X., Zhang, D. (2022). Single-cell RNA sequencing profiles of stem-differentiating xylem in poplar. Plant Biotechnol. J. 20, 417–419. doi: 10.1111/pbi.13763
Xu, X., Crow, M., Rice, B. R., Li, F., Harris, B., Liu, L., et al. (2021b). Single-cell RNA sequencing of developing maize ears facilitates functional analysis and trait candidate gene discovery. Dev. Cell 56, 557–568 e556. doi: 10.1016/j.devcel.2020.12.015
Xu, N., Wei, S. F., Xu, H. J. (2021a). Transcriptome analysis of the regulatory mechanism of FoxO on wing dimorphism in the brown planthopper, Nilaparvata lugens (Hemiptera: delphacidae). Insects 12, 413. doi: 10.3390/insects12050413
Xue, Y., Muhammad, S., Yang, J., Wang, X., Zhao, N., Qin, B., et al. (2023). Comparative transcriptome-wide identification and differential expression of genes and lncRNAs in rice near-isogenic line (KW-Bph36-NIL) in response to BPH feeding. Front. Plant Sci. 13. doi: 10.3389/fpls.2022.1095602
Xue, J., Zhou, X., Zhang, C. X., Yu, L. L., Fan, H. W., Wang, Z., et al. (2014). Genomes of the rice pest brown planthopper and its endosymbionts reveal complex complementary contributions for host adaptation. Genome Biol. 15, 521. doi: 10.1186/s13059-014-0521-0
Zhang, T. Q., Chen, Y., Liu, Y., Lin, W. H., Wang, J. W. (2021a). Single-cell transcriptome atlas and chromatin accessibility landscape reveal differentiation trajectories in the rice root. Nat. Commun. 12, 2053. doi: 10.1038/s41467-021-22352-4
Zhang, T. Q., Chen, Y., Wang, J. W. (2021b). A single-cell analysis of the Arabidopsis vegetative shoot apex. Dev. Cell 56, 1056–1074 e1058. doi: 10.1016/j.devcel.2021.02.021
Zhang, A. W., O’flanagan, C., Chavez, E. A., Lim, J. L. P., Ceglia, N., Mcpherson, A., et al. (2019). Probabilistic cell-type assignment of single-cell RNA-seq for tumor microenvironment profiling. Nat. Methods 16, 1007–1015. doi: 10.1038/s41592-019-0529-1
Zheng, X., Zhu, L., He, G. (2021). Genetic and molecular understanding of host rice resistance and Nilaparvata lugens adaptation. Curr. Opin. Insect Sci. 45, 14–20. doi: 10.1016/j.cois.2020.11.005
Zhou, S., Chen, M., Zhang, Y., Gao, Q., Noman, A., Wang, Q., et al. (2019). OsMKK3, a stress-responsive protein kinase, positively regulates rice resistance to Nilaparvata lugens via phytohormone dynamics. Int. J. Mol. Sci. 20, 3023. doi: 10.3390/ijms20123023
Zhu, Z., Guo, H. (2008). Genetic basis of ethylene perception and signal transduction in Arabidopsis. J. Integr. Plant Biol. 50, 808–815. doi: 10.1111/j.1744-7909.2008.00710.x
Zong, Y., Chen, Z., Innes, J. B., Chen, C., Wang, Z., Wang, H. (2007). Fire and flood management of coastal swamp enabled first rice paddy cultivation in east China. Nature 449, 459–462. doi: 10.1038/nature06135
Keywords: single-cell, RNA-seq, leaf sheath cells, resistant, Nilaparvata lugens, rice
Citation: Zha W, Li C, Wu Y, Chen J, Li S, Sun M, Wu B, Shi S, Liu K, Xu H, Li P, Liu K, Yang G, Chen Z, Xu D, Zhou L and You A (2023) Single-Cell RNA sequencing of leaf sheath cells reveals the mechanism of rice resistance to brown planthopper (Nilaparvata lugens). Front. Plant Sci. 14:1200014. doi: 10.3389/fpls.2023.1200014
Received: 04 April 2023; Accepted: 26 April 2023;
Published: 19 June 2023.
Edited by:
Shengli Jing, Xinyang Normal University, ChinaReviewed by:
Peiying Hao, China Jiliang University, ChinaHao Zhou, Sichuan Agricultural University, China
Copyright © 2023 Zha, Li, Wu, Chen, Li, Sun, Wu, Shi, Liu, Xu, Li, Liu, Yang, Chen, Xu, Zhou and You. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Deze Xu, ZGV6ZXh1QDE2My5jb20=; Lei Zhou, eXV0aWFuX3pob3U4M0AxNjMuY29t; Aiqing You, YXFfeW91QDE2My5jb20=
†These authors have contributed equally to this work