Differential Expression of miRNAs Involved in Response to Candidatus Liberibacter asiaticus Infection in Mexican Lime at Early and Late Stages of Huanglongbing Disease

Huanglongbing (HLB) is one of the most destructive diseases threatening citriculture worldwide. This disease has been associated with α-proteobacteria species, namely Candidatus Liberibacter. Due to the unculturable nature of the causal agent, it has been difficult to mitigate the disease, and nowadays a cure is not available. MicroRNAs (miRNAs) are key regulators of gene expression, playing an essential role in abiotic and biotic stress in plants including antibacterial responses. However, knowledge derived from non-model systems including Candidatus Liberibacter asiaticus (CLas)-citrus pathosystem remains largely unknown. In this study, small RNA profiles from Mexican lime (Citrus aurantifolia) plants infected with CLas at asymptomatic and symptomatic stages were generated by sRNA-Seq, and miRNAs were obtained with ShortStack software. A total of 46 miRNAs, including 29 known miRNAs and 17 novel miRNAs, were identified in Mexican lime. Among them, six miRNAs were deregulated in the asymptomatic stage, highlighting the up regulation of two new miRNAs. Meanwhile, eight miRNAs were differentially expressed in the symptomatic stage of the disease. The target genes of miRNAs were related to protein modification, transcription factors, and enzyme-coding genes. Our results provide new insights into miRNA-mediated regulation in C. aurantifolia in response to CLas infection. This information will be useful to understand molecular mechanisms behind the defense and pathogenesis of HLB.


Introduction
Huanglongbing (HLB), or citrus greening, is considered the most devastating and threatening disease in citriculture worldwide, causing a reduction in fruit quantity and quality. The causal agent is a phloem-limited, Gram-negative, α-proteobacterium belonging to the genus Candidatus Liberibacter. It is associated with three species, namely Candidatus Liberibacter asiaticus (CLas), Ca. L. africanus, and Ca. L. americanus [1,2]. The dominant form is associated with CLas, which is mainly transmitted by the Asian citrus psyllid, Diaphorina citri [3]. In Mexico, the first report of HLB was in 2009 in the state of Yucatan; since then, CLas has been detected in 24 of 28 citrus growing states [4]. Typical HLB symptoms include blotchy, mottled, and asymmetrical patterns in leaves, yellow shoots, and malformed fruits with aborted seeds [5,6].

Data Analysis of Small RNA Sequencing
Illumina sequencing generated eight small RNA libraries ranging from 34 to 54 million raw reads (Table 1). After removing adapters and low-quality reads, and size selection, a total of 22 to 37 million clean reads per library were obtained and used for further analysis. Length distribution analysis was performed in clean reads, and as expected, the most abundant reads were 24 nt in length representing 35.40% to 44.64% of the total clean reads, followed by 21,22, and 23 nt (Figure 1). The eight sRNA libraries had a similar length distribution and showed no obvious differences between the asymptomatic (8 wpi) and symptomatic stage (16 wpi). This length distribution pattern is consistent with the biological nature of sRNAs in other plants. Finally, 48.86% to 55.75% of the clean reads were successfully mapped to the Citrus aurantifolia transcriptome. known and novel miRNAs and their target genes at the early and late stages of HLB disease by NGS. These results facilitate the understanding of the role of miRNAs in the bipartite Mexican lime-CLas interaction to elucidate the molecular mechanisms, which will contribute to developing efficient alternatives for disease management.

Data Analysis of Small RNA Sequencing
Illumina sequencing generated eight small RNA libraries ranging from 34 to 54 million raw reads (Table 1). After removing adapters and low-quality reads, and size selection, a total of 22 to 37 million clean reads per library were obtained and used for further analysis. Length distribution analysis was performed in clean reads, and as expected, the most abundant reads were 24 nt in length representing 35.40% to 44.64% of the total clean reads, followed by 21,22, and 23 nt (Figure 1). The eight sRNA libraries had a similar length distribution and showed no obvious differences between the asymptomatic (8 wpi) and symptomatic stage (16 wpi). This length distribution pattern is consistent with the biological nature of sRNAs in other plants. Finally, 48.86% to 55.75% of the clean reads were successfully mapped to the Citrus aurantifolia transcriptome.

Identification of Known and Novel miRNAs in Mexican Lime
ShortStack was used to identify miRNAs in Citrus aurantifolia by considering the fulfillment of the sixteen criteria (Supplementary Table S1). Clusters ranging from 5336 to

Identification of Known and Novel miRNAs in Mexican Lime
ShortStack was used to identify miRNAs in Citrus aurantifolia by considering the fulfillment of the sixteen criteria (Supplementary Table S1). Clusters ranging from 5336 to 5956 sRNAs were identified in the eight libraries mapped to the C. aurantifolia transcriptome. According to ShortStack, approximately between 99.45% and 99.60% of the sRNA clusters were not determined as miRNAs (Supplementary Table S2). Nevertheless, ShortStack is a tool that follows rigorous criteria to determine miRNA annotation. To identify known miRNAs, all the sRNAs were aligned to the database of Plant small RNA genes (https: //plantsmallrnagenes.science.psu.edu accessed on 24 February 2022), allowing no more than two mismatches.
A total of 46 miRNAs from the eight citrus samples were identified by ShortStack, including 29 known miRNAs and 17 novel miRNAs, corresponding to 35 unique sequences ( Table 2). The known miRNAs were classified into 15 families. MiR482 and miR164 were the most abundant families with four members each. MiR156, miR403, and miR391 had three members each and miR399 and miR166 contained two members; furthermore, eight families (miR159, miR172, miR167, miR396, miR162, miR393, miR168, and miR160) had only one member. The size of the majority of the mature miRNAs was 21 nt with 26 sequences, followed by 22, 24, and 20 nt with 12, 6, and 2 mature sequences, respectively. Previously, Axtell and Meyers [23] published the criteria for annotating novel miRNAs, and ShortStack follows those recommendations. According to the length, miRNAs of 23 or 24 nucleotides are scarce, requiring their accumulation in a minimum of four libraries. Here, three of the six 24-nt miRNAs fulfill these criteria. Even so, the minimum free energy (MFE) of the novel miRNA precursors were superior to −53.62 kcal/mol, except cau-miR008, suggesting a strong stability in their secondary structure. The normalized expression of miRNAs in reads per million (RPM) varied significantly, with a range from 0.598 to 69,297.24 (Supplementary Table S3). The known miRNAs showed a higher expression level compared with novel miRNAs. MiR166 and miR159 represented the top two families with a higher number of read counts, with more than 9681.477 RPM per library. Among novel miRNAs, cau-miR002, cau-miR003, and cau-miR009 showed more abundance, ranging from 128.21 to 633.544 RPM. A few miRNAs were identified in only one or two libraries, probably because some miRNAs were expressed at a specific stage of CLas infection or because of low expression levels.  Additionally, the remaining sRNA sequences were compared with the Plant small RNA genes database. We found 27 sRNA sequences that perfectly match with mature miRNAs in other plants (such as miR3951, miR393, and miR399); however, these sequences did not fulfill the criteria evaluated by ShortStack (Supplementary Table S4). Considering ShortStack is a stringent tool, these miRNAs could be false negatives, thus their biological validation and further studies are required.

Differential Expression of miRNAs in Response to CLas in Citrus aurantifolia
To identify miRNAs that respond to CLas infection, differential expression analysis was performed at the early (8 wpi) and late (16 wpi) stages of HLB disease using the read counts obtained by ShortStack. The results showed that 386 and 566 sRNA clusters were differentially expressed in the asymptomatic and symptomatic stages, respectively (Supplementary Table S5). Of the 46 miRNAs identified, at an early stage of the disease, four known (cau-miR403a, cau-miR403b, cau-miR403b, cau-miR172a) and two novel miRNAs (cau-miR006a and cau-miR011) were differentially expressed ( Figure 2a). Additionally, miR3951 did not fulfill the criteria by ShortStack; however, it was identified as a miRNA in Plant small RNA genes and miRBase databases and was thus included for further analysis. On the other hand, eight known miRNAs were differentially expressed at the late stage, corresponding to cau-miR156a, cau-miR156b, cau-miR156c, cau-miR160a, cau-miR482a, cau-miR482c, cau-miR399a, and cau-miR399b. We include the sequences miR399, miR393, and miR3951 that were identified as miRNAs in databases ( Figure 2b). Overall, three miRNAs were upregulated and four were downregulated during the early stage in response to CLas, whereas at the late stage, six miRNAs were upregulated and five were downregulated. Interestingly, during the initial stage of the disease, CLas induced the expression of the two novel miRNAs (cau-miR006a and cau-miR011), being two of the most upregulated miRNAs. Contrariwise, the miR399 family was the most downregulated, with a more than three-fold change. These miRNA expression levels suggested an important role of miRNAs in bacterial infections.

Prediction of Potential Target Genes of the Differentially Expressed miRNAs
To gain further insights into the biological role of miRNAs during CLas infection, we used the transcriptome of C. aurantifolia and the psRNATarget server to identify putative target genes. A total of 436 and 645 target genes were predicted from differentially expressed miRNAs for the asymptomatic and symptomatic stages, respectively. The top 10 candidate target genes based in expectation value of psRNATarget were selected for each miRNA. The seed region between miRNA and target was evaluated, allowing only one mismatch. The hybridization free energy in the asymptomatic stage ranged from −8.36 kcal/mol to −25.08 kcal/mol; meanwhile, in the symptomatic stage, it was −4.58 kcal/mol to −34.76 kcal/mol. Several genes showed unfavorable binding sites (Supplementary Table  S6). This is because the energy required to open nucleotide bonds in the interaction was high, generating a decrease in free energy, causing less stability in miRNA-gene binding. According to the differential expression analysis of target genes obtained by Arce-Leal et al. [15], 16 and 44 genes were discarded for presenting the same expression pattern of their respective miRNA in the asymptomatic and symptomatic stages, respectively (Supplementary Table S6).
Consistent with previous studies, the target genes for most known miRNAs were transcriptional factors, signal transduction, and resistance proteins such as TIR-NBS-LRR (target of cau-miR482), involved in the first line in the detection of pathogens, including bacteria, and other proteins. The selected targets of the novel miRNAs were mainly related to the degradation, modification, and transport of proteins. A total of two target genes were chosen for each differentially expressed miRNA (Table 3). On the other hand, eight known miRNAs were differentially expressed at the late stage, corresponding to cau-miR156a, cau-miR156b, cau-miR156c, cau-miR160a, cau-miR482a, cau-miR482c, cau-miR399a, and cau-miR399b. We include the sequences miR399, miR393, and miR3951 that were identified as miRNAs in databases ( Figure 2b). Overall, three miRNAs were upregulated and four were downregulated during the early stage in response to CLas, whereas at the late stage, six miRNAs were upregulated and five were downregulated. Interestingly, during the initial stage of the disease, CLas induced the expression of the two novel miRNAs (cau-miR006a and cau-miR011), being two of the most upregulated miRNAs. Contrariwise, the miR399 family was the most downregulated, with a more than three-fold change. These miRNA expression levels suggested an important role of miRNAs in bacterial infections.

Prediction of Potential Target Genes of the Differentially Expressed miRNAs
To gain further insights into the biological role of miRNAs during CLas infection, we used the transcriptome of C. aurantifolia and the psRNATarget server to identify putative target genes. A total of 436 and 645 target genes were predicted from differentially expressed miRNAs for the asymptomatic and symptomatic stages, respectively. The top 10 candidate target genes based in expectation value of psRNATarget were selected for each miRNA. The seed region between miRNA and target was evaluated, allowing only one mismatch. The hybridization free energy in the asymptomatic stage ranged from −8.36 kcal/mol to −25.08 kcal/mol; meanwhile, in the symptomatic stage, it was −4.58 kcal/mol to −34.76 kcal/mol. Several genes showed unfavorable binding sites (Supplementary Table S6). This is because the energy required to open nucleotide bonds in the interaction was high, generating a decrease in free energy, causing less stability in miRNA-gene binding. According to the differential expression analysis of target genes obtained by Arce-Leal et al. [15], 16 and 44 genes were discarded for presenting the same expression pattern of their respective miRNA in the asymptomatic and symptomatic stages, respectively (Supplementary Table S6).
Consistent with previous studies, the target genes for most known miRNAs were transcriptional factors, signal transduction, and resistance proteins such as TIR-NBS-LRR (target of cau-miR482), involved in the first line in the detection of pathogens, including bacteria, and other proteins. The selected targets of the novel miRNAs were mainly related to the degradation, modification, and transport of proteins. A total of two target genes were chosen for each differentially expressed miRNA (Table 3).

Validation of miRNAs and Target Genes
To validate the expression of differentially expressed miRNAs due to CLas infection, five miRNAs and six target genes were evaluated by RT-qPCR. The consistency obtained in the expression levels of the selected miRNAs and target genes between the RT-qPCR and the RNA-seq data is shown in Figure 3. In both cases (symptomatic and asymptomatic), the expression pattern of the target genes was opposite to the corresponding miRNA. For example, cau-miR172a was downregulated with a log 2 fold-change of −0.433 in the asymptomatic stage. The result showed by RT-qPCR was −0.389. Meanwhile, its target gene, Citrus_au_0000046354 (transcription factor MYB105), was induced with a log 2 foldchange of 0.446 and 0.787 in sRNA-Seq and RT-qPCR, respectively. This confirms the reliability of the data analysis.

Validation of miRNAs and Target Genes
To validate the expression of differentially expressed miRNAs due to CLas infection, five miRNAs and six target genes were evaluated by RT-qPCR. The consistency obtained in the expression levels of the selected miRNAs and target genes between the RT-qPCR and the RNA-seq data is shown in Figure 3. In both cases (symptomatic and asymptomatic), the expression pattern of the target genes was opposite to the corresponding miRNA. For example, cau-miR172a was downregulated with a log2 fold-change of −0.433 in the asymptomatic stage. The result showed by RT-qPCR was −0.389. Meanwhile, its target gene, Citrus_au_0000046354 (transcription factor MYB105), was induced with a log2 fold-change of 0.446 and 0.787 in sRNA-Seq and RT-qPCR, respectively. This confirms the reliability of the data analysis.

Discussion
Huanglongbing (HLB), or greening disease, is the biggest threat to citriculture, and it spreads widely. When CLas is transmitted to the plant, the bacteria mainly replicate in the roots, while in foliar tissue, the titer is very low and cannot be detected during the first months of infection [32]. Nowadays, all citrus varieties are susceptible to the disease, leading to tree death a few years after infection [33,34]. Despite the efforts to study the disease and the impact of CLas in citrus crops, at present, there is no effective control or cure available. Recently, high-throughput sequencing technologies have developed rapidly,

Discussion
Huanglongbing (HLB), or greening disease, is the biggest threat to citriculture, and it spreads widely. When CLas is transmitted to the plant, the bacteria mainly replicate in the roots, while in foliar tissue, the titer is very low and cannot be detected during the first months of infection [32]. Nowadays, all citrus varieties are susceptible to the disease, leading to tree death a few years after infection [33,34]. Despite the efforts to study the disease and the impact of CLas in citrus crops, at present, there is no effective control or cure available. Recently, high-throughput sequencing technologies have developed rapidly, achieving the identification of known and novel transcripts with high sensitivity in plants [35,36]. Numerous miRNAs have been identified in many plants in response to pathogen attacks. Among the few studies related to miRNA identification in CLas infection [30,31], none of them studied the disease in C. aurantifolia, even though it is one of the most relevant citrus crops. In this study, we obtained sRNA-Seq libraries of C. aurantifolia at the asymptomatic (early, 8 wpi) and symptomatic (late, 16 wpi) stages of HLB development. These libraries were compared with their controls (HLB-) to detect miRNAs associated with CLas infection. The sequencing data showed that the most abundant sRNAs were 24 nt, followed by 21 nt. The lengths between 21 and 24 are consistent with the biological nature of sRNAs [37,38]. The size distributions agree with previous studies, such as C. sinensis in leaves [39], roots [40,41], and fruits [42]; C. trifoliata in fruits and flowers [43]; and C. Junos in roots [44]. Therefore, the sRNA-Seq data obtained here was reliable for further analysis.
More than 22 million filtered reads per library were obtained and then successfully aligned to Mexican lime transcriptome with a range of 48.86-55.77%. The identification of miRNAs by ShortStack provided relevant and reliable information of miRNAs annotation by reducing false positives and increasing precision and accuracy. Previously, the transcriptome of Paspalum notatum was used as a reference in ShortStack, obtaining a 30.3% reads alignment [45]. In addition, 24 miRNAs were detected in L. campestre small RNA data using the L. campestre genome; meanwhile, 29 miRNAs were identified in L. appelianum sRNA data using the L. appelianum transcriptome with ShortStack [46]. We detected 46 miRNAs in C. aurantifolia, corresponding to 29 known miRNAs and 17 novel miRNAs. Known miRNAs were identified in several plant species, being the most abundant sequences. Previous works have reported this association between the abundance of miRNAs and their conservation across the plant kingdom [37,47,48]. Our results show that cau-miR166, cau-miR159, and cau-miR396 were the most abundant families, similar to other works in citrus plants [30,49,50]. Of the novel miRNAs, most had low expression, including 12 miRNAs with a range of 0.598 to 11.627 RPM. In previous studies, novel sequences represent a small percentage of all miRNAs, and most of these miRNAs are present in low abundance [37,51,52].
The differentially expressed miRNAs at the asymptomatic stage showed the repression of cau-miR403 and cau-miR172. MiR403 has been reported in soybean plants (Glycine max) infected with the oomycete Phytophthora sojae [53]. Similarly, Zhang et al. [25] reported miR403 repression in Arabidopsis plants infected with three different strains of Pseudomonas syringae (Pst DC3000 EV, Pst DC3000 hrcC, Pst DC3000 avrRpt2), suggesting a potential role in plant immunity. This is because the target gene of miR403 is Argonaute 2 (AGO2), in association with AGO1 [54]. The regulatory network miR403/AGO2/AGO1 has been validated with viruses, where repression of AGO1 by viral suppressors leads to negative regulation of miR403 and subsequent induction of AGO2 [55]. However, AGO2 induction by bacteria could be mechanistically different from virus-triggered induction [25]. Hence, further studies are required to elucidate the relationship between miR403 and bacterial infection in plants. In our work, AGO2 was one of the targets with the best rate. Nevertheless, the gene had the same expression pattern as the miRNA for which it is a target, thus it was discarded. Previously, miR172 has been shown to play an important role in plant growth and development [56]. Overexpression of this miRNA disrupts normal leaf and flower development [57], associating miR172 with the occurrence of disease symptoms in plants [58]. The target gene of miR172 is MYB105. MYB genes play an important function in plant defense against pathogens and plant hormone responses [59]. A previous work reported that MYB gene expression contributed to defense against Pyricularia oryzae and Xanthomonas oryzae in rice [60].
Among novel miRNAs, cau-miR011 and cau-miR006 represent the upregulated miR-NAs in the early stage of infection. Cau-miR011 represses the expression of SUC3 gene, a protein responsible for sucrose transport in plants. Like most fruit species, sucrose is the main photoassimilate that is synthesized in leaves and then exported through the phloem sieve elements to sink tissue in citrus [61,62]. In a previous study, Fan et al. [63] reported accumulated sucrose in leaves of C. sinensis infected by CLas, concluding that photoassimilate translocation is affected by HLB infection. On the other hand, Martinelli et al. [64] argued that sucrose transporter genes were induced in peel tissue in symptomatic fruits (HLB+). Likewise, sucrose also acts as a signaling compound that can alter gene expression in plant growth and development [65]. It has been well-established that HLB affects carbohydrate metabolism, which suggests that cau-miR011, by regulating the expression of sucrose transporter proteins and inhibiting their action, could have generated the accumulation of sucrose in CLas-infected leaves, resulting in the impairment of photoassimilate translocation. This could contribute to small fruits and physiological disorders in plants, which are typical symptoms of HLB. However, it is necessary to study the biological role of cau-miR011. Cau-miR006 regulates the expression of the FtsH extracellular protease family. This target is an ATP-dependent zinc metalloprotease involved in plant response to stress, maintaining a normal progress and stability in the chloroplast [66]. Loss of FtsH leads to different phenotypes, including pale seedlings, albino seeds, ROS accumulation, and disrupted thylakoid formation [67][68][69].
Regarding miRNAs differentially expressed at the symptomatic stage (16 wpi), the same expression pattern of miR393 and miR160 has already been described in previous works in C. sinensis infected with CLas [30] and Arabidopsis-Pseudomonas interaction [29,70]. Interestingly, in our study, cau-miR399 was the most downregulated miRNA and is considered a specific miRNA of HLB infection. Contrastingly, miR399 was significantly upregulated in C. sinensis, causing the degradation of mRNA encoding UBC24 (pho2) enzyme responsible for phosphorus translocation and remobilization [30]. MiR399 induces its expression in response to phosphorus starvation in HLB-positive samples; however, in our study, miR399 was highly downregulated in the symptomatic stage. Additionally, its target gene (UBC24) was differentially expressed with a 1.277 log 2 fold-change, suggesting that miR399 response is differential among citrus cultivars.
It is known that the ubiquitin-proteosome system is the main protein turnover pathway responsible for protein degradation and modification in eukaryotic cells, contributing to many aspects of cellular processes [71]. Ubiquitination is a process that requires ubiquitinactivating enzyme (E1), ubiquitin-conjugating enzyme (E2), and ubiquitin ligase (E3). Specifically, E2 regulates the interaction, function, and destination of target genes [72,73]. Cau-miR403 induces the expression of ubiquitin-conjugating enzyme E2 (UBC38). Likewise, cau-miR399 induced the expression of ubiquitin-conjugating enzyme E2 (UBC24). This ubiquitin has been described in previous works [30,[74][75][76]. Ubiquitin E3 are classified in specific domains such as U-box, specifically PUB40 degraded BZR1; the functions are related to biotic and abiotic stress [77]. Cau-miR006 downregulated the expression of PUB40. Despite the importance of E2 and E3 enzymes, information and studies are largely lacking [78].
The first defense mechanism of plants is based on the initial recognition of pathogen by pattern recognition receptors (PRRs) receiving signals about the attack to activate the immune response denominated PAMP-triggered immunity (PTI) through protein interactions. Important signaling factors involve mitogen-activated protein kinase (MAPK) cascades [79,80]. The novel miRNA cau-miR011 regulated MAPKKK10 in the early stage of the disease. In addition, the plant hormone auxin is known to be involved in many plant processes; the repression of auxin increases plant resistance to pathogens [27]. Cau-miR160 regulated the expression of ARF17, playing an important role in the infection and suggesting that suppression of ARF could lead to a partial resistance to the bacteria. A second pathway to detect pathogen attack is through effector recognition, leading to the activation of effector-triggered immunity (ETI). MiRNAs regulate hundreds of disease resistance genes (R-genes) by targeting sites with similar nucleotide motifs [81]. Among the predicted target genes, one of the most reported gene was LRR (miR3951), including TIR-NBS-LRR (miR482). In a previous work, the genes LRR1 and LRR2 were induced in tomato plants during infection by Pst DC3000, while the expression of miR482 was downregulated [82]. Our results showed the downregulation of R-genes in both stages of HLB disease.
At the symptomatic stage, cau-miR156a and cau-miR156b were downregulated, while cau-miR156c was upregulated. Precursors and mature miRNA sequences are different. However, the miRNAs share the same family of target genes, corresponding to Squamosa promoter-binding-like protein (SPL). SPL plays an important role in plant development of flowers and fruits. The overexpression of miR156 in Arabidopsis infected with Pst DC3000 induced susceptibility to the bacteria, whereas overexpression of SPL9 increased resistance in the plant [83]. Similarly, a recent report demonstrated that overexpression of OsSPL4 in transgenic rice plants enhanced disease resistance against Magnaporthe oryzae [84].
In the present study we generated a homogeneous and accurate experimental C. aurantifolia-CLas system, by selecting 8 and 16 wpi as asymptomatic and early symptomatic points with the idea to avoid the pleiotropic effects observed in very late infection stages. In both asymptomatic and symptomatic stages, all detected miRNAs are present; however, in order to understand the miRNAs-mediated response during infection, we analyzed only the differential expressed miRNAs. For this reason, only a few of them are common in both stages. Our data are in agreement with a previous study of miRNA differential response in a C. sinensis-CLas system, where not all detected miRNAs are differentially expressed during different disease stages [30]. Our results provide putative functions of miRNAs in response to CLas infection. This information would allow the development of sustainable management strategies, since in the end, the most desirable strategy for disease control in plants will always be the use of natural defense mechanisms.

Plant Material and Experimental Design
Mexican lime plants on alemow (C. macrophylla) rootstock (nine months after grafted) placed in 40 L pots with 20 kg substrate (Coconut coir dust, vermicompost, and perlite, 1:1:1) were used in this study. Plants were irrigated two or three times per week and fertilized every three weeks with a water-soluble fertilizer mixture (20-0-30 N-P-K) and micronutrients (Microfol ® Combi P.S., Biolchim, Bologna, Italy). Inoculation was carried out in February 2018. The inoculum source was symptomatic budwoods with the bacterium Candidatus Liberibacter asiaticus. The experimental design was carried out as previously described by Arce-Leal et al. [15]. A total of 45 plants were inoculated with budwood from CLas-infected Mexican lime trees, and 15 plants were inoculated with CLas-free budwood as controls. Plants were maintained in a 250 m 2 size shaded greenhouse with an anti-aphid insect screen at the INIFAP Experimental Station in Tecomán, Colima, Mexico, with an average annual temperature of 26 • C (range 16-34 • C) and a mean humidity of 64% (range 46-100%).
For this study, HLB disease progress was designated as early or late stage based on bacterial titer and plant symptoms. The early (asymptomatic) stage was designated for plants at eight weeks post-inoculation (8 wpi) of CLas and containing a homogeneous bacterial titer of~2.5 × 10 2 bacterial cells/100 ng of total DNA. The late stage (symptomatic) of HLB disease progress was considered at 16 wpi because plant symptoms at this point included leaf yellowing and asymmetric blotchy mottling [15]

sRNA Library Construction and Sequencing
Total RNA was extracted from the full leaves of the selected plants using the TRIzol ® protocol (Thermo Fisher Scientific, Carlsbad, CA, USA). The concentration and integrity of RNA samples were verified using a Nanodrop 2000 and capillary electrophoresis by the 2100 Bioanalyzer RNA Nano Chip (Agilent Technologies, Inc., Santa Clara, CA, USA). For each condition, total RNA from the 10 selected plants were divided in two groups, each one containing five plants and pooled in an equimolar ratio to construct each cDNA library. A total of eight small RNA libraries were generated, four at the early (two CLas inoculated and two negative controls) and four at the late (two CLas inoculated and two negative controls) stage of HLB disease progress. We generated libraries LM8wpiHLB

Identification of Known and Novel miRNAs in Mexican Lime
Raw data were processed using Atropos v1.1.28 [85] to remove adaptor sequences, low-quality reads, and contaminated reads. Reads smaller than 18 nt or longer than 27 nt were also removed. The quality control standards were evaluated with FASTQC software (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/, accessed on 17 January 2022). Because the genome of Citrus aurantifolia is not available, clean reads were mapped (up to one mismatch) to the reference transcriptome [86] using Bowtie v1.2.3 software [87]. ShortStack v3.8.5 package (https://github.com/MikeAxtell/ShortStack, accessed on 21 January 2022) [21] was used to annotate and quantify the sRNAs previously aligned to the reference. Then, the sRNA sequences were aligned against the databases of Plant small RNA genes (https://plantsmallrnagenes.science.psu.edu, accessed on 24 February 2022) [88] and miRBase (https://www.mirbase.org, accessed on 7 March 2022), with no more than two mismatches to identify known mature miRNAs. For the known miRNAs, the miRNA family was assigned by considering the best match with the database of Plant small RNA genes, but with a new letter suffix to differentiate miRNA sequences. The names of novel miRNAs were assigned sequentially.

Differential Expression Analysis of miRNAs
Shortstack-derived read counts were used to identify differentially expressed miRNAs. The expression of sRNAs in each library was normalized and analyzed by DESeq2 v1.36 package (https://bioconductor.org/packages/release/bioc/html/DESeq2.html, accessed on 15 August 2022) [89]. Normalized data in asymptomatic and symptomatic samples were compared with their corresponding control samples. Those miRNAs with an adjusted p-value ≤ 0.05 were considered as differentially expressed.

Prediction and Annotation of Target Genes
The psRNATarget server (https://www.zhaolab.org/psRNATarget/analysis, accessed on 22 August 2022) [90] was used to predict the potential target genes of the differentially expressed miRNAs at early and late stages using the default parameters and a maximum expectation value of four. The transcriptome of C. aurantifolia [86] was used as a reference for the target search. The top 10 rate targets were chosen, and the criteria used to select the best two targets for each miRNA were as follows: (1) complementarity in the seed region (2-13 nt), (2) the binding properties miRNA-target with RNAUp server (http://rna.tbi. univie.ac.at/cgi-bin/RNAWebSuite/RNAup.cgi, accessed on 21 September 2022) [91], and (3) mRNA expression obtained by RNA-Seq [15]. Putative functions of the final targets were annotated in the NR protein database using BLASTx with default settings.

Validation of miRNAs and Target Genes Expression by RT-qPCR
To validate the expression of differentially expressed miRNAs and target genes, stemloop RT-qPCR [92] and RT-qPCR were performed, respectively. Five candidate miRNAs (two from early stage and three from late stage) and six target genes were selected. Total RNA (one microgram) from leaves of five plants was pooled in equimolar concentration to obtain eight samples, and M-MLV Reverse Transcriptase (Invitrogen) used to produce the cDNA according to the manufacturer's instructions. Each condition was represented by two biological replicates.
Specific primers for the miRNAs and target genes (Supplementary Table S7) were designed using Primer Select software (DNASTAR Lasergene, Madison, WI). qPCR reactions were performed on a CFX96TM real-time PCR system (Bio-Rad) and SsoFast EvaGreen ® Supermix (Bio-Rad, Foster City, CA, USA) according to the manufacturer's instructions. The cycling program included incubation at 95 • C for 30 s followed by 40 cycles of denaturation at 95 • C for 5 s and annealing at 58 • C for 10 s. All reactions were analyzed using two technical replicates. Relative expression levels were calculated using the 2 −∆∆Ct method [93]. COX [94] and U6 [95] genes were used as internal controls for the target genes and the miRNAs, respectively.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12051039/s1, Table S1: Criteria evaluated by ShortStack for miRNA analysis; Table S2: Clusters sequences identified in the eight libraries by ShortStack; Table S3: Abundance of identified miRNAs in Citrus aurantifolia; Table S4: sRNAs sequences identified as miRNAs in Plant small RNA genes; Table S5: sRNAs differentially expressed in asymptomatic and symptomatic stages; Table S6: Expression statistics and free energy of binding between miRNAs and target genes; Table S7: Specific oligonucleotides of miRNAs and target genes for RT-qPCR.

Data Availability Statement:
The raw data presented in this study are openly available in NCBI database, reference number PRJNA574168.