Next-Generation Sequencing Identification and Characterization of MicroRNAs in Dwarfed Citrus Trees Infected With Citrus Dwarfing Viroid in High-Density Plantings

Citrus dwarfing viroid (CDVd) induces stunting on sweet orange trees [Citrus sinensis (L.) Osbeck], propagated on trifoliate orange rootstock [Citrus trifoliata (L.), syn. Poncirus trifoliata (L.) Raf.]. MicroRNAs (miRNAs) are a class of non-coding small RNAs (sRNAs) that play important roles in the regulation of tree gene expression. To identify miRNAs in dwarfed citrus trees, grown in high-density plantings, and their response to CDVd infection, sRNA next-generation sequencing was performed on CDVd-infected and non-infected controls. A total of 1,290 and 628 miRNAs were identified in stem and root tissues, respectively, and among those, 60 were conserved in each of these two tissue types. Three conserved miRNAs (csi-miR479, csi-miR171b, and csi-miR156) were significantly downregulated (adjusted p-value < 0.05) in the stems of CDVd-infected trees compared to the non-infected controls. The three stem downregulated miRNAs are known to be involved in various physiological and developmental processes some of which may be related to the characteristic dwarfed phenotype displayed by CDVd-infected C. sinensis on C. trifoliata rootstock field trees. Only one miRNA (csi-miR535) was significantly downregulated in CDVd-infected roots and it was predicted to target genes controlling a wide range of cellular functions. Reverse transcription quantitative polymerase chain reaction analysis performed on selected miRNA targets validated the negative correlation between the expression levels of these targets and their corresponding miRNAs in CDVd-infected trees. Our results indicate that CDVd-responsive plant miRNAs play a role in regulating important citrus growth and developmental processes that may participate in the cellular changes leading to the observed citrus dwarf phenotype.

MicroRNAs have essential functions in plant development and are involved in regulating a myriad of plant processes such as leaf, root, stem, and floral organ morphogenesis and development, biosynthesis, metabolism, homeostasis, vegetative to reproductive growth transition, senescence, signal transduction, and response to biotic and abiotic stress. Upon expression and processing, plant miRNAs are incorporated into the activated RNA-induced silencing complex (RISC) to target RNAs which are complementary to the miRNA guide strand. Once the activated miRNA-RISC complex finds the complementary plant mRNA, it silences the target via RNA degradation or translational repression   (Figure 1A).
The Citrus genus (family Rutaceae), includes several cultivars of high economic value including oranges, mandarins, grapefruits and lemons (2018-2019 US citrus crop packinghousedoor equivalent $3.35 billion) (USDA-NASS, 2019). Citrus flavors and aromas are among the most recognizable and preferred worldwide. In addition, citrus fruits are a rich source of vitamins, antioxidants, minerals, and dietary fiber essential for overall nutritional wellbeing (Van Duyn and Pivonka, 2000;Yao et al., 2004). Citrus trees are produced by grafting a desired scion variety onto a suitable rootstock species that then are planted in commercial citrus orchards. Tree spacing in citrus orchards has varied depending upon the cultivated species and a variety of factors such as soil type, climatic conditions and available farming equipment (Platt, 1973;Tucker and Wheaton, 1978). The historical global trend of citrus orchard spacing has been toward higher tree densities to maintain yield on the reduced available agricultural land and to increase economic returns.
Citrus dwarfing viroid (CDVd) infection of navel orange trees [Citrus sinensis (L.) Osb.] propagated on trifoliate orange [Citrus trifoliata (L.), syn. Poncirus trifoliata (L.) Raf.] rootstock has been previously reported to reduce canopy volume by approximately 50% (Figure 1B) (Vidalakis et al., 2011) and we recently demonstrated that the observed reduction in tree size results from a > 20% reduction in the apical growth of individual shoots within the tree canopy (Lavagi-Craddock et al., 2020). Understanding the molecular mechanism of the CDVd-induced citrus tree size reduction, will be most valuable as it could provide information on how to systematically produce dwarf trees for high density plantings.
To date, very few studies exist of miRNAs in citrus and even fewer in navel orange trees (Lu et al., 2015;Ma et al., 2016;Liang et al., 2017;Xie et al., 2017;Huang et al., 2019) and to our knowledge, there are no published studies on citrus miRNAs in response to viroid infection of citrus field trees. To explore the effect of CDVd-infection on citrus miRNAs and gain insight into the symptom development mechanism leading to the dwarfed phenotype observed in field plantings, we analyzed the effect of CDVd infection using next-generation sequencing (NGS) approach. The increasing number of miRNAs deposited in the miRBase database (Kozomara and Griffiths-Jones, 2014;Kozomara et al., 2019) from a wide range of species (<200), including C. sinensis, enables the discovery of novel miRNAs and their responses to pathogen infection, which may account for the observed species specific reactions and symptom development. Many plant miRNAs are conserved (Axtell and Bartel, 2005) but some are species specific (Moxon et al., 2008) and expressed at lower levels, thus making NGS the ideal approach to discover them and study their expression profiles (Jagadeeswaran et al., 2010;Motameny et al., 2010). Indeed, miRNAs from different plant species such as maize (Zhang et al., 2009), potato (Zhang et al., 2013), peanut , barley (Ferdous et al., 2017), soybean (Zhang et al., 2008), and hop (Mishra et al., 2016), have been identified using NGS approaches.
In this study, we analyzed sRNA libraries prepared from field grown CDVd-infected navel orange and non-infected control trees to characterize miRNAs in the C. sinensis (stems) and C. trifoliata (roots) genomes and their expression profile in response to CDVd infection. This work provides valuable information at the molecular level and establishes Osbeck] on trifoliate orange [Citrus trifoliata (L.) Raf.] rootstock, citrus tree and regulation of target messenger RNA (mRNA) by host miRNA. miRNA genes are transcribed by RNA polymerase II into primary miRNA (pri-miRNA) that folds into self-complementary stem-loop structures. The pri-RNA is then processed by the dicer-like (DCL) RNAse III endonucleases to form precursor miRNA (pre-miRNA), and the pre-miRNA is cleaved again to generate a miRNA/miRNA duplex. One strand of the duplex is degraded and the guide strand is incorporated into the RNA-inducing silencing effector complex (RISC). RISC, with RNAse H activity, targets RNA complementary to the miRNA guide strand. Once the complex finds the complementary plant mRNA, cleaves it, thus regulating its gene expression levels (for a review see Wang et al., 2019). (B) Viroid-infected, dwarf phenotype (∼2 m height) of navel orange trees on trifoliate orange rootstock and the production of viroid derived small RNAs (vdsRNAs). The highly structured CDVd RNA molecule and double stranded intermediate RNAs produced during the viroid rolling-circle replication, are processed by DCL into 21-22 nt vdsRNAs. vdsRNAs are incorporated into RISC and can be guided to form either vdsRNA/host mRNA or vdsRNA/host miRNA complexes. As a result, CDVd vdsRNAs could alter the expression levels of plant mRNA (directly or indirectly via cleavage of plant miRNAs), therefore resulting in the observed citrus tree dwarf phenotype in the field (for a review see Dadami et al., 2017). the foundational framework that is necessary to dissect the subcellular mechanisms responsible for the observed citrus dwarf phenotype in the field.
Stem and root samples were processed in the field and immediately frozen in liquid nitrogen. For each tree, eight stem samples from around the canopy were collected. Leaves and petioles were removed, the stems were roughly chopped into approximately 0.5-1 cm pieces, placed into 50 ml conical tubes, and flash frozen. Root samples were collected from around the tree, at approximately 1 m away from the trunk and 20 cm deep, near the irrigation emitters, using a corer. The roots from eight core soil samples were washed thoroughly with water, gently blotted dry with paper towels, chopped into 0.5-1 cm pieces, placed into 50 ml conical tubes, and flash frozen. In between each sample collection and processing, cutting tools, and working surfaces were sanitized with 10% bleach solution (0.5-1% sodium hypochlorite) and rinsed with water and new sterile disposable plasticware and razor blades were used. Samples were transported into the Citrus Clonal Protection Program (CCPP), Citrus Diagnostic Therapy and Research Laboratory at the UC Riverside (Riverside, CA, United States) on dry ice and stored at −80 • C until analysis.
Total RNA was isolated using the Invitrogen TM TRIzol TM (Thermo Fisher Scientific, Waltham, MA, United States) reagent. For each sample, 300 mg of frozen tissue were ground in liquid nitrogen with mortar and pestle. The ground material was transferred to a 5 ml Eppendorf tube and 3 ml of TRIzol TM reagent was added immediately. RNA extraction was performed according to the manufacturer's instructions. The eluted RNA was aliquoted into four 1.5 ml microcentrifuge tubes to prevent freezing-thawing cycles during downstream analysis. The RNA concentration and quality was assessed with a spectrophotometer and the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, United States) using the Plant RNA Nano assay (RIN values were between 7.9 and 8.6).
Next-Generation Sequencing, sRNA Library Preparation and Sequencing Analysis to Identify Conserved and Novel miRNAs and Their Predicted Targets The sRNA libraries were prepared using the Illumina TruSeq Small RNA Kit (San Diego, CA, United States) following the manufacturer's recommended protocol. The libraries were sequenced using an Illumina HiSeq TM 2500 instrument with single-end 50 bp reads (SeqMatic, Fremont, CA, United States). Raw reads were trimmed to remove low quality bases and adapters using cutadapt v. 1.15 (Martin, 2011) to generate clean sRNAs reads ranging from 18 to 28 nt in length.
The clean reads were then filtered for rRNA, tRNA, snRNA, snoRNA, repeat sequences, and other ncRNAs, using Rfam v.13.0 (Kalvari et al., 2018) with default parameters. The remaining reads were mapped to known miRNAs from the miRBase database (release 21, June 2014) to identify conserved miRNAs Griffiths-Jones, 2011, 2014). The reads were further analyzed to predict potential novel miRNAs using miR-PREFer v. 0.24 using default parameters (Lei and Sun, 2014).
The conserved and novel miRNA sequences were analyzed against C. sinensis mRNA transcripts and C. trifoliata coding sequences (CDS) using psRNATarget v. 2.0 (Dai et al., 2018) to predict potential miRNA-mRNA interactions. DESeq2 v. 1.18 (Love et al., 2014) was used for the differential expression analyses of the miRNAs. The annotated mRNA targets were identified using the Blast2Go (Götz et al., 2008)

Expression Analysis of Citrus miRNAs and miRNA Target Genes Using RT-qPCR
To validate the expression levels of conserved and novel miRNAs, custom stem-loop RT-qPCR assays (catalog number: 4398987) were designed by Thermo Fisher Scientific based on the sequences provided in Supplementary Table 1. For the relative expression quantification, U6 spliceosomal RNA was used as an internal control gene to normalize the efficiency between the target and internal control using the comparative Cq method (Schmittgen and Livak, 2008;Kou et al., 2012). The assay was carried out based on the manufacturer's recommended protocol and all samples were standardized to the same concentration to ensure equal representation. The reverse transcription reactions were performed in a total volume of 15 µL with the TaqMan TM MicroRNA Reverse Transcription Kit (Thermo Fisher) which contained 0.15 µL of 100 mM dNTP, 1 µL of MultiScribe Reverse transcriptase, 1.5 µL of 10x RT Buffer, 0.19 µL of RNase Inhibitor, 4.16 µL of nuclease-free water, 5 µL of total RNA, and 3 µL of 5x RT primer. The reverse transcription reactions were performed with the ProFlex PCR System (Thermo Fisher) as follows: 16 • C for 30 min, 42 • C for 30 min, 85 • C for 5 min, and 4 • C hold. The endpoint qPCR was performed in triplicates, according to the MIQE guidelines (Bustin et al., 2009), on a QuantStudio 12K Flex Real-Time PCR System (Thermo Fisher) with the TaqMan TM Fast Advanced Master Mix (Thermo Fisher) in a total of 20 µL reactions which included: 10 µL of master mix, 7.67 µL of nuclease-free water, 1 µL of TaqMan Small RNA Assay, 1.33 µL of the cDNA template. The endpoint PCR conditions were as follows: 50 • C for 2 min, 95 • C for 20 s, followed by 40 cycles of 95 • C for 1 s, and 60 • C for 20 s.
To verify the relative expression levels of the miRNA target genes, primers for the predicted target genes of miRNAs, were designed for RT-qPCR (Supplementary Table 2). Actin2 was used as an internal control gene to determine the relative abundance of the target mRNA expression levels by the comparative Cq method (Schmittgen and Livak, 2008;Mafra et al., 2012). Reverse transcription was performed using the Invitrogen TM SuperScript TM II Reverse Transcriptase (RT) (Carlsbad, CA, United States). The reaction was performed using the manufacturer's recommended protocol as follows: 1µL of olig (dT) (500 µg/mL), 1 µL of dNTP (10 mM), 2 µL of total RNA and 8 µL of nuclease-free water. The mixture was incubated for 5 min at 65 • C and subsequently chilled on ice. The reaction was prepared with 4 µL of 5x First-Strand Buffer, 2 µL of 0.1M DTT, and 1 µL of RNaseOUT (40 units/µL) and then incubated for 2 min at 42 • C. Finally, 1 µL of SuperScript TM II RT (200 units) was added and the reaction was incubated at 42 • C for 50 min followed by 70 • C for 15 min. Downstream qPCR was also performed in triplicates, according to the MIQE guidelines, using the iTaq Universal SYBR Supermix (Bio-Rad): 10 µL of iTaq Universal SYBR Supermix, 1 µL of cDNA, 0.6 µL of each forward and reverse primers and 7.8 µL of nuclease-free water. The qPCR was performed on the Bio-Rad CFX-96 with the following conditions: 95 • C for 1 min, followed by 40 cycles of 95 • C for 10 s, and 60 • C for 15 s.

Next-Generation Sequencing and Characterization of Potential Citrus miRNAs
To characterize citrus miRNAs and their expression profile in response to CDVd infection, we prepared and analyzed two sRNA libraries from stems and root samples of CDVd-infected and non-infected controls of navel orange citrus trees on trifoliate orange rootstock (Table 1).
From the non-infected trees, 8.1% of the stem and 6.1% of the root were classified as miRNAs. Similarly, for the CDVdinfected stems and roots, 7.2 and 5.8% of the reads, respectively, were classified as miRNAs. The unique unannotated sequences in both the non-infected and CDVd-infected stems represented at least 90% of the total reads while for the roots they represented over 82% ( Table 1). The total unique miRNA reads for both non-infected and CDVd-infected stems represented 0.05% of the reads, while both non-infected and CDVd-infected roots represented 0.0006% of the reads ( Table 1).
The most common size among the total mapped miRNAs sequences ranged between 20 and 24 nt in length, with 21-nt being the predominant miRNA class across different treatments and tissue types. This is consistent with plant antiviral RNAi responses and DCL-mediated processing of dsRNA producing 21 nt siRNAs.

Identification of Conserved miRNAs and Their Expression Profiles
The miRNA sequencing from non-infected control and CDVdinfected stems and roots identified 60 unique conserved miRNAs that ranged from 20 to 24 nt (Supplementary Tables 3, 4). Based on differential expression analysis, four conserved miRNAs (three in the stems and one in the roots) were found to be significantly altered in response to CDVd-infection (P-value and adjusted P-value < 0.05) ( Table 2). Our results indicated that different members of the three miRNA families of interest had different expression levels between the non-infected and the CDVd-infected trees. The conserved miRNA families in the stems included csi-miR156, csi-miR171b, and csi-miR479, while csi-miR535 was the only conserved miRNA found in the roots. The conserved stem miRNAs were moderately more abundant compared to the conserved root miRNAs ( Table 2). All four conserved miRNAs had higher expression levels in the noninfected control than the CDVd-infected trees ( Table 2).
Five miRNA families present in both stem and root tissues were identified: miR166, miR171b, miR399, miR477, and miR482. In stems, the highest represented miRNA families were miR166, and miR399, with five members each, followed by miR171b with four members, miR396, miR477, and miR482 with three members and the remaining 40 miRNAs were represented by a single member (Supplementary Table 3). In the roots, two miRNA families (miR166 and miR399) were represented by five members, five miRNA families (miR167, miR172, miR396, miR477, and miR482) were represented by three members, two miRNA families (miR530 and miR171b) contained two members, and the remaining 40 root miRNA families were represented by a single member (Supplementary Table 4).
Stem-loop RT-qPCR analysis was performed on the four conserved miRNAs in root and stem tissues (csi-miR479, csi-miR156, csi-miR171b, and csi-miR535) from non-infected and CDVd-infected trees to determine their relative abundance. The expression levels of the four conserved miRNAs were significantly altered as a result of CDVd infection. In the stems, csi-miR479's expression decreased 3.55 fold and csi-miR171b had a fold decrease of 2.24, while csi-miR156 had the smallest negative fold change (0.11) (Figure 2). In the roots, csi-miR535 s expression decreased 1.12 fold (Figure 2). The results obtained from the RT-qPCR analysis (Figure 2) were consistent with the NGS read

Identification of Novel miRNAs and Their Expression Profiles
The lengths of the predicted novel miRNAs from stems and roots ranged between 19 and 24 nt. A total of 646 stem and 108 root novel miRNAs were identified (Supplementary Tables 5, 6). No novel root miRNAs had significant differential expression levels.
On the other hand, three novel stem miRNAs (csi-miRNA-75, csi-miRNA-114, and csi-miRNA-435) had significantly different expression levels in response to CDVd infection (P-value and adjusted P-value < 0.05). All three novel stem miRNAs had higher expression levels in the non-infected trees than in the CDVd-infected trees ( Table 3). Stem-loop RT-qPCR was performed to confirm the NGS read frequency results of the identified novel miRNAs in response to CDVd-infection. csi-miRNA-114 showed the largest expression fold change (−9.23), while csi-miRNA-75 and csi-miRNA-435 showed similar fold changes (−1.44 and −1.22, respectively) (Figure 3). These results also support the reliability of RT-qPCR (Figure 3) and NGS read frequencies ( Table 3) as measures of the expression levels of miRNAs.

Citrus miRNA-Target Prediction and Functional Analysis
To understand the function of the identified citrus miRNAs, host target genes were analyzed using the psRNATarget program by cross referencing the results against the C. sinensis genome for the stems (ref: GCF_000317415.1) and the C. trifoliata CDS for the roots (Kawahara et al., 2020). For both conserved and novel stem miRNAs, 83.1% of the miRNA targets were predicted to be regulated by cleavage and 16.9% by translational inhibition (Supplementary Tables 7, 8). Similarly, in the roots, 86.4% of the miRNA targets were predicted to be regulated by cleavage and 13.6% by translational inhibition (Supplementary Tables 9, 10).
Based on the extent of sequence complementarity between miRNAs and their targets, a total of 5,542 potential targets were predicted for the conserved and novel stem and root miRNAs (conserved: stem 63 and root 64; novel: stem 647 and root 109). Of the 5,542 potential miRNA target genes, 494 and 3,926 were targets of the conserved and novel stem miRNAs while 495 and 627 were targets of the conserved and novel root miRNAs, respectively.
Clusters of orthologous groups (COG) functional classification of the targets of conserved and novel miRNAs revealed that the highest proportion of the genes were associated with (i) the nucleus (21% conserved and 11% novel); (ii) the integral component of membrane (14% conserved and 22% novel); and (iii) ATP binding (11% conserved and 10% novel) (Figures 4A,B).
Our data were further annotated based on ontological definitions of the gene ontology (GO) terms, which categorized the predicted targets of the conserved miRNAs differentially expressed in response to CDVd infection into various biological, Frontiers in Microbiology | www.frontiersin.org FIGURE 2 | Differential expression analysis of four identified conserved miRNAs in response to CDVd infection. The relative abundance of each analyzed miRNA in CDVd-infected and non-infected trees was determined using the comparative Cq method by normalization to the U6 spliceosomal RNA. Conserved miRNAs with a significant change (P-value and adjusted P-value < 0.05) were considered to be differentially expressed. The bar graph shows the log 2 fold change of expression levels of the miRNAs in CDVd-infected samples relative to non-infected samples in stem and root tissues. molecular and cellular processes ( Figure 5A). Under the biological process, the predicted targets of conserved miRNA responsive to CDVd infection were subcategorized to (i) metabolic process; (ii) cellular process; (iii) biological regulation; and (iv) regulation of biological processes. The number of sequences associated with these four biological process subcategories had similar values in the stems and roots with the exception of the metabolic process that was higher in the roots ( Figure 5A). For the molecular process, the majority of the predicted target genes of conserved miRNAs responsive to CDVd infection in the roots were subcategorized to catalytic activity while the targets in the stems were mostly subcategorized to binding. Targets belonging to the cellular components category were subcategorized to (i) cellular anatomical entity; and (ii) intracellular subcategories ( Figure 5A) (Supplementary Table 11).
The targets of the predicted novel miRNAs displaying significant differential expression in response to CDVd infection  were categorized to (i) biological process; (ii) molecular process; and (iii) cellular components. The cellular components subcategories (i) cellular anatomical entity; and (ii) intracellular contained most of the sequences (Figure 5B) (Supplementary Table 11).

Expression Profiles and Experimental Validation of miRNA Target Transcripts
The expression levels of eight predicted mRNA targets of the conserved CDVd-responsive miRNAs was determined via RT-qPCR. The results indicate that the expression of the miRNA target genes correlates negatively with the expression of their corresponding miRNA (Figure 6), thus confirming the relationship between CDVd-infection and altered expression levels of specific miRNA targets. Targets of miRNAs belonging to the same miRNA family showed variable results. For example, orange1.1g011651m (Figure 6, bar #5), orange1.1g032310m (Figure 6, bar #6), and orange1.1g008776m (Figure 6, bar #7), which are all targets of different members of the csi-miR156 family, were not uniform and showed fold-change differences.

DISCUSSION
Citrus production exceeded 157.9 million tons in over 9.8 million hectares worldwide for 2019 1 , while in California alone the citrus industry is valued at $3.4 billion dollars with an estimated total economic impact of $7.1 billion (Babcock, 2018). Global decrease in farmland availability, increasing land, water and labor costs, and the continued spread of the deadly Huanglongbing (HLB) disease of citrus, make it imperative to develop tools that allow for high-density citrus plantings for maximization of yields and economic returns per land surface unit. In addition, these factors have forced the citrus industry toward the implementation of novel cultivation practices that would allow for mechanized citrus production under protective structures (Gottwald, 2010;Vidalakis et al., 2011;Lambin, 2012;Verburg et al., 2013;da Graça et al., 2016). The idea of using "graft-transmissible dwarfing agents" for tree size reduction has been continuously investigated in citrus since originally proposed by Cohen (1968) and Mendel (1968) and countries such as Australia and Israel have explored the application of such technology in commercial settings (Broadbent et al., 1992;Bar-Joseph, 1993;Hutton et al., 2000;Semancik, 2003). The observation that CDVd significantly reduced C. sinensis canopy volume on C. trifoliata rootstock (Semancik et al., 1997;Vidalakis et al., 2011); by reducing vegetative growth (Lavagi-Craddock et al., 2020) indicated that CDVd may be used as a possible tool for high-density plantings of citrus, and provided key information on the possible biological mechanism through which CDVd affects specific rootstockscion combinations to reduce tree canopy volume (Vidalakis et al., 2011). Furthermore, understanding the detailed molecular regulatory mechanisms that lead to a reduction in tree canopy volume in response to CDVd infection would provide the necessary knowledge to produce reduced-size citrus trees without the need of a -graft-transmissible viroid agent.
Small RNAs play an essential regulatory role in cellular and plant development functions including antiviral host responses and potentially viroid pathogenesis (Borges and Martienssen, 2015;Dadami et al., 2017;Flores et al., 2017;Czech et al., 2018;Zhu et al., 2018;Treiber et al., 2019;Wang et al., 2019;Ramesh et al., 2020). Two models have been proposed to explain the involvement of RNAi in the pathogenic process induced by viroid infections and both involve vdsRNAs. In the first model, vdsRNAs might act as miRNAs, downregulating the expression of physiologically important host genes, thus inducing disease associated symptoms. vdsRNAs are expected to contain significant identity to a region of the host genome for this model to work and resistance of viroids to RNAi is a feature of the viroid genome (Wang et al., 2004). According to the second model, disease symptoms caused by the nucleus replicating pospiviroids might result from the incorporation of viroid replication intermediates into the transacting small interfering RNA (ta-siRNA) biogenesis pathway. The nucleolus is a ta-siRNA free zone, and mature viroid forms produced in the nucleolus are resistant to degradation. In contrast, (vd)ta-siRNA produced in the nucleus from replication intermediates can then translocate to the cytoplasm where they guide the cleavage of target host mRNA leading to observed symptoms (Gómez et al., 2009). Both models involve viroid secondary structures as a key element that can therefore be interpreted as the evolutionary compromise between the need to interact with host factors and the necessity to survive RNAi. Regardless of whether vdsRNAs are produced according to the first or second model, it is also important to point out that rather than acting directly on host mRNA, vdsRNAs may affect host mRNA targets genes via host miRNAs as previously described (Mishra et al., 2016;Dadami et al., 2017).
Plant disease resistance gene families are typically very large with thousands of members and are commonly considered the putative targets of sRNAs , thus making the study of sRNAs in response to viroid infection a valid approach to investigate the biological mechanisms associated with symptoms. FIGURE 3 | Differential expression analysis of three predicted novel miRNAs in response to CDVd infection. The relative abundance of each analyzed miRNA in CDVd-infected compared to non-infected trees was determined using the comparative Cq method by normalization to the U6 spliceosomal RNA (U6). Novel miRNAs with significant fold changes (P-value and adjusted P-value < 0.05) were considered to be differentially expressed. The bar graph shows log 2 fold changes in expression levels of miRNAs in CDVd-infected samples relative to non-infected samples from stems.
The systematic profiling of sRNAs in CDVd-infected trees, using NGS technologies, was the next logical step to gain insight into the function and regulatory mechanisms of miRNAs through which CDVd may reduce tree canopy size. In this study, we identified conserved and novel miRNAs in citrus and their response to CDVd infection. Consistently with the distribution patterns of sRNAs in other plant species, most sRNAs from both the CDVd-infected and non-infected libraries were found in the 21 and 24 nt classes (Jia et al., 2014;Gao et al., 2015;Mishra et al., 2016;Farooq et al., 2017;Zhang et al., 2018). CDVd-infected stems produced higher frequencies of the 21-nt class than their non-infected counterparts, indicating a CDVd induction of the 21-nt class since the other sRNA classes remained at comparable levels with the non-infected libraries. The increased abundance of the 21-nt class of sRNAs in response to viroid infection observed here is in agreement with previous reports for viral (viroid and virus) infections (Minoia et al., 2014;Zavallo et al., 2015).
All identified differentially expressed conserved and novel miRNAs, in this study, displayed overall reduced expression levels in response to CDVd-infection (Figures 2, 3). Several evolutionary deeply-conserved miRNAs have been shown to retain homologous targets across plant phyla (Axtell et al., 2007) and these include miR156 (stem), miR535 (roots), and miR171b (stem), which represent three out of the four conserved miRNA with differential expression levels in response to CDVd infection identified in this study. In agreement with previous studies, miR156 was shown here to direct the cleavage of squamosa-promoter binding-like protein (SBP) box genes (orange1.1g011651m; orange1.1g032310m; and orange1.1g008776m) (Cardon et al., 1999;Rhoades et al., 2002;Wu and Poethig, 2006;Xie et al., 2006;Gandikota et al., 2007;Riese et al., 2007) (Supplementary Table 7 and Figure 6). Members of this transcription factor family are known to play important roles in flower and fruit development, plant architecture, and in the transitions from juvenile to adult stages and to flowering Yu et al., 2012). Even though miR535 is also known to target squamosa promoter-binding-like protein 3 (Shi et al., 2017;Zhou et al., 2020), our study identified rhomboid-like protein 1 as the target of miR535 in the roots (P_trifoliata_00066_mRNA_51.1). In Arabidopsis, a rhomboidlike protein was identified, providing evidence for the existence of regulated intramembrane proteolysis (RIP), a fundamental mechanism for controlling a wide range of cellular functions, in plants ( Kanaoka et al., 2005). miR171b directs the cleavage of GRAS domain transcription factor genes . However, in our study, we found that miR171b's target, a probable glutathione-S-transferase, was altered in response to CDVd infection (orange1.1g033674m). Glutathione-S-transferases are ubiquitous and multifunctional enzymes encoded by large gene families that can be highly induced by biotic stress including bacterial, fungal, and viral infection (Gullner et al., 2018). We found that the less conserved miR479 cleaves the UDPglucose flavonoid glucosyl-transferase (orange1.1g033614m) (Vogt and Jones, 2000;Offen et al., 2006) which is involved in the process of conjugating hormones, stabilizing secondary metabolites, solubility, transport, and regulating bioavailability of compounds for other metabolic process in Arabidopsis thalianaand DEAD/DEAH box helicases which play an important in regulatory events such as organ maturation and cellular growth and differentiation (orange1.1g028826m and  6 | . Differential expression profile of selected miRNA target genes. The relative gene expression was evaluated by the comparative Cq method using actin2 as a reference gene. The bar graph shows log 2 fold changes of expression levels of target genes in CDVd-infected stems and roots relative to non-infected tissues. The predicted target genes used in the analysis were (1) UDP-glucose flavonoid glucosyl-transferase (orange1.1g033614m, target of csi-miR479-1-stem); (2) DEAD/DEAH box helicase (orange1.1g028826m, target of csi-miR479-2-stem); (3) DEAD/DEAH box helicase (orange1.1g026925m, target of csi-miR479-3-stem); (4) glutathione S-transferase (orange1.1g033674m, target of csi-miR171b-stem); (5) squamosa promoter binding protein-like 2 (orange1.1g011651m, target of csi-miR156-1-stem); (6) squamosa promoter binding protein-like 3 (orange1.1g032310m, target of csi-miR156-2-stem); (7) squamosa promoter-binding protein-like transcription factor family protein (orange1.1g008776m, target of csi-miR156-3-stem); and (8) RHOMBOID-like protein, P_trifoliata_00066_mRNA_51.1, target of csi-miR535-root). orange1.1g026925m) (Macovei et al., 2012). Both target genes may be related to the observed dwarf phenotype (Supplementary Table 11).
The predicted target genes of the novel miRNAs identified in this study  include proteins were associated with target genes including SBP proteins (orange1.1g030599m, orange1.1g029650m, orange1.1g021420m, orange1.1g017256m, orange1.1g016971m, orange1.1g032310m orange1.1g046416m, orange1.1g011662m, orange1.1g010865m, and orange1.1g010591m; target of csi-miRNA-75); plastid-lipid associated protein (orange1.1g030218m, orange1.1g025746m, orange1.1g030180m, and orange1.1g020639m; target of csi-miRNA-114) which are structures that contain lipids and proteins that sequester the overaccumulation of carotenoids during flower development and fruit ripening (Moriguchi et al., 1998;Leitner-Dagan et al., 2006); and vacuolar protein sorting-associated proteins (orange1.1g021304m and orange1.1g017530m; target of csi-miR435), which direct protein cargo from the Golgi apparatus to the vacuoles (Xiang et al., 2013) and have been shown to be important in plant development  (Supplementary Table 11). The results suggest that CDVdinfection affects a wide range of biological functions via different miRNAs. In addition, the GO distribution analysis performed in this study, identified targets of conserved and novel CDVd-responsive miRNAs involved in various processes (Figure 5). Taken together, our findings might suggest that CDVd-infection could lead to developmental reprogramming and growth alterations of citrus trees, leading to the observed symptoms of reduced vegetative growth and overall smaller tree size. Future transcriptome studies could provide additional evidence to elucidate the molecular details in support of this hypothesis.
In this study, the overall miRNA profile of roots (trifoliate orange, C. trifoliata) was not altered in response to CDVd infection to the same extent demonstrated by the stem (navel orange, C. sinensis) miRNA profile. Although CDVd-derived sRNAs were detected in the roots, the trifoliate orange rootstock does not display major symptoms in response to CDVd infection (Vidalakis et al., 2004;Vernière et al., 2006;Murcia et al., 2009). This observation is consistent with the findings that the striking dwarfed citrus tree phenotype caused by CDVd infection results from the reduced vegetative growth of the stems, supporting the hypothesis that the molecular mechanisms responsible for this reprograming must be primarily active in the stems. Future greenhouse studies could provide additional evidence to support this hypothesis.
Finally, the term "Transmissible small nuclear ribonucleic acids" (TsnRNAs) was coined to identify those viroids that do not express a disease syndrome, but rather act as modifying agents of tree performance that result in desirable agronomic traits with potential economic advantages (Semancik et al., 1997;Semancik, 2003). In our study, we did not observe induction and regulation of defense genes via the sRNA pathway in response to CDVd infection, a finding that concurs with the hypothesis that some viroid species might be considered as RNAs modifying cellular functions and plant performance.

CONCLUSION
The expression profile of CDVd-responsive miRNAs indicates that these miRNAs play a role in regulating important citrus tree growth and development processes that may participate in the cellular changes leading to the observed C. sinensis on C. trifoliata rootstock dwarf phenotype.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in the NCBI Sequence Read Archive under accession numbers SRS8100788-SRS8100791 and PRJNA693870.

AUTHOR CONTRIBUTIONS
TD, IL-C, and GV conceived and designed the experiments, analyzed the data, and wrote the manuscript. TD and IL-C performed the experiments with the assistance of SB. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We would like to acknowledge the assistance of the dedicated personnel of the CCPP and the UC-ANR Lindcove Research and Extension Center. We would also like to acknowledge Prof. Danelle Seymour (University of California, Riverside) for providing advice regarding Citrus trifoliata genome. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the University of California. UC is an equal opportunity provider and employer.