Role of long non-coding RNA in regulatory network response to Candidatus Liberibacter asiaticus in citrus

Abstract

Long non-coding RNAs (lncRNAs) serve as crucial regulators in plant response to various diseases, while none have been systematically identified and characterized in response to citrus Huanglongbing (HLB) caused by Candidatus Liberibacter asiaticus (CLas) bacteria. Here, we comprehensively investigated the transcriptional and regulatory dynamics of the lncRNAs in response to CLas. Samples were collected from leaf midribs of CLas- and mock-inoculated HLB-tolerant rough lemon (Citrus jambhiri) and HLB-sensitive sweet orange (C. sinensis) at week 0, 7, 17, and 34 following inoculation using CLas+ budwood of three biological replicates in the greenhouse. A total of 8,742 lncRNAs, including 2,529 novel lncRNAs, were identified from RNA-seq data with rRNA-removed from strand-specific libraries. Genomic variation analyses of conserved lncRNAs from 38 citrus accessions showed that 26 single nucleotide polymorphisms (SNPs) were significantly correlated with HLB. In addition, lncRNA-mRNA weighted gene co-expression network analysis (WGCNA) showed a significant module correlated with CLas-inoculation in rough lemon. Notably, the most significant LNC_28805 and multiple co-expressed genes related to plant defense in the module were targeted by miRNA5021, suggesting that LNC28805 might compete with endogenous miR5021 to maintain the homeostasis of immune gene expression levels. Candidate WRKY33 and SYP121 genes targeted by miRNA5021 were identified as two key hub genes interacting with bacteria pathogen response genes based on the prediction of protein-protein interaction (PPI) network. These two genes were also found within HLB-associated QTL in linkage group 6. Overall, our findings provide a reference for a better understanding of the role of lncRNAs involved in citrus HLB regulation.

1 Introduction

Transcripts with a length of more than 200 nt and lower protein-coding potential are operationally termed long non-coding RNAs (lncRNAs), which are widespread non-coding RNAs (ncRNAs) in eukaryotes. In animals and plants, lncRNA can function as important and versatile regulators in a variety of cellular and biological processes (Kim and Sung, 2012). Biochemical experiments and genetics studies have demonstrated that lncRNAs are associated with chromatin modification (Rinn and Chang, 2012), mRNA splicing (Bardou et al., 2014), transcriptional gene silencing (Wierzbicki, 2012), and posttranscriptional gene regulation (Yoon et al., 2013). Studies indicate that lncRNAs can help the host to prevent pathogen replication or be used by pathogens to promote pathogen proliferation (Li et al., 2016; Zaynab et al., 2018; Shirahama et al., 2020). In tomato, yellow leaf curl virus intergenic siRNAs target a host long noncoding RNA to modulate disease symptoms (Yang et al., 2019); also, tomato lncRNA23468 functions as a competing endogenous miR482b to enhance the defense against Phytophthora infestans (Jiang et al., 2019). It was also found that lncRNAs function as precursors of miRNAs having stable hairpin structures in wheat (Xin et al., 2011). In a word, lncRNAs are involved in pathogenic infection by acting as miRNA targets, miRNA precursors, or endogenous target mimics (eTMs) to regulate the expression of their target genes (Shirahama et al., 2020; Song et al., 2021).

Recent studies have revealed a set of the important regulatory functions of lncRNAs in response to pathogen infection. Transcriptome analyses revealed that a large number of lncRNAs were differentially expressed in response to pathogen infection in plants. For instance, lncRNA16397 was involved in resistance to P. infestans infection by co-expressing glutaredoxin in tomato (Cui et al., 2017). In melon, lncRNAs function as miRNA precursors and are involved in the response of powdery mildew fungi (Gao et al., 2020). Also, lncRNAs are involved in the response of Arabidopsis thaliana to Fusarium oxysporum infection (Zhu et al., 2014) and cotton to Verticillium dahlia infection (Zhang et al., 2018). In addition, a study also found that interaction between Sl-lncRNA15492 and Sl-miR482a can affect Solanum lycopersicum immunity against P. infestans (Jiang et al., 2020). Thus, lncRNAs are important components of regulated networks in response to pathogen infection. Although many lncRNAs have been identified from transcriptome data in diverse plant species, most of them are not well characterized.

Huanglongbing (HLB), a disease caused by the phloem-limited bacterium Candidatus Liberibacter asiaticus (CLas), is the most prevalent and destructive citrus epidemic. It has devastated the citrus industry in Florida and is threatening the global citrus industries (Bové, 2006; Graham et al., 2020). Thus far, citrus HLB has not been controlled effectively, and some research directions are precluded because of the challenge of CLas’ unculturable and phloem-limited nature. HLB and its vector, Asian citrus psyllid (ACP, Diaphorina citri) is still rapidly spreading in citrus-producing areas, which leads to billions of dollars in annual economic loss (Alvarez et al., 2016; Wang, 2019; Monzó and Stansly, 2020). Current strategies for insecticide and antibiotics application are limited and unsustainable. One of the most effective and eco-friendly strategies is strengthening host plant defense and immunity. Usually, the plant innate immune response can be triggered when they are infected by the pathogen (Nobori and Tsuda, 2019). In citrus, CLas-triggered plant immune responses are delayed 5–9 weeks after inoculation (Albrecht and Bowman, 2008; Yu et al., 2017). Traditional molecular biology, genetic, and multi-omics analyses also incompletely revealed the nature of pathogenesis of citrus HLB (da Graça et al., 2016; Wang, 2019). A study indicated that Citrus tristeza virus (CTV) can produce LMT1 lncRNA to suppress salicylic acid (SA) accumulation and mitigate reactive oxygen species (ROS) accumulation (Kang et al., 2019). These cases related to lncRNA involved in plant disease regulation bring us a promising direction to explore the role of lncRNA against citrus HLB disease (Jiang et al., 2020; Hong et al., 2022; Sharma et al., 2022). Thus far, no studies have explored such roles and characteristics in citrus.

Several citrus relatives such as Poncirus trifoliata and Microcitrus were considered as tolerant or resistant (Ramadugu et al., 2016; Godfrey et al., 2017), however, they are genetically distant from commercial citrus varieties which have originated mainly from three common ancestors, wild mandarin (C. reticulata), pummelo (C. maxima), and citron (C. medica) through a long domestication evolution (Wu et al., 2018). Rough lemon shares wild mandarin as a common ancestor with sweet orange. Our previous study showed that rough lemon is HLB tolerant compared with sweet orange (Fan et al., 2012). Once rough lemon trees are infected and symptomatic, they can be rejuvenated by the continued growth of new shoots with few or no foliar symptoms of the disease, and they repeat this cycle for many growing seasons; in contrast, sweet orange exhibits continuous growth inhibition and eventual dieback (Fan et al., 2012). By comparative transcriptional and anatomical analysis of rough lemon and sweet orange in response to CLas, phloem transport activity and the expression of defense-related genes are much greater in rough lemon than in sweet orange (Fan et al., 2012; Yu et al., 2017), suggesting the ability to maintain good phloem transport with extensive CLas titer is likely critical to good HLB tolerance. To further explore the contributions of lncRNA in response to HLB, we systematically identified lncRNAs from rough lemon and sweet orange at four different time points of a greenhouse experiment and characterized their genomic transcriptional and regulatory dynamics. We predicted their potential regulatory genes and functions and constructed a co-expression network. Our study provides valuable information and expands the knowledge of the role of citrus lncRNA in HLB disease expression.

2 Materials and methods

2.1 Plant materials

The plant inoculation was performed using the method as previously described (Fan et al., 2012). Briefly, two-year-old HLB-sensitive sweet orange (C. sinensis L Osb.) and HLB-tolerant rough lemon (C. jambhiri Lush.) plants were graft-inoculated with CLas positive bud wood collected from Carrizo citrange (C. sinensis × P. trifoliata L. Raf.) grown in a protected greenhouse, and mock-inoculated controls used bud wood from pathogen tested and healthy Carrizo trees. Each treatment had three biological replicates. All these plants were kept in a state-certified disease-free greenhouse (a United States Department of Agriculture Animal and Plant Health Inspection Service and Center for Disease Control-approved and secured greenhouse) at the University of Florida, Citrus Research and Education Center, Lake Alfred. Midribs of mature leaves were sampled from CLas-inoculated and mock-inoculated trees every two weeks after inoculation (WAI) at early stages (before ten weeks) and every one week at later stages; quantitative real-time PCR (qRT-PCR) was performed to test for the presence of CLas as previously described (Li et al., 2006). Plants were considered HLB positive when PCR cycle threshold (CT) values were below 30 (Yu et al., 2017). Positive plants were not detected until 17 WAI. The typical blotchy mottled HLB symptom was observed around 34 WAI in rough lemon and sweet orange. Based on the presence of positive samples and HLB symptoms, midribs of mature leaf from CLas-inoculated and mock-inoculated rough lemon and sweet orange trees at 0, 7, 17, and 34 WAI were collected. A total of 48 samples from four different time-points were used for RNA-seq in this study. The information about plant materials and the HLB test is shown in Table S1. Midribs from the mature leaves of rough lemon and sweet orange, which were under Huanglongbing (HLB) disease stress for more than ten years in the field, were also used for qPCR validation to explore whether the regulatory relationship between candidate lncRNAs and genes also existed in the plants from different growth condition.

2.2 Strand-specific RNA sequencing

Total RNA was extracted using TRIzol^® Reagent and purified using a RNeasy Mini Kit (Valencia, CA, United States) following the manufacturer’s protocol. Ribosomal RNA was removed from the total RNA using a Ribo-Zero rRNA removal kit (Epicenter-Illumina, Madison, WI, United States) following the manufacturer’s protocol. High-quality RNA was used to construct strand-specific RNA (ssRNA) libraries at the Interdisciplinary Center for Biotechnology Research (ICBR) Gene Expression Core, University of Florida (UF) described by Yu et al. (Yu et al., 2017). The prepared libraries were sequenced on an Illumina HiSeq 2000 (Illumina Inc., San Diego, CA, USA) producing paired-end 100 bp reads.

2.3 Transcript assembly and lncRNA identification

Low-quality reads and the adaptor sequences were removed using the fastp tool (Chen et al., 2018). After filtering, the clean reads were mapped to the Citrus clementina genome v1.0 (JGI) (Wu et al., 2014) using HISAT2 software (Kim et al., 2019). Next, StringTie was used to assemble transcripts of each sample, merge transcripts to get a consensus transcriptome assembly, and compute the abundance of these transcripts (Pertea et al., 2016). Subsequently, the newly assembled transcripts were compared with the C. clementina reference genome annotations using the GffCompare program (Pertea and Pertea, 2020). Transcripts overlapped with the known genes were discarded. The resulting transcripts with length ≥ 200 bp and fragments per kilobase of transcript per million fragments (FPKM) ≥ 1 in more than three samples were extracted, and then the tRNAs and rRNA were removed from the extracted transcripts using tRNAccan-SE (Lowe and Eddy, 1997) and RNAmmer (Lagesen et al., 2007), respectively. To reduce the noise of transcripts encoding proteins, TBtools software (Chen et al., 2020) was used to identify the open reading frames (ORFs) of these transcripts. Transcripts with significant ORFs were aligned to the Swiss-Prot, Nr, and Pfam databases, and the transcripts with E-value ≤ 10^-5 were excluded. Finally, we further evaluated the coding ability of the remaining transcripts using the Coding Potential Calculator version2 (CPC2) (Kang et al., 2017). We searched the lncRNAs in the PLncDB V2.0 database (Jin et al., 2021) and CANTATAdb (Szcześniak et al., 2016) to identify novel lncRNAs. The information on lncRNAs is listed in Table S2. Based on the genomic coordinates of protein-coding genes, the lncRNAs were divided into five groups (Scheuermann and Boyer, 2013; Ransohoff et al., 2018): Intergenic lncRNAs (LINC), intronic lncRNA (INTRONIC), natural antisense transcripts (NAT), genic lncRNA (GENIC), and exonic lncRNA (EXON).

2.4 Prediction of lncRNAs targets, precursors, and eTMs

Mature miRNA and miRNA precursors were downloaded from the miRbase database (Release 22.1) (Kozomara et al., 2019). The psRNATarget (Dai et al., 2018) was used to predict lncRNAs acting as putative miRNA targets with the default settings. All 8742 identified lncRNAs were used to predict eTMs using the psMimic software with the default settings (Wu et al., 2013). The lncRNAs transcripts were aligned to the miRNA precursors sequences from the miRbase database to predict the precursors using BLASTN software based on the best hits with E-value < 1e-10 and query identify > 80%. The miRNA target of candidate genes was predicted using TAPIR (http://bioinformatics.psb.ugent.be/webtools/tapir/) online tool.

2.5 Differential expression analysis of lncRNAs and mRNAs

The expression level of lncRNAs and mRNAs were quantified based on the position of mapped reads using StringTie software (Pertea et al., 2016) and evaluated by FPKM. The sample biological replicates were examined using principal component analysis (PCA) and correlation analysis, and the differentially expressed (DE) analysis was performed using the DEseq2 package in R software (Love et al., 2014). DE mRNAs and lncRNAs were determined with false discovery rate (FDR) < 0.01 and fold change (FC) |log2FC| ≥ 1. Protein coding (PC) genes were annotated using Swiss-Prot (Boeckmann et al., 2003) and non-redundant (NR) database (Pruitt et al., 2005). Gene functional enrichment was analyzed using Metascape (Zhou et al., 2019) and MapMan (Thimm et al., 2004).

2.6 Orthologous identification and phylogenetic analysis of lncRNAs

The orthologous lncRNAs were identified based on reciprocal best blast hits (RHB) using the Basic Local Alignment Search Tool (BLAST) with query coverage > 80% and E-value < 1e-10 (Moreno-Hagelsieb and Latimer, 2008). The orthologous lncRNA sequences of Arabidopsis thaliana, Oryza sativa, and Populus trichocarpa were downloaded from the PLncDB V2.0 database (Jin et al., 2021). A total of 36 orthologues were identified. For phylogenetic analysis, sequences of 38 accessions (Table S3) were from previous publications except for Citrus latipes (Wu et al., 2014; Wu et al., 2018). Variant calling and filtering were according to the method previously described (Peng et al., 2020). Briefly, sequences were mapped to the Citrus clementina genome v1.0 (JGI) (Wu et al., 2014) via BWA-MEM (Li, 2013). Raw aligned reads were sorted and duplicate reads removed via samtools V1.7 (Li et al., 2009) and sambamba V0.6.7 (Tarasov et al., 2015), respectively. Variant calling and filtering were performed using The Genome Analysis Toolkit (GATK) v4.1.2 (McKenna et al., 2010) and VCFtools (Danecek et al., 2011), respectively. Finally, 1,658 bi-allelic variants derived from 36 conserved orthologous lncRNAs genomic loci were extracted and used to construct a maximum-likelihood (ML) phylogenetic tree (Data S1). The best substitution model general time-reversible (GTR) for the ML tree was inferred using Smart Model Selection (SMS) web server (Lefort et al., 2017).

2.7 Construction of co-expression network for lncRNAs and mRNAs

First, we excluded the transcripts that had similar expression patterns between mock-inoculated and CLas-inoculated plants to reduce the noise caused by gene spatiotemporal-specific expression using the ‘Mfuzz’ package in R software with the k-nearest neighbor method (Peterson, 2009). Co-expression analysis was performed using the weighted co-expression network analysis (WGCNA) package in R software with threshold power = 6, minimum module size = 3, and a branch merge cut height = 0.25 (Langfelder and Horvath, 2008). The co-expression network was plotted using Gephi software (Bastian et al., 2009). Proteins of Arabidopsis thaliana were used as model to infer the protein-protein interaction (PPI) of co-expressed genes in citrus, based on STRING database (Szklarczyk et al., 2015).

2.8 cDNA synthesis and qRT-PCR

Twelve genes with high amplification efficiency and primer specificity were selected to validate the RNA-seq data of rough lemon and sweet orange using qRT-PCR. In addition, we also validate LNC28805, WRKY33, and SYP121 and their targeting miRNA5021 in rough lemon and sweet orange, which were under Huanglongbing (HLB) disease stress for more than ten years in the field in Lake Alfred. The STEM-LOOP RT-qPCR method was carried out for miRNA5021 based on the method of Kramer (2011) (Kramer, 2011). The qPCR method for mRNA and lncRNA was described by Yu et al. (2017). Briefly, the first strand cDNA was synthesized using an Affinityscript qPCR cDNA Synthesis Kit (Agilent Technologies, Santa Clara, CA, USA), and RT-qPCR was performed using SYBR Green qPCR Master Mix (Agilent Technologies) in a 20-μl volume. 18S rRNA gene was used for an internal reference according to previous studies (Yan et al., 2012; Wei et al., 2021). The primers are listed in Table S4.

3 Results

3.1 Identification and characterization of lncRNAs in citrus

Approximately 2,523 million paired-end reads from 48 sample libraries were produced and mapped to the C. clementina v1.0 reference genome (Table S5). After a comprehensive pipeline of filtering, a total of 8,742 lncRNAs with FPKM > 1 in at least three samples were identified in both sweet orange and rough lemon, including 2,529 novel lncRNAs (Figure 1A and Table S2). Most of lncRNAs belong to LINC and EXON groups (Figure 1C and Table S2). The EXON group has a greater exon number and transcript length in comparison to other groups (Figure S1A, B).

Figure 1

To investigate the characteristics of citrus lncRNAs response to HLB, we analyzed the correlation coefficients between different samples using the expression profiles of lncRNAs. A relatively high correlation between biological replicates was observed (Figure S2). Analysis of lncRNAs distribution in the citrus genome showed that lncRNAs were widely expressed across all citrus chromosomes, and the highest number of expressed lncRNA was around the 40 Mb region on scaffold_3 (Figure 1B). We then investigated the transcript length and exon number distribution of lncRNAs and mRNAs. The exon number of lncRNA and mRNA had similar distribution patterns in rough lemon and sweet orange (Figure S3), while mRNA had a greater exon number and sequence length than lncRNA (Figures 1D, E; Figure S3). Meanwhile, we also investigated the difference in the overall gene expression levels of lncRNA and mRNA in four different periods. The results show that expression levels of mRNA were always higher than lncRNAs in both mock- and CLas-inoculated plants (Figure S1C-D).

3.2 Differential expression dynamics of lncRNAs and mRNA after CLas-inoculation

The repeatability of samples was evaluated using PCA analysis. The VH1_W0 sample was eliminated due to poor replication (Figure S4). After that, twelve genes were selected to further validate the RNA-seq data, indicating good reliability of the data (Figure S5). Based on the analysis, we systematically compared the expression levels of lncRNAs and mRNAs in rough lemon and sweet orange plants at different time points of mock- and CLas-inoculated conditions. A total of 1,943 and 25,118 differentially expressed (DE) lncRNA and mRNAs were identified from 14 pairwise comparisons of rough lemon and sweet oranges (Data S2-S3). The percentage of DE lncRNAs and mRNAs was similar in both mock-inoculated sweet orange and rough lemon plants, but there was a prominent difference between them in comparisons of the CLas-inoculated groups (Figure 2A). A larger percentage of DE lncRNAs and mRNAs was found at week 7 and week 34 in CLas compared to groups of sweet orange, whereas there was no prominent difference in rough lemon (Figure 2A). These results indicate CLas altered the spatiotemporal expression pattern of plants and reflect the significant difference of genes in response to Clas infection between sweet orange and rough lemon. We further investigated the number of up- and down-regulated DE lncRNA in sweet orange and rough lemon after Clas-inoculation (Figures 2B, C). Comparing W0 and different time points after CLas inoculation, we found that the number of DE lncRNAs significantly increased at week 7 (W7) and week 34 (W34) in sweet orange (Figure 2B), and the largest number of DE lncRNA was at W7 (down/up, 103/176). However, the largest number of DE lncRNA was found at W34 in rough lemon, and only included 75 down-regulated and 77 up-regulated DE lncRNAs (Figure 2B). Comparing healthy plants with CLas infected plants at the same time point, we found that the number of upregulated DE lncRNAs significantly increased at W34 in rough lemon (down/up, 16/76), while it was the opposite in sweet orange at the same points (down/up, 36/9) (Figure 2C).

Figure 2

Heatmaps showed completely different expression patterns of lncRNAs between rough lemon and sweet orange species but a similar expression pattern in different time points of mock- and CLas-inoculated plants of the same species (Figure 2D), indicating lncRNAs have higher species specificity in response to Clas infection. In addition, we further investigated the dynamic changes of specific DE lncRNAs and mRNAs in rough lemon and sweet orange (Figure S6). The number of specific DE lncRNAs in the pairwise groups of mock-inoculation vs. CLas-inoculation at the same time point was higher in rough lemon than in sweet orange at W17 and W34 groups (Figure S6E). However, in the pairwise groups of W0 -W7, W0-W17, and W0-W34 of CLas-inoculation, we found the number of specific DE lncRNAs of sweet orange was near two-fold higher than rough lemon (Figure S6A), and a mass of specific DE lncRNAs were found at week 7 and week 34 after CLas-inoculation, and most of them were upregulated in both rough lemon and sweet orange. However, a larger number of DE lncRNAs were downregulated in the same pairwise mock inoculation groups (Figure S6B). In addition, we also found the specific DE lncRNAs were overall upregulated in rough lemon or sweet orange in the other compared groups (Figure S6C, D).

3.3 Evolution and function of specific DE lncRNAs in citrus

To further explore the characteristics of specific DE lncRNAs, we analyzed their genomic origins from rough lemon and sweet orange. The result shows that a total of 35 conserved lncRNAs, 166 specific lncRNAs in rough lemon, and 170 specific lncRNAs in sweet orange, were identified, and the distribution landscape of these specific DE lncRNAs was different between rough lemon and sweet orange on chromosome-level scaffolds (Figure 3A). Compared with the W0 time point (W7-W0, W17-W0, and W34-W0 pairwise groups of CLas-inoculated plants), most of DE lncRNAs were distributed in scaffold_1, 2, 3, 5, and their numbers exhibited significant differences between rough lemon and sweet orange (Figure 3B). The largest number was located on scaffold_5 and scaffold_3 in rough lemon and sweet orange, respectively (Figure 3B). In the pairwise groups between mock- and CLas-inoculated plants at the same time point, the largest percentage of DE lncRNAs (31.6%) was found in unanchored scaffolds in rough lemon, while the largest percentage was found on scaffold_3 in sweet orange (Figure 3C). We further investigated the number of common and unique genomic origin DE lncRNAs in rough lemon and sweet orange (Figure 3D). Only 43 specific DE lncRNAs derived from the same genomic loci of these two citrus species were identified, indicating that lncRNAs have high species-specific expression profiles in response to CLas.

Figure 3

Analysis of the evolutionary conservation of lncRNAs showed that thirty-one sequences of conserved lncRNAs located in seven different scaffolds were identified with e-value > 1×10⁻¹⁰ (Figure 3A and Table S6), and the expression levels of them represented specifically spatiotemporal different expression patterns between mock- and CLas-inoculated plants (Figures 3E, F). For instance, three lncRNAs (LNC_57342, LNC_20398, and LNC_61191) were significantly upregulated at W7 stage in rough lemon, while they were downregulated in sweet orange; the other three lncRNAs (LNC_61191, LNC_155540, and LNC_52993) exhibited reverse expression patterns in sweet orange (Figure 3E). Compared with W0, we found that seven lncRNAs exhibited higher differentially expressed levels in W7 and W17 time points after CLas inoculation in rough lemon and sweet orange. which might be involved in response to HLB in early phases after CLas inoculation. We further identified 1,658 bi-allelic variants derived from these homologous lncRNAs based on resequencing data of 38 citrus accessions, and a phylogenetic tree was constructed (Table S3 and Data S1). The result showed that citrus species with close relatives were clustered together with higher support values (Figure 3G). Based on the HLB symptom evaluation of citrus relatives (Folimonova et al., 2009; Ramadugu et al., 2016; Godfrey et al., 2017), correlation analysis between bi-allelic variants and HLB symptom evaluation showed that 26 variants were significantly correlated with the HLB traits (Figure 3H). Most of the significant SNPs were located in 41.64-41.65 Mb region on scaffold_3 and 28.45-29.95 Mb region on scaffold_9. Six conserved DE lncRNAs with these correlated SNPs derived from four genomic loci (XLOC_007316, XLOC_030383, XLOC_072976, and XLOC_037329) were identified (Table S7).

In this study, we also predicted the miRNA targets, precursors, and eTMs of lncRNAs (Data S4). A total of 133 lncRNAs were identified as precursors for 33 miRNA families, 116 lncRNAs were identified to be targeted by 35 miRNA families, and 40 lncRNAs were predicted as eTMs for 14 miRNA families. Notably, 16 lncRNAs predicted as miRNA targets were simultaneously acting as miRNA precursors or eTMs, and nine lncRNAs were differentially expressed in sweet orange (Table S8), suggesting that these DE lncRNAs may be involved in response to HLB by interacting with miRNAs in sweet orange.

3.4 Identification of lncRNA-mRNA co-expression modules related to HLB response

To identify the DE lncRNAs and mRNAs potential response to HLB and to reduce the noise caused by gene spatiotemporal-specific expression, we first excluded the DE lncRNAs or mRNAs that had similar expression patterns in different time points of mock- and CLas-inoculation plants by using the k-nearest neighbor method (Figure 4A and Figure S7, 8). For instance, 53 and 49 similar expression pattern lncRNAs were excluded in the comparison of cluster1 (CLas) vs. cluster7 (Mock) and cluster3 (CLas) vs. cluster3 (Mock) in rough lemon (Figure 4A), respectively. Finally, a total of 2133 (including 246 lncRNAs and 1887 mRNAs) and 2863 (including 295 lncRNAs and 2568 mRNAs) transcripts from rough lemon and sweet orange, respectively, were potentially responsible for HLB response and were used for co-expression analysis. This analysis resulted in 12 and 7 distinct modules, which are clusters of highly interconnected genes (Langfelder and Horvath, 2008) in rough lemon (Figure 4B) and sweet orange (Figure S9A), respectively. Module eigengene is considered a representative gene expression profile in a module and correlated with the corresponding tissue type or trait (Langfelder and Horvath, 2008). We found that 3 out of 12 co-expression modules are comprised of genes that are significantly correlated with Mock- or CLas-inoculation, and 7 out of 12 modules are significantly correlated with time points (0, 7, 17, 34 WAI) (P ≤ 0.05; Figure 4C) in rough lemon. However, no significant modules were correlated with CLas inoculation in sweet orange except for three significant modules correlated with the WAI (Figure S9B).

Figure 4

Modules with high trait significance may represent potential pathways associated with the trait, and genes with high module membership in modules significantly correlated with traits are hub candidate genes (Langfelder and Horvath, 2008). In this study, scatterplots showed that genes with high module membership in 8 modules in rough lemon (Figure 4D and Figure S10) and three modules in sweet orange (Figure S9C) were identified as potential hub genes with high significance (P < 0.01), including a total of 130 lncRNAs and 61 lncRNAs in rough lemon and sweet orange, respectively (Table S9). Among these lncRNAs, 23 lncRNAs were predicted to function as miRNA targets, miRNA precursors, or eTM (Table S10).

A heatmap showing the relative normalized FPKM (NFPKM) of lncRNAs revealed that most of the co-expression lncRNAs were from brown module correlated with CLas-inoculation and significantly upregulated at 7 WAI in rough lemon, and lncRNAs in different modules exhibited highly specific expression (Figure 4E). However, few co-expression lncRNAs showed spatiotemporal-specific expression in sweet orange (Figure 4F). Heatmap of co-expressed mRNAs also showed prominently different expression patterns between rough lemon and sweet orange (Figure S11). Genes in the pink and blue modules were specifically upregulated at 17 WAI and 34 WAI after CLas-inoculation in rough lemon, respectively (Figure S11A), while genes in the blue module of sweet orange were significantly downregulated at 7 WAI and 34 WAI (Figure S11B). These specific DE lncRNAs and mRNAs might specifically contribute to the response to CLas infection.

3.5 Functional annotation and enrichment of genes in significant modules

A total of 1146 and 879 were included in the eight significant modules in rough lemon and three significant modules in sweet orange, respectively. GO enrichment indicated that stress response-related biological process terms (such as response to reactive oxygen species, plant-pathogen interaction, and response to stimulus processes) were significantly over-represented in rough lemon (Figure 5A and Figure S12A). However, the development and growth of biological processes, negatively regulated cell proliferation, and response to stimulus terms were significantly over-represented in sweet orange (Figure 5B and Figure S12B). MapMan functional categories related to biotic stress show that redox state, glutathione-S-transferase, and secondary metabolites were significantly enriched in rough lemon (Figure S12C). Notably, genes enriched in the Glutathione-S-transferase category were mainly involved in pathogenic effector-triggered immunity, systemic acquired resistance, and WRKY33-dependent plant immunity in rough lemon, and most of these genes were significantly upregulated (Data S5 and Figure S12C). In sweet orange, most genes were mainly enriched in the redox state and secondary metabolites categories related to cell wall organization, but few genes were enriched in the Glutathione-S-transferase category related to pathogen response (Data S5 and Figure S12D). This result indicates pathogenic response genes in the Glutathione-S-transferase category play an important role in HLB tolerance of rough lemon. In addition, we also found that the functional roles of genes in the redox state and secondary metabolite categories were quite different between rough lemon and sweet orange. Compared with rough lemon, many cell wall related genes were enriched in sweet orange, including genes related to arabinogalactan, callose, and lipid biosynthesis (Data S5). According to the previous study, callose deposition can play a role as a defensive fortification in response to CLas bacteria in citrus (Achor et al., 2010). Based on the results of gene function enrichment, we suggest that the mechanisms of HLB response are quite different between rough lemon and sweet orange.

Figure 5

3.6 Co-expression regulatory network of lncRNA-mRNA in response to CLas infection

Based on the results of GO and MapMan functional enrichment, co-expression genes associated with pathogen stimuli response in rough lemon and genes associated with callose synthase in sweet orange were of particular interest (Table S11). Multiple genes (such as NLR, RIN4, and CBP60/SARD) were suggested to play an important role in plant innate immunity (Abramovitch et al., 2006; da Graça et al., 2016). Genes showing the most connections with those immunities might be mutually involved in the pathway of pathogen response.

Therefore, we further identified genes first connected to stimuli response and callose synthase related genes and constructed their potential co-expression regulatory network. A total of 205 and 73 first neighboring genes with weight > 0.2 were identified in rough lemon and sweet orange, respectively (Figure S13A, B and Data S6). Multiple defense-related genes were identified in the network, such as BOTRYTIS-INDUCED KINASE 1 (BIK1)/MSTRG.37943.1 and NECROTIC SPOTTED LESIONS 1 (NLS1)/MSTRG.17571.7 in rough lemon; peroxidase superfamily protein (PRX52)/MSTRG.57748.1 and subtilisin-like protease (SBT1.5)/MSTRG.58400.1 in sweet orange. GO enrichment showed that the neighboring genes in rough lemon and sweet orange were significantly enriched in plant-type hypersensitive response and cutin biosynthetic process, respectively (Figure S13C-D). In addition, seven pathogen-related genes derived from MapMan pathogen stimuli response category showed a strong association with each other (Figure S13A). Among these seven genes, effector-triggered immunity related RPM1-interacting factor 4 (RIN4)/MSTRG.29058.1 and WRKY33-dependent plant immunity related SIGMA FACTOR BINGD PROTEIN (SIB)/Ciclev10002803m.v1.0 with the highest edge weight are suggested to be two key hub genes in the co-expression network.

To know more about the protein-protein interaction (PPI) relationships between these connected genes, proteins of A. thaliana were used as model to infer the potential interaction network (Szklarczyk et al., 2015). A total of 43 and 13 PPI genes (edge confidence > 0.4) were identified in rough lemon and sweet oranges, respectively (Figure S14 and Data S7). Most of them in the PPI network were related to stress response (Figure 6A; Figure S14). For instance, WRKY33 (Ciclev10011386m.v1.0) is important for plant pathogen immunity and interacted with multiple disease related proteins in rough lemon. ACS6 is involved in bacterial flagellin-induced ethylene production (Gravino et al., 2015), SYP121 is involved in ABA signaling (Eisenach et al., 2012), NECROTIC SPOTTED LESIONS (NSL1) negatively regulated cell death programs, and salicylic acid-related defense responses; nsl1 mutants exhibited higher levels of salicylic acid (SA) and callose deposition (Noutoshi et al., 2006). As to the first neighboring genes of callose synthase in sweet orange, LONG-CHAIN ACYL-COA SYNTHETASE 2 (LACS2) is required for the accumulation of cuticular wax to enhance plant stress resistance (Figure S14B), which strongly interacts with CYP86A2 and is involved in the processes of cuticle development and repression of bacterial type III gene (Xiao et al., 2004). Notably, WRKY33 and its interacted SIB (Ciclev10002803m.v1.0) had the highest confidence in the PPI and co-expression networks (Figure 6A and Figure S13A), implying an extent regulatory relationship between them in response to CLas infection. Based on the PPI network, genes with high confidence (edge confidence > 0.7) were identified, and a set of them related to pathogenic response, such as effector receptor, systemic acquired resistance protein, callose synthase, and calcium-dependent protein et al., are identified as key genes in response to CLas infection (Table 1). The analysis of gene expression showed that pathogen response genes were distinctly upregulated in rough lemon at 17 weeks after CLas inoculation, while callose synthase related CALS1 was not expressed in rough lemon (Figure 6E). Interestingly, effector-triggered immunity related leucine-rich repeat receptor-like protein kinase family protein (AT4G08850), pathogen-associated molecular pattern (PAMP)-triggered immunity (PTI) signaling related CPK7, and ethylene signaling response gene ERF-1 were not expressed in sweet orange, which partially explained greater HLB tolerance.

Figure 6

Next, we further investigated the co-expression relationship between pathogen stimuli response genes and lncRNAs. We found three lncRNAs (LNC_28805, LNC_24384, and LNC_128248) that were co-expressed with pathogen response genes in rough lemon (Figure 6B). In this network, LNC_28805 was a hub lncRNA first connecting to all five pathogen stimuli response genes, suggesting it might be a key regulator involved in response to CLas. LncRNAs could regulate gene transcription levels and cellular processes by acting as miRNA precursors or miRNA mimics (Chekanova, 2015). Among the three lncRNAs, LNC28805 was solely found to act as miR167 and miR5021 targets in both rough lemon and sweet orange and exhibited similar expression patterns after CLas inoculation (Table S10 and Figure 6C). LNC28805 was significantly upregulated at 17 WAI in rough lemon. To realize whether other pathogenic response genes could be co-expressed with LNC28805, we further identified the first neighboring genes of LNC28805 with edge weight > 0.2 from rough lemon and sweet orange. A total of 145 differentially expressed transcripts first connect to LNC28805 in rough lemon (Data S8). However, none of the genes co-expressed with LNC28805 were found in the co-expression network of sweet orange. Functional enrichment analysis of the 145 transcripts showed that plant defense-related genes were the most significantly enriched (Figure 6D).

3.7 A hypothetical regulatory pathway of LNC28805 involved in HLB regulation

LNC28805 was not predicted to act as an eTM or miRNA precursor, but it did appear to act as a target for miR167 and miR5021, thus suggesting a role in the regulation of HLB response by interacting with miRNAs. To validate this point, we analyzed the miR5021 and miR167 potential targets among the co-expression genes of LNC28805, as well as pathogen stimuli response genes and callose synthase genes. A total of 25 transcripts were identified to be targeted by miR5021 in rough lemon and sweet orange, respectively (Table S12 and Data S9). However, none of the transcripts were identified to be targeted by miR167. Most of these genes were involved in disease regulation, such as WRKY33, SYP121, and NB-ARC domain-containing disease resistance protein (Table S12). Interestingly, we found eight genes targeted by miR5021 were found with the QTLs identified in our previous study (Huang et al., 2018) (Table S13). Intriguingly, QTLs on scaffold_6 were simultaneously detected for foliar symptoms (FS) and canopy damage (CD) in two different years, and there two key hub genes (WRKY33 and SYP121) in the PPI network, strongly linked to the QTL peak markers, were found (Figure 7A). This result further supports a putatively important regulatory relationship between LNC28805 and these disease response genes. It is tempting to speculate that LNC28805 might act as competing endogenous miR5021 to regulate HLB response genes by attenuating their cleavage or translated inhibition caused by miR5021. The mechanism of miRNA-LncRNA interactions regulating host immunity-related genes has been reported in tomato (Jiang et al., 2020). Overall, a hypothetical model of LNC28805 and its potential role in regulatory processes was proposed (Figure 7B). This model showed that WRKY33 might regulate cross-talk between jasmonate-, abscisic acid (ABA)-, and salicylic acid (SA)-regulated disease response pathways (Zheng et al., 2006; Birkenbihl et al., 2012; Liu et al., 2015).

Figure 7

Birkenbihl et al., 2012; Liu et al., 2015). LNC28805 might importantly compete with endogenous miR5021 to maintain the homeostasis of these immune gene expression levels. Red asterisks indicate miR5021 targets. (C) Relative expression of miRNA5021 and its targets in rough lemon and sweet orange under Huanglongbing (HLB) stress more than ten years in the field (shown in Figure S15). Asterisks indicate significant differences (Students’ t-test, *p < 0.05, **p < 0.01, ***p < 0.001).

To explore if the regulatory relationship also exists in the plants under HLB stress for a long term in the field, we further tested the expression levels of miR5021 and its targets (WRKY33, SYP121, and LNC28805) using midribs of rough lemon and sweet orange under HLB stress more than ten years in the field (Figure S15). We found LNC28805 was significantly upregulated, and miR5021 was significantly downregulated in rough lemon, which supports the hypothesis that LNC28805 is involved in regulating the expression level of pathogenic response genes by competing for endogenous miR5021.

4 Discussion

It has been known that lncRNAs act as versatile regulators involved in diverse biological processes in plants. In citrus, thousands of lncRNAs from different species have been identified in flowers, leaves, seeds, and fruits based on RNA-seq transcriptome data (Xu et al., 2018). However, lncRNAs have not been systematically identified in response to citrus HLB, and little is known about their biological function during the response to the disease. Here, we present a comprehensive picture of lncRNAs and mRNAs transcription patterns in different time points after CLas-inoculation in HLB-tolerant rough lemon and HLB-sensitive sweet orange. A set of DE lncRNAs exhibited stage-specific and species-specific expression, implying specific roles for lncRNAs in response to infection of citrus with CLas, leading to HLB.

4.1 Identification, conservation, and phylogenetic relationship of lncRNAs

In this study, we have identified 8742 high confident lncRNAs from mock- and CLas-inoculated rough lemon and sweet orange. There are more than 50% lncRNAs from intergenic regions, which is consistent with a study of lncRNAs and flowering in trifoliate orange (Wang et al., 2017). We also identified 2529 novel lncRNAs by aligning lncRNA sequences to the annotated lncRNA database (Szcześniak et al., 2016; Jin et al., 2021). These lncRNAs might specifically respond to HLB disease and need more investigation. Therefore, our study has expanded the information on disease-related lncRNAs in citrus. Compared to protein-coding genes, lncRNAs identified from rough lemon and sweet orange shared most of the common features of lncRNAs reported in other citrus species and model plants (Di et al., 2014; Yuan et al., 2018; Ke et al., 2019), such as low conservation, fewer exons, shorter transcript length, and lower expression levels. Similar features in different species might be explained by the high rates of origin sequence evolution (Mercer et al., 2009). These results are consistent with previous studies in other species (Ke et al., 2019; Tang et al., 2021), also supporting the reliability of lncRNAs identified in this study. Analysis of the primary sequence evolutionary conservation by comparing with A. thaliana, O. sativa, and P. trichocarpa lncRNAs showed that only 33 lncRNAs derived from 22 genomic loci were identified. Compared with Arabidopsis and rice, a relatively high proportion of conserved lncRNAs was found in P. trichocarpa, indicating high sequence similarity in a more closely related taxon. The phylogenetic tree of 38 citrus relatives also displayed high conservation of these lncRNAs genomic sequence at the intra-species level (Figure 3G), suggesting that these conserved lncRNAs might share a common ancestor and ancient evolutionary origin at genomic levels (Ulitsky, 2016). Similar results were also observed in animals and other plants (Pang et al., 2006; Di et al., 2014). Intriguingly, correlation analysis between SNPs derived from the conserved lncRNAs indicated 26 SNPs were significantly associated with HLB tolerance, and six conserved lncRNAs contained these SNPs were significantly upregulated at W7 time point after CLas-inoculation (Figure 3E and Table S7). Some studies reported that conserved lncRNAs might play an important role in conserved functions, such as interacting with RNA-binding proteins by conserving secondary structures (Wang et al., 2014; Ulitsky, 2016). Although experimental evidence about the biogenesis and functions of conserved lncRNAs is limited, most current studies suggest specific secondary structures of conserved lncRNAs might be important for a conserved function (Johnsson et al., 2014; Ulitsky, 2016). More extensive experiments are needed to validate whether these conserved lncRNAs have important contributions to HLB response at an early stage in our future studies.

4.2 Tissue-specific and species-specific expression pattern of lncRNAs

It is difficult to deduce and validate the lncRNA function due to its low conservation and non-coding property. The tissue-specific expression may help to understand the potential functions underlying lncRNAs. In the current study, lncRNAs exhibited more prominent species-specific expression patterns compared with spatiotemporal-specific expression patterns (Figure 2D). Pairwise comparisons between mock- or CLas-inoculated tissues or different time points revealed that the number of DE lncRNAs and mRNAs distinctly increased at week 7 and decreased at week 17 after CLas inoculation in sweet orange (Figure 2A). However, it did not change significantly in rough lemon. The greater numbers of DE lncRNAs and mRNAs in sweet orange may imply greater sensitivity and activity in response to external stimuli at transcriptional regulation levels. Changes in the environment and rhythmic plant growth also can affect plant gene expression levels (López-Maury et al., 2008; Covington et al., 2008). The pairwise comparison of lncRNAs between mock- and CLas-inoculated plants was further performed to reduce the environmental noise, presenting obvious differences of lncRNA number at 34 WAI between rough lemon and sweet orange (79 upregulated lncRNAs in rough lemon and nine upregulated lncRNAs in sweet orange). Moreover, most of these lncRNAs presented high species-specificity. These results indicated that species-specific DE lncRNA might play an essential role in the regulation of HLB tolerance in rough lemon. High tissue- and species-specific expression of lncRNAs, such as LNC_57342 and LNC_61191 specifically expressed in rough lemon and sweet orange at W7 time point, respectively, also make them potentially useful as biomarkers for HLB detection or screening HLB tolerant species in early stages. Based on the aforementioned results, we suggest these specific DE lncRNAs might contribute to citrus HLB tolerance or sensitivity and imply important dynamics of lncRNAs response to HLB in citrus.

4.3 Putative interaction between DE lncRNAs and miRNAs

The lncRNA-miRNA interaction is a vital regulatory mechanism of lncRNAs. Several studies have shown that lncRNAs can serve as miRNA and siRNA precursors to assist in target gene cleavage or translation inhibition, function as traps for miRNA binding, or are directly targeted by miRNA to attenuate miRNA presence in plant immunity (Song et al., 2021). For instance, lncRNA MuLnc1 was identified to be cleaved by mul-miR3954 to produce si16157 and negatively regulated Botrytis cinerea and Pseudomonas syringae resistance by inhibiting the functions of calmodulin-like protein 27 (CML27) in mulberry (Gai et al., 2018). In tomato, lncRNA42705 and lncRNA08711 increase the expression levels of MYB genes by acting as decoys for miR159 and enhance resistance against P. infestans (Cui et al., 2019); and, Sl-lncRNA15492 was targeted by SI-miR482a to maintain Sl-NBS-LRR1 at an appropriate expression level during the immune response to P. infestans (Jiang et al., 2020). In this study of citrus and HLB, we identified a set of high confidence lncRNAs that potentially serve as precursors, eTMs, and targets of miRNAs. Among them, it is noteworthy that five lncRNAs simultaneously serve as eTMs and targets of miRNA, and eleven lncRNAs simultaneously serve as precursors and miRNA targets (Table S8). Notably, nine DE lncRNAs were specifically identified in sweet orange. Several miRNAs that interacted with these lncRNAs have been reported to be involved in disease resistance, such as miR858 negatively regulated Arabidopsis immunity (Camargo-Ramírez et al., 2018), miR477 enhanced the susceptibility of the tea plants to Pseudopestalotiopsis species infection (Wang et al., 2020), and miRNA482 suppressed the expression of NBS-LRR defense genes in cotton (Zhu et al., 2014). Interestingly, LNC_40405 and LNC_69103 were predicted to act as eTMs of miR477 and miR482, respectively, and were significantly downregulated in CLas-inoculated sweet orange at 7 WAI and 17 WAI. In addition, multiple DE lncRNAs were also targeted by miR2111 and miR5021. A recent study indicated that miR2111 positively regulates shoot-to-root systemic effectors of rhizobia and promotes nodule formation (Moreau et al., 2021). Surprisingly, miR5021 was also predicted to be one having maximal matches against 19910 ESTs from periwinkle (Catharanthus roseus) (Pani and Mahapatra, 2013), which is an alternate host of the HLB bacteria (Zhang et al., 2010). Most of the predicted miR5021 targets were the genes involved in cell growth and development, signaling, and metabolism in periwinkle. In citrus, we also found multiple genes associated with disease response were targeted by miR5021 (Table S12). LNC_69103 was not only acting as an eTM of miR482 but also acting as a target of miR5021. Compared to healthy plants, LNC_69103 was extremely downregulated (log2FC < -9) in CLas-infected sweet orange plants. These findings indicated the expression level of these lncRNAs might be tightly related to the sensitivity of citrus plants to HLB.

4.4 A potential regulatory model of LNC_28805 in response to CLas

One of the main objectives of this study was to understand the expression dynamics and co-expression networks responding to citrus HLB in tolerant rough lemon and sensitive sweet orange. The WGCNA co-expression network revealed inoculation and temporal specific modules and important hub genes involved in citrus HLB response. Brown module was a unique significant module associated with CLas-inoculation in rough lemon and included a large group of co-expressed lncRNAs. Cluster analysis of expression pattern showed that lncRNAs displayed more temporal specific expression patterns than mRNAs and were significantly upregulated at 7 WAI (Figure 4E), suggesting that the functions of these lncRNAs might be closely related to HLB response at early stages after CLas-inoculation. We speculate that HLB-tolerant genes in rough lemon might be mediated by these early response lncRNAs. Evidence from most disease response genes exhibiting strong protein-protein interactions was identified in the blue module, but not the brown module, and these genes were significantly upregulated at 17 WAI, though not at 7 WAI (Figure S11 and Figure 6E). Functional enrichment analysis indicated that co-expressed genes in rough lemon were mainly responsible for the response to reactive oxygen species, programmed cell death, and plant pathogen, which are highly related to disease defense (Figure 6E). Although genes in response to bacteria were also enriched in sweet orange, most of them were responsible for growth and development, fatty acids biosynthesis, and reproductive processes (Figure 5B). A fraction of genes was even associated with negative regulation of response to stimulus in sweet orange. These results also reflect the greater sensitivity of sweet orange to HLB than rough lemon.

Additionally, MapMan analysis identified multiple genes involved in pathogenic effector-triggered immunity, systemic acquired resistance, WRKY33-dependent plant immunity, and callose synthase. Their first neighboring genes in the co-expression network were highly correlated with pathogen response and exhibited strong interactions with each other. Interestingly, WRKY33 and SYP121 were two key hub genes with the highest edge weight in the PPI network (Figure 6A), and both of them are located in significant QTLs identified by the previous study (Huang et al., 2018), suggesting that these two genes might play important roles in HLB regulation. In addition, we found a hub LNC_28805 was the most probable lncRNA involved in HLB regulation, which was co-expressed with multiple genes associated with the effector-triggered immunity (ETI) network, including RIN4, CBP60, SIB, and NLR. Moreover, LNC_28805 has 145 first neighboring DE mRNAs (including WRKY33 and SYP121) in the co-expression network. These neighboring genes were the most significantly enriched in the plant defense component (Figure 6D). Meanwhile, we also found that LNC_28805, belonging to intergenic lncRNAs type, was predicted to be targeted by the homologous miR5021 of Arabidopsis with one mismatch, indicating that LNC_28805 may be an evolutionarily conserved and functionally maintained lncRNA. We notably found that multiple disease resistance genes were targeted by miR5021, such as WRKY33, SYP121, and NB-ARC domain-containing disease resistance genes (Table S12), suggesting LNC28805 might compete with endogenous miR5021 to maintain the homeostasis of expression levels between immune-related genes and growth genes. All these results further indicate that LNC28805 might be an important regulator involved in the HLB tolerance of rough lemon.

Plant WRKY transcription factors (TF) play key roles in plant responses to microbial infection. The PPI network showed WRKY33 was involved in multiple disease response pathways (Figure 6A and Figure S14). WRKY33 can regulate the expression of defense-related genes toward the necrotrophic fungus B. cinerea, but WRKY33 has also been shown to regulate cross-talk between jasmonate-, abscisic acid (ABA)-, and salicylic acid (SA)-regulated disease response pathways (Zheng et al., 2006; Birkenbihl et al., 2012; Liu et al., 2015). In Arabidopsis, the ectopic expression of WRKY33 results in enhanced susceptibility to the bacterial pathogen P. syringae caused by the reduced expression of the salicylate-regulated PR-1 gene (Zheng et al., 2006). Loss of WRKY33 function results in activation of the ABA- and salicylic acid (SA)-related host response (Birkenbihl et al., 2012; Liu et al., 2015). The evidence indicated that high expression WRKY33 might suppress the expression of ABA- or SA-regulated genes response to HLB bacteria (Figure 7B). WRKY33 also can activate the expression of RING-type ubiquitin ligase ATL31 involved in vesicle trafficking with PEN1/SYP121 SNARE protein (Reyes et al., 2015), which function to guard cell membrane transport and stomatal control (Eisenach et al., 2012). Some others candidate genes in the PPI networks, such as ACS6, ERF-1, CML37, CALS1, NSL1, and MPK1 involved in bacteria-induced ethylene, ABA, and JA signaling, callose synthesis, and SA-related defense (Tsuda and Katagiri, 2010; Luna et al., 2010; Fukunaga et al., 2017), were also putatively involved the defense response of CLas invasion (Figures 6A and 7B). Based on the regulatory relationship of the genes in the networks, we suggested that LNC_28805 may play an important role in maintaining the homeostasis of antagonistic relationship between defense pathways mediating WRKY33 associated with ABA- or SA-regulated genes involved CLas response (Zheng et al., 2006; Liu et al., 2015).

In addition, WRKY33 also can activate the expression of Arabidopsis RING-type ubiquitin ligase ATL31 involved in vesicle trafficking with PEN1/SYP121 SNARE protein (Reyes et al., 2015), which functions to guard cell membrane transport and stomatal control (Eisenach et al., 2012). Because LNC28805 and some of its co-expressed genes were targeted by miRNA5021, we suggest that LNC28805 is probably involved in regulating the expression level of pathogenic response genes by competing for endogenous miR5021. If this hypothetical mechanism for HLB tolerance is correct, it should also be observed in plants grown in the field. The expression levels of miR5021 targets (WRKY33, SYP121, LNC28805) presented significantly higher in rough lemon than in sweet orange under HLB stress for more than ten years in the field (Figure 7C). However, the expression levels of miRNA5021 were reversed between them. It indicates that this relationship also exists in the naturally infected plants grown in the field. Thus, we suggest that miR5021 targeting WRKY33 and SYP121 might promote the expression of genes responding to CLas. Dynamic expression of WRKY33 might be required to balance the expression levels between immune-related genes and growth genes. LNC28805 probably plays an important role in regulating the expression level of these pathogenic response genes by competing for endogenous miR5021.

Though callose deposition plays a role in defense against the pathogen, overaccumulation of callose inhibits phloem transport activities in Clas-infected citrus (Achor et al., 2010). According to the previous study, callose-plugged phloem sieve elements were less serious in HLB-diseased rough lemon than in HLB-diseased sweet orange (Fan et al., 2012). A potential possibility of HLB tolerance mechanisms of rough lemon might be that LNC28805 is involved in competing for endogenous miR5021 to promote the expression of WRKY33 and SYP121, which might function to suppress the immune-related genes overresponse to CLas infection and enhance the activity of phloem transport by reducing callose deposition (Fan et al., 2012; Deng et al., 2019). Taken together, our results not only represent the gene modules of lncRNAs and mRNAs related to pathogenic response but also bring new insights into the roles of lncRNAs acting as potential regulatory factors for citrus HLB tolerance.

5 Conclusion

To conclude, we systematically identified and characterized 8,742 lncRNAs among HLB-tolerant rough lemon and HLB-sensitive sweet orange from different time points after CLas-inoculation. Based on the integrated analysis of sequence conservation and variation, spatiotemporal-specific expression, functional enrichment, and lncRNA-mRNA co-expression networks with WGCNA, we identified a fraction of lncRNAs and mRNAs that were potentially responsive to CLas bacterium infection in citrus. LNC_28805 was identified as one of the most important candidate lncRNAs involved in citrus HLB regulation. Two key candidates (WRKY33 and SYP121) in the PPI network are known to negatively regulate bacteria pathogen responses and were found within overlapping QTLs identified in our previous study. Based on the reported studies and PPI network and gene co-expression networks in this study, a putative hypothesis for the regulatory pathway of LNC_28805 is proposed (Figure 7B). This study will be useful in understanding the role of lncRNAs involved in citrus HLB regulation and provide a foundation for further investigation of their regulatory functions.

References

• 1

  AbramovitchR. B.AndersonJ. C.MartinG. B. (2006). Bacterial elicitation and evasion of plant innate immunity. Nat. Rev. Mol. Cell Biol.7 (8), 601–611. doi: 10.1038/nrm1984

  • CrossRef
  • Google Scholar

• 2

  AchorD. S.EtxeberriaE.WangN.FolimonovaS. Y.ChungK. R.AlbrigoL. G. (2010). Sequence of anatomical symptom observations in citrus affected with huanglongbing disease. Plant Pathol. J.9 (2), 56–64. doi: 10.3923/ppj.2010.56.64

  • CrossRef
  • Google Scholar

• 3

  AlbrechtU.BowmanK. D. (2008). Gene expression in citrus sinensis (L.) osbeck following infection with the bacterial pathogen candidatus liberibacter asiaticus causing huanglongbing in Florida. Plant Sci.175 (3), 291–306. doi: 10.1016/j.plantsci.2008.05.001

  • CrossRef
  • Google Scholar

• 4

  AlvarezS.RohrigE.SolísD.ThomasM. H. (2016). Citrus greening disease (Huanglongbing) in Florida: Economic impact, management and the potential for biological control. Agric. Res.5 (2), 109–118. doi: 10.1007/s40003-016-0204-z

  • CrossRef
  • Google Scholar

• 5

  BardouF.ArielF.Simpson CraigG.Romero-BarriosN.LaporteP.BalzergueS.et al. (2014). Long noncoding RNA modulates alternative splicing regulators in arabidopsis. Dev. Cell30 (2), 166–176. doi: 10.1016/j.devcel.2014.06.017

  • CrossRef
  • Google Scholar

• 6

  BastianM.HeymannS.JacomyM. (Eds.) (2009). Gephi: an open source software for exploring and manipulating networks. Proceedings of the international AAAI conference on web and social media. 3(1), 361–36. doi: 10.1609/icwsm.v3i1.13937

  • CrossRef
  • Google Scholar

• 7

  BirkenbihlR. P.DiezelC.SomssichI. E. (2012). Arabidopsis WRKY33 is a key transcriptional regulator of hormonal and metabolic responses toward Botrytis cinerea infection. Plant Physiol.159 (1), 266–285. doi: 10.1104/pp.111.192641

  • CrossRef
  • Google Scholar

• 8

  BoeckmannB.BairochA.ApweilerR.BlatterM.-C.EstreicherA.GasteigerE.et al. (2003). The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res.31 (1), 365–370. doi: 10.1093/nar/gkg095

  • CrossRef
  • Google Scholar

• 9

  BovéJ. M. (2006). Huanglongbing: A destructive, newly-emerging, century-old disease of citrus. J. Plant Pathology.88 (1), 7–37. Available at: http://www.jstor.org/stable/41998278

  • Google Scholar

• 10

  Camargo-RamírezR.Val-TorregrosaB.San SegundoB. (2018). MiR858-mediated regulation of flavonoid-specific MYB transcription factor genes controls resistance to pathogen infection in arabidopsis. Plant Cell Physiol.59 (1), 190–204. doi: 10.1093/pcp/pcx175

  • CrossRef
  • Google Scholar

• 11

  ChekanovaJ. A. (2015). Long non-coding RNAs and their functions in plants. Curr. Opin. Plant Biol.27, 207–216. doi: 10.1016/j.pbi.2015.08.003

  • CrossRef
  • Google Scholar

• 12

  ChenC.ChenH.ZhangY.ThomasH. R.FrankM. H.HeY.et al. (2020). TBtools: An integrative toolkit developed for interactive analyses of big biological data. Mol. Plant13 (8), 1194–1202. doi: 10.1016/j.molp.2020.06.009

  • CrossRef
  • Google Scholar

• 13

  ChenS.ZhouY.ChenY.GuJ. (2018). Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics34 (17), i884–ii90. doi: 10.1093/bioinformatics/bty560

  • CrossRef
  • Google Scholar

• 14

  CovingtonM. F.MaloofJ. N.StraumeM.KayS. A.HarmerS. L. (2008). Global transcriptome analysis reveals circadian regulation of key pathways in plant growth and development. Genome Biol.9 (8), R130. doi: 10.1186/gb-2008-9-8-r130

  • CrossRef
  • Google Scholar

• 15

  CuiJ.JiangN.HouX.WuS.ZhangQ.MengJ.et al. (2019). Genome-wide identification of lncRNAs and znalysis of ceRNA networks during tomato resistance to phytophthora infestans. Phytopathology®110 (2), 456–464. doi: 10.1094/PHYTO-04-19-0137-R

  • CrossRef
  • Google Scholar

• 16

  CuiJ.LuanY.JiangN.BaoH.MengJ. (2017). Comparative transcriptome analysis between resistant and susceptible tomato allows the identification of lncRNA16397 conferring resistance to Phytophthora infestans by co-expressing glutaredoxin. Plant J.89 (3), 577–589. doi: 10.1111/tpj.13408

  • CrossRef
  • Google Scholar

• 17

  da GraçaJ. V.DouhanG. W.HalbertS. E.KeremaneM. L.LeeR. F.VidalakisG.et al. (2016). Huanglongbing: An overview of a complex pathosystem ravaging the world's citrus. J. Integr. Plant Biol.58 (4), 373–387. doi: 10.1111/jipb.12437

  • CrossRef
  • Google Scholar

• 18

  DaiX.ZhuangZ.ZhaoP. X. (2018). psRNATarget: A plant small RNA target analysis server (2017 release). Nucleic Acids Res.46 (W1), W49–W54. doi: 10.1093/nar/gky316

  • CrossRef
  • Google Scholar

• 19

  DanecekP.AutonA.AbecasisG.AlbersC. A.BanksE.DePristoM. A.et al. (2011). The variant call format and VCFtools. Bioinformatics.27 (15), 2156–2158. doi: 10.1093/bioinformatics/btr330

  • CrossRef
  • Google Scholar

• 20

  DengH.AchorD.ExteberriaE.YuQ.DuD.StantonD.et al. (2019). Phloem regeneration is a mechanism for huanglongbing-tolerance of “Bearss” lemon and “LB8-9” sugar belle mandarin. Front. Plant Science.10, 277. doi: 10.3389/fpls.2019.00277

  • CrossRef
  • Google Scholar

• 21

  DiC.YuanJ.WuY.LiJ.LinH.HuL.et al. (2014). Characterization of stress-responsive lncRNAs in Arabidopsis thaliana by integrating expression, epigenetic and structural features. Plant J.80 (5), 848–861. doi: 10.1111/tpj.12679

  • CrossRef
  • Google Scholar

• 22

  EisenachC.ChenZ.-H.GrefenC.BlattM. R. (2012). The trafficking protein SYP121 of Arabidopsis connects programmed stomatal closure and k+ channel activity with vegetative growth. Plant J.69 (2), 241–251. doi: 10.1111/j.1365-313X.2011.04786.x

  • CrossRef
  • Google Scholar

• 23

  FanJ.ChenC.YuQ.KhalafA.AchorD. S.BrlanskyR. H.et al. (2012). Comparative transcriptional and anatomical analyses of tolerant rough lemon and susceptible sweet orange in response to ‘Candidatus liberibacter asiaticus’ infection. Mol. Plant-Microbe Interactions®.25 (11), 1396–1407. doi: 10.1094/MPMI-06-12-0150-R

  • CrossRef
  • Google Scholar

• 24

  FolimonovaS. Y.RobertsonC. J.GarnseyS. M.GowdaS.DawsonW. O. (2009). Examination of the responses of different genotypes of citrus to huanglongbing (citrus greening) under different conditions. Phytopathology®.99 (12), 1346–1354. doi: 10.1094/PHYTO-99-12-1346

  • CrossRef
  • Google Scholar

• 25

  FukunagaS.SogameM.HataM.Singkaravanit-OgawaS.Piślewska-BednarekM.Onozawa-KomoriM.et al. (2017). Dysfunction of arabidopsis MACPF domain protein activates programmed cell death via tryptophan metabolism in MAMP-triggered immunity. Plant J.89 (2), 381–393. doi: 10.1111/tpj.13391

  • CrossRef
  • Google Scholar

• 26

  GaiY.-P.YuanS.-S.ZhaoY.-N.ZhaoH.-N.ZhangH.-L.JiX.-L. (2018). A novel lncRNA, muLnc1, associated with environmental stress in mulberry (Morus multicaulis). Front. Plant Science.9, 669. doi: 10.3389/fpls.2018.00669

  • CrossRef
  • Google Scholar

• 27

  GaoC.SunJ.DongY.WangC.XiaoS.MoL.et al. (2020). Comparative transcriptome analysis uncovers regulatory roles of long non-coding RNAs involved in resistance to powdery mildew in melon. BMC Genomics21 (1), 125. doi: 10.1186/s12864-020-6546-8

  • CrossRef
  • Google Scholar

• 28

  GodfreyP. M.EdS.ChandrikaR.ManjunathL. K.RichardF. L. (2017). Apparent tolerance to huanglongbing in citrus and citrus-related germplasm. HortScience horts.52 (1), 31–39. doi: 10.1111/jipb.12437

  • CrossRef
  • Google Scholar

• 29

  GrahamJ.GottwaldT.SetamouM. (2020). Status of huanglongbing (HLB) outbreaks in Florida, California and Texas. Trop. Plant Pathology.45 (3), 265–278. doi: 10.1007/s40858-020-00335-y

  • CrossRef
  • Google Scholar

• 30

  GravinoM.SavatinD. V.MaconeA.De LorenzoG. (2015). Ethylene production in Botrytis cinerea- and oligogalacturonide-induced immunity requires calcium-dependent protein kinases. Plant J.84 (6), 1073–1086. doi: 10.1111/tpj.13057

  • CrossRef
  • Google Scholar

• 31

  HongY.ZhangY.CuiJ.MengJ.ChenY.ZhangC.et al. (2022). The lncRNA39896–miR166b–HDZs module affects tomato resistance to phytophthora infestans. J. Integr. Plant Biol.64 (10), 1979–1993. doi: 10.1111/jipb.13339

  • CrossRef
  • Google Scholar

• 32

  HuangM.RooseM. L.YuQ.DuD.YuY.ZhangY.et al. (2018). Construction of high-density genetic maps and detection of QTLs associated with huanglongbing tolerance in citrus. Front. Plant Science.9, 1694. doi: 10.3389/fpls.2018.01694

  • CrossRef
  • Google Scholar

• 33

  JiangN.CuiJ.HouX.YangG.XiaoY.HanL.et al. (2020). Sl-lncRNA15492 interacts with sl-miR482a and affects solanum lycopersicum immunity against phytophthora infestans. Plant J.103 (4), 1561–1574. doi: 10.1111/tpj.14847

  • CrossRef
  • Google Scholar

• 34

  JiangN.CuiJ.ShiY.YangG.ZhouX.HouX.et al. (2019). Tomato lncRNA23468 functions as a competing endogenous RNA to modulate NBS-LRR genes by decoying miR482b in the tomato-phytophthora infestans interaction. Horticulture Res.6, 28. doi: 10.1038/s41438-018-0096-0

  • CrossRef
  • Google Scholar

• 35

  JinJ.LuP.XuY.LiZ.YuS.LiuJ.et al. (2021). PLncDB V2.0: a comprehensive encyclopedia of plant long noncoding RNAs. Nucleic Acids Res.49 (D1), D1489–D1D95. doi: 10.1093/nar/gkaa910

  • CrossRef
  • Google Scholar

• 36

  JohnssonP.LipovichL.GrandérD.MorrisK. V. (2014). Evolutionary conservation of long non-coding RNAs; sequence, structure, function. Biochim. Biophys. Acta (BBA) - Gen. Subjects.1840 (3), 1063–1071. doi: 10.1016/j.bbagen.2013.10.035

  • CrossRef
  • Google Scholar

• 37

  KangS.-H.SunY.-D.AtallahO. O.Huguet-TapiaJ. C.NobleJ. D.FolimonovaS. Y. (2019). A long non-coding RNA of Citrus tristeza virus: Role in the virus interplay with the host immunity. Viruses.11 (5), 436. doi: 10.3390/v11050436

  • CrossRef
  • Google Scholar

• 38

  KangY.-J.YangD.-C.KongL.HouM.MengY.-Q.WeiL.et al. (2017). CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res.45 (W1), W12–WW6. doi: 10.1093/nar/gkx428

  • CrossRef
  • Google Scholar

• 39

  KeL.ZhouZ.XuX.-W.WangX.LiuY.XuY.et al. (2019). Evolutionary dynamics of lincRNA transcription in nine citrus species. Plant J.98 (5), 912–927. doi: 10.1111/tpj.14279

  • CrossRef
  • Google Scholar

• 40

  KimD.PaggiJ. M.ParkC.BennettC.SalzbergS. L. (2019). Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol.37 (8), 907–915. doi: 10.1038/s41587-019-0201-4

  • CrossRef
  • Google Scholar

• 41

  KimE.-D.SungS. (2012). Long noncoding RNA: unveiling hidden layer of gene regulatory networks. Trends Plant Science.17 (1), 16–21. doi: 10.1016/j.tplants.2011.10.008

  • CrossRef
  • Google Scholar

• 42

  KozomaraA.BirgaoanuM.Griffiths-JonesS. (2019). miRBase: From microRNA sequences to function. Nucleic Acids Res.47 (D1), D155–DD62. doi: 10.1093/nar/gky1141

  • CrossRef
  • Google Scholar

• 43

  KramerM. F. (2011). Stem-loop RT-qPCR for miRNAs. Curr. Protoc. Mol. Biol.95 (1), 15.0.1–15.0.0. doi: 10.1002/0471142727.mb1510s95

  • CrossRef
  • Google Scholar

• 44

  LagesenK.HallinP.RødlandE. A.StaerfeldtH.-H.RognesT.UsseryD. W. (2007). RNAmmer: Consistent and rapid annotation of ribosomal RNA genes. Nucleic Acids Res.35 (9), 3100–3108. doi: 10.1093/nar/gkm160

  • CrossRef
  • Google Scholar

• 45

  LangfelderP.HorvathS. (2008). WGCNA: An r package for weighted correlation network analysis. BMC Bioinf.9 (1), 559. doi: 10.1186/1471-2105-9-559

  • CrossRef
  • Google Scholar

• 46

  LefortV.LonguevilleJ.-E.GascuelO. (2017). SMS: Smartmodel selection in PhyML. Mol. Biol. Evolution.34 (9), 2422–2424. doi: 10.1093/molbev/msx149

  • CrossRef
  • Google Scholar

• 47

  LiH. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv. arXiv:1303.3997v2 [q-bio.GN]. doi: 10.48550/arXiv.1303.3997

  • CrossRef
  • Google Scholar

• 48

  LiJ.ChenC.MaX.GengG.LiuB.ZhangY.et al. (2016). Long noncoding RNA NRON contributes to HIV-1 latency by specifically inducing tat protein degradation. Nat. Commun.7 (1), 11730. doi: 10.1038/ncomms11730

  • CrossRef
  • Google Scholar

• 49

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al. (2009). The sequence alignment/map format and SAMtools. Bioinf. (Oxford England).25 (16), 2078–2079. doi: 10.1093/bioinformatics/btp352

  • CrossRef
  • Google Scholar

• 50

  LiW.HartungJ. S.LevyL. (2006). Quantitative real-time PCR for detection and identification of Candidatus liberibacter species associated with citrus huanglongbing. J. Microbiological Methods66 (1), 104–115. doi: 10.1016/j.mimet.2005.10.018

  • CrossRef
  • Google Scholar

• 51

  LiuS.KracherB.ZieglerJ.BirkenbihlR. P.SomssichI. E. (2015). Negative regulation of ABA signaling by WRKY33 is critical for Arabidopsis immunity towards Botrytis cinerea 2100. eLife.4, e07295. doi: 10.7554/eLife.07295

  • CrossRef
  • Google Scholar

• 52

  López-MauryL.MargueratS.BählerJ. (2008). Tuning gene expression to changing environments: From rapid responses to evolutionary adaptation. Nat. Rev. Genet.9 (8), 583–593. doi: 10.1038/nrg2398

  • CrossRef
  • Google Scholar

• 53

  LoveM. I.HuberW.AndersS. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.15 (12), 550. doi: 10.1186/s13059-014-0550-8

  • CrossRef
  • Google Scholar

• 54

  LoweT. M.EddyS. R. (1997). tRNAscan-SE: A program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res.25 (5), 955–964. doi: 10.1093/nar/25.5.955

  • CrossRef
  • Google Scholar

• 55

  LunaE.PastorV.RobertJ.FlorsV.Mauch-ManiB.TonJ. (2010). Callose deposition: A multifaceted plant defense response. Mol. Plant-Microbe Interactions®.24 (2), 183–193. doi: 10.1094/MPMI-07-10-0149

  • CrossRef
  • Google Scholar

• 56

  McKennaA.HannaM.BanksE.SivachenkoA.CibulskisK.KernytskyA.et al. (2010). The genome analysis toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res.20 (9), 1297–1303. doi: 10.1101/gr.107524.110

  • CrossRef
  • Google Scholar

• 57

  MercerT. R.DingerM. E.MattickJ. S. (2009). Long non-coding RNAs: Insights into functions. Nat. Rev. Genet.10 (3), 155–159. doi: 10.1038/nrg2521

  • CrossRef
  • Google Scholar

• 58

  MonzóC.StanslyP. A. (2020). Economic value of conservation biological control for management of the Asian citrus psyllid, vector of citrus huanglongbing disease. Pest Manage. Science.76 (5), 1691–1698. doi: 10.1002/ps.5691

  • CrossRef
  • Google Scholar

• 59

  MoreauC.GautratP.FrugierF. (2021). Nitrate-induced CLE35 signaling peptides inhibit nodulation through the SUNN receptor and miR2111 repression. Plant Physiol.185 (3), 1216–1228. doi: 10.1093/plphys/kiaa094

  • CrossRef
  • Google Scholar

• 60

  Moreno-HagelsiebG.LatimerK. (2008). Choosing BLAST options for better detection of orthologs as reciprocal best hits. Bioinformatics.24 (3), 319–324. doi: 10.1093/bioinformatics/btm585

  • CrossRef
  • Google Scholar

• 61

  NoboriT.TsudaK. (2019). The plant immune system in heterogeneous environments. Curr. Opin. Plant Biol.50, 58–66. doi: 10.1016/j.pbi.2019.02.003

  • CrossRef
  • Google Scholar

• 62

  NoutoshiY.KuromoriT.WadaT.HirayamaT.KamiyaA.ImuraY.et al. (2006). Loss of NECROTIC SPOTTED LESIONS 1 associates with cell death and defense responses in arabidopsis thaliana. Plant Mol. Biol.62 (1), 29–42. doi: 10.1007/s11103-006-9001-6

  • CrossRef
  • Google Scholar

• 63

  PangK. C.FrithM. C.MattickJ. S. (2006). Rapid evolution of noncoding RNAs: lack of conservation does not mean lack of function. Trends Genet.22 (1), 1–5. doi: 10.1016/j.tig.2005.10.003

  • CrossRef
  • Google Scholar

• 64

  PaniA.MahapatraR. K. (2013). Computational identification of microRNAs and their targets in Catharanthus roseus expressed sequence tags. Genomics Data.1, 2–6. doi: 10.1016/j.gdata.2013.06.001

  • CrossRef
  • Google Scholar

• 65

  PengZ.BredesonJ. V.WuG. A.ShuS.RawatN.DuD.et al. (2020). A chromosome-scale reference genome of trifoliate orange (Poncirus trifoliata) provides insights into disease resistance, cold tolerance and genome evolution in citrus. Plant J.104 (5), 1215–1232. doi: 10.1111/tpj.14993

  • CrossRef
  • Google Scholar

• 66

  PerteaM.KimD.PerteaG. M.LeekJ. T.SalzbergS. L. (2016). Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and ballgown. Nat. Protoc.11 (9), 1650–1667. doi: 10.1038/nprot.2016.095

  • CrossRef
  • Google Scholar

• 67

  PerteaG.PerteaM. (2020). GFF utilities: GffRead and GffCompare. F1000Research9 (304), 304. doi: 10.12688/f1000research.23297.1

  • CrossRef
  • Google Scholar

• 68

  PetersonL. E. (2009). K-Nearest neighbor. Scholarpedia4 (2), 1883. doi: 10.4249/scholarpedia.1883

  • CrossRef
  • Google Scholar

• 69

  PruittK. D.TatusovaT.MaglottD. R. (2005). NCBI reference sequence (RefSeq): A curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res.33 (suppl_1), D501–D5D4. doi: 10.1093/nar/gki025

  • CrossRef
  • Google Scholar

• 70

  RamaduguC.KeremaneM. L.HalbertS. E.DuanY. P.RooseM. L.StoverE.et al. (2016). Long-term field evaluation reveals huanglongbing resistance in citrus relatives. Plant Disease.100 (9), 1858–1869. doi: 10.1094/PDIS-03-16-0271-RE

  • CrossRef
  • Google Scholar

• 71

  RansohoffJ. D.WeiY.KhavariP. A. (2018). The functions and unique features of long intergenic non-coding RNA. Nat. Rev. Mol. Cell Biol.19 (3), 143–157. doi: 10.1038/nrm.2017.104

  • CrossRef
  • Google Scholar

• 72

  ReyesT. H.MaekawaS.SatoT.YamaguchiJ. (2015). The arabidopsis ubiquitin ligase ATL31 is transcriptionally controlled by WRKY33 transcription factor in response to pathogen attack. Plant Biotechnol. 32 (1), 11–19. doi: 10.5511/plantbiotechnology.14.1201b

  • CrossRef
  • Google Scholar

• 73

  RinnJ. L.ChangH. Y. (2012). Genome regulation by long noncoding RNAs. Annu. Rev. Biochem.81 (1), 145–166. doi: 10.1146/annurev-biochem-051410-092902

  • CrossRef
  • Google Scholar

• 74

  ScheuermannJ. C.BoyerL. A. (2013). Getting to the heart of the matter: Long non-coding RNAs in cardiac development and disease. EMBO J.32 (13), 1805–1816. doi: 10.1038/emboj.2013.134

  • CrossRef
  • Google Scholar

• 75

  SharmaY.SharmaA.SinghK.UpadhyayS. K. (2022). Long non-coding RNAs as emerging regulators of pathogen response in plants. Non-coding RNA.8 (1), 4. doi: 10.3390/ncrna8010004

  • CrossRef
  • Google Scholar

• 76

  ShirahamaS.MikiA.KaburakiT.AkimitsuN. (2020). Long non-coding RNAs involved in pathogenic infection. Front. Genet.11 (454). doi: 10.3389/fgene.2020.00454

  • CrossRef
  • Google Scholar

• 77

  SongL.FangY.ChenL.WangJ.ChenX. (2021). Role of non-coding RNAs in plant immunity. Plant Commun.2 (3), 100180. doi: 10.1016/j.xplc.2021.100180

  • CrossRef
  • Google Scholar

• 78

  SzcześniakM. W.RosikiewiczW.MakałowskaI. (2016). CANTATAdb: A collection of plant long non-coding RNAs. Plant Cell Physiol.57 (1), e8–ee. doi: 10.1093/pcp/pcv201

  • CrossRef
  • Google Scholar

• 79

  SzklarczykD.FranceschiniA.WyderS.ForslundK.HellerD.Huerta-CepasJ.et al. (2015). STRING v10: Protein–protein interaction networks, integrated over the tree of life. Nucleic Acids Res.43 (D1), D447–DD52. doi: 10.1093/nar/gku1003

  • CrossRef
  • Google Scholar

• 80

  TangY.QuZ.LeiJ.HeR.AdelsonD. L.ZhuY.et al. (2021). The long noncoding RNA FRILAIR regulates strawberry fruit ripening by functioning as a noncanonical target mimic. PloS Genet.17 (3), e1009461. doi: 10.1371/journal.pgen.1009461

  • CrossRef
  • Google Scholar

• 81

  TarasovA.VilellaA. J.CuppenE.NijmanI. J.PrinsP. (2015). Sambamba: fast processing of NGS alignment formats. Bioinformatics.31 (12), 2032–2034. doi: 10.1093/bioinformatics/btv098

  • CrossRef
  • Google Scholar

• 82

  ThimmO.BläsingO.GibonY.NagelA.MeyerS.KrügerP.et al. (2004). Mapman: A user-driven tool to display genomics data sets onto diagrams of metabolic pathways and other biological processes. Plant J.37 (6), 914–939. doi: 10.1111/j.1365-313X.2004.02016.x

  • CrossRef
  • Google Scholar

• 83

  TsudaK.KatagiriF. (2010). Comparing signaling mechanisms engaged in pattern-triggered and effector-triggered immunity. Curr. Opin. Plant Biol.13 (4), 459–465. doi: 10.1016/j.pbi.2010.04.006

  • CrossRef
  • Google Scholar

• 84

  UlitskyI. (2016). Evolution to the rescue: Using comparative genomics to understand long non-coding RNAs. Nat. Rev. Genet.17 (10), 601–614. doi: 10.1038/nrg.2016.85

  • CrossRef
  • Google Scholar

• 85

  WangN. (2019). The citrus huanglongbing crisis and potential solutions. Mol. Plant12 (5), 607–609. doi: 10.1016/j.molp.2019.03.008

  • CrossRef
  • Google Scholar

• 86

  WangY.FanX.LinF.HeG.TerzaghiW.ZhuD.et al. (2014). Arabidopsis noncoding RNA mediates control of photomorphogenesis by red light. Proc. Natl. Acad. Sci.111 (28), 10359. doi: 10.1073/pnas.1409457111

  • CrossRef
  • Google Scholar

• 87

  WangS.LiuS.LiuL.LiR.GuoR.XiaX.et al. (2020). miR477 targets the phenylalanine ammonia-lyase gene and enhances the susceptibility of the tea plant (Camellia sinensis) to disease during Pseudopestalotiopsis species infection. Planta.251 (3), 59. doi: 10.1007/s00425-020-03353-x

  • CrossRef
  • Google Scholar

• 88

  WangC.-Y.LiuS.-R.ZhangX.-Y.MaY.-J.HuC.-G.ZhangJ.-Z. (2017). Genome-wide screening and characterization of long non-coding RNAs involved in flowering development of trifoliate orange (Poncirus trifoliata l. raf.). Sci. Rep.7 (1), 43226. doi: 10.1038/srep43226

  • CrossRef
  • Google Scholar

• 89

  WeiX.MiraA.YuQ.GmitterF. G. (2021). The mechanism of citrus host defense response repression at early stages of infection by feeding of diaphorina citri transmitting candidatus liberibacter asiaticus. Front. Plant Sci.12. doi: 10.3389/fpls.2021.635153

  • CrossRef
  • Google Scholar

• 90

  WierzbickiA. T. (2012). The role of long non-coding RNA in transcriptional gene silencing. Curr. Opin. Plant Biol.15 (5), 517–522. doi: 10.1016/j.pbi.2012.08.008

  • CrossRef
  • Google Scholar

• 91

  WuG. A.ProchnikS.JenkinsJ.SalseJ.HellstenU.MuratF.et al. (2014). Sequencing of diverse mandarin, pummelo and orange genomes reveals complex history of admixture during citrus domestication. Nat. Biotechnol.32 (7), 656–662. doi: 10.1038/nbt.2906

  • CrossRef
  • Google Scholar

• 92

  WuG. A.TerolJ.IbanezV.López-GarcíaA.Pérez-RománE.BorredáC.et al. (2018). Genomics of the origin and evolution of citrus. Nature.554 (7692), 311–316. doi: 10.1038/nature25447

  • CrossRef
  • Google Scholar

• 93

  WuH.-J.WangZ.-M.WangM.WangX.-J. (2013). Widespread long noncoding RNAs as endogenous target mimics for microRNAs in plants. Plant Physiol.161 (4), 1875–1884. doi: 10.1104/pp.113.215962

  • CrossRef
  • Google Scholar

• 94

  XiaoF.Mark GoodwinS.XiaoY.SunZ.BakerD.TangX.et al. (2004). Arabidopsis CYP86A2 represses Pseudomonas syringae type III genes and is required for cuticle development. EMBO J.23 (14), 2903–2913. doi: 10.1038/sj.emboj.7600290

  • CrossRef
  • Google Scholar

• 95

  XinM.WangY.YaoY.SongN.HuZ.QinD.et al. (2011). Identification and characterization of wheat long non-protein coding RNAs responsive to powdery mildew infection and heat stress by using microarray analysis and SBS sequencing. BMC Plant Biol.11 (1), 61. doi: 10.1186/1471-2229-11-61

  • CrossRef
  • Google Scholar

• 96

  XuW.YangT.WangB.HanB.ZhouH.WangY.et al. (2018). Differential expression networks and inheritance patterns of long non-coding RNAs in castor bean seeds. Plant J.95 (2), 324–340. doi: 10.1111/tpj.13953

  • CrossRef
  • Google Scholar

• 97

  YangY.LiuT.ShenD.WangJ.LingX.HuZ.et al. (2019). Tomato yellow leaf curl virus intergenic siRNAs target a host long noncoding RNA to modulate disease symptoms. PloS Pathogens.15 (1), e1007534. doi: 10.1371/journal.ppat.1007534

  • CrossRef
  • Google Scholar

• 98

  YanJ.YuanF.LongG.QinL.DengZ. (2012). Selection of reference genes for quantitative real-time RT-PCR analysis in citrus. Mol. Biol. Rep.39 (2), 1831–1838. doi: 10.1007/s11033-011-0925-9

  • CrossRef
  • Google Scholar

• 99

  YoonJ.-H.AbdelmohsenK.GorospeM. (2013). Posttranscriptional gene regulation by long noncoding RNA. J. Mol. Biol.425 (19), 3723–3730. doi: 10.1016/j.jmb.2012.11.024

  • CrossRef
  • Google Scholar

• 100

  YuanJ.LiJ.YangY.TanC.ZhuY.HuL.et al. (2018). Stress-responsive regulation of long non-coding RNA polyadenylation in oryza sativa. Plant J.93 (5), 814–827. doi: 10.1111/tpj.13804

  • CrossRef
  • Google Scholar

• 101

  YuQ.ChenC.DuD.HuangM.YaoJ.YuF.et al. (2017). Reprogramming of a defense signaling pathway in rough lemon and sweet orange is a critical element of the early response to ‘Candidatus liberibacter asiaticus’. Horticulture Res.4 (1), 17063. doi: 10.1038/hortres.2017.63

  • CrossRef
  • Google Scholar

• 102

  ZaynabM.FatimaM.AbbasS.UmairM.SharifY.RazaM. A. (2018). Long non-coding RNAs as molecular players in plant defense against pathogens. Microbial Pathogenesis.121, 277–282. doi: 10.1016/j.micpath.2018.05.050

  • CrossRef
  • Google Scholar

• 103

  ZhangM.DuanY.ZhouL.TurechekW. W.StoverE.PowellC. A. (2010). Screening molecules for control of citrus huanglongbing using an optimized regeneration system for ‘Candidatus liberibacter asiaticus’-infected periwinkle (Catharanthus roseus) cuttings. Phytopathology.100 (3), 239–245. doi: 10.1094/PHYTO-100-3-0239

  • CrossRef
  • Google Scholar

• 104

  ZhangL.WangM.LiN.WangH.QiuP.PeiL.et al. (2018). Long noncoding RNAs involve in resistance to Verticillium dahliae, a fungal disease in cotton. Plant Biotechnol. J.16 (6), 1172–1185. doi: 10.1111/pbi.12861

  • CrossRef
  • Google Scholar

• 105

  ZhengZ.QamarS. A.ChenZ.MengisteT. (2006). Arabidopsis WRKY33 transcription factor is required for resistance to necrotrophic fungal pathogens. Plant J.48 (4), 592–605. doi: 10.1111/j.1365-313X.2006.02901.x

  • CrossRef
  • Google Scholar

• 106

  ZhouY.ZhouB.PacheL.ChangM.KhodabakhshiA. H.TanaseichukO.et al. (2019). Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat. Commun.10 (1), 1523–. doi: 10.1038/s41467-019-09234-6

  • CrossRef
  • Google Scholar

• 107

  ZhuQ.-H.FanL.LiuY.XuH.LlewellynD.WilsonI. (2014). miR482 regulation of NBS-LRR defense genes during fungal pathogen infection in cotton. PloS One8 (12), e84390. doi: 10.1371/journal.pone.0084390

  • CrossRef
  • Google Scholar

• 108

  ZhuQ.-H.StephenS.TaylorJ.HelliwellC. A.WangM.-B. (2014). Long noncoding RNAs responsive to Fusarium oxysporum infection in arabidopsis thaliana. New Phytol.201 (2), 574–584. doi: 10.1111/nph.12537

  • CrossRef
  • Google Scholar