Transcriptome Analysis of Citrus Dwarfing Viroid Induced Dwarfing Phenotype of Sweet Orange on Trifoliate Orange Rootstock

Dwarfed citrus trees for high-density plantings or mechanized production systems will be key for future sustainable citrus production. Citrus trees consist of two different species of scion and rootstock. Therefore, any observed phenotype results from gene expression in both species. Dwarfed sweet orange trees on trifoliate rootstock have been produced using citrus dwarfing viroid (CDVd). We performed RNA-seq transcriptome analysis of CDVd-infected stems and roots and compared them to non-infected controls. The identified differentially expressed genes validated with RT-qPCR corresponded to various physiological and developmental processes that could be associated with the dwarfing phenotype. For example, the transcription factors MYB13 and MADS-box, which regulate meristem functions and activate stress responses, were upregulated in the stems. Conversely, a calcium-dependent lipid-binding protein that regulates membrane transporters was downregulated in the roots. Most transcriptome reprogramming occurred in the scion rather than in the rootstock; this agrees with previous observations of CDVd affecting the growth of sweet orange stems while not affecting the trifoliate rootstock. Furthermore, the lack of alterations in the pathogen defense transcriptome supports the term “Transmissible small nuclear ribonucleic acid,” which describes CDVd as a modifying agent of tree performance with desirable agronomic traits rather than a disease-causing pathogen.


Introduction
Citrus (family Rutaceae) is an economically important crop with several high-value cultivars, including oranges, mandarins, grapefruits, and lemons. The citrus industry was recently valued at USD3.63 billion in California alone, with an estimated economic impact of USD7.6 billion [1]. Citrus fruits and flavors have a wide range of usage and are important sources of vitamins, antioxidants, minerals, and dietary fiber essential for overall nutritional wellbeing [2,3]. Citrus trees are produced by grafting the desired scion variety onto a suitable rootstock species and then planting it in commercial citrus orchards. Tree spacing in citrus orchards has reduced over time in favor of higher tree densities to maintain yield on the reduced available land for agriculture, increase economic returns, and combat the growing threat of citrus huanglongbing (HLB). However, high-density citrus plantings cannot be achieved without the use of dwarfed citrus trees [4][5][6][7].
In addition, CDVd-induced symptoms on citrus hosts depend on species, variety, and rootstock. CDVd infection of navel orange trees [Citrus sinensis (L.) Osbeck] propagated on trifoliate rootstock [C. trifoliata (L.), syn. Poncirus trifoliata (L.) Raf.] has shown that the observed stunting phenotype induced by CDVd-infection reduces canopy volume by approximately 50% (Figure 1), and the reduction in the canopy volume is a result of a >20% decrease in the apical growth of the shoots [19][20][21]. Therefore, there is a need to elucidate the molecular mechanism of the dwarfing phenotype, so the valuable information can be used to produce dwarfed trees for highdensity plantings.
Currently, identifying potential targets of CDVd is challenging due to the nature of viroids, which do not encode for proteins, limiting the experimental approaches. In addition, most of the data on symptom induction by single viroids refer to experimental tests on the biological indicator 'Etrog' citron (C. medica L.). Previously, studies on CDVdinduced differentially expressed genes were performed on 'Etrog' citron in growth chamber conditions; microarray analysis identified mainly genes related to the cell wall structure, amino acid transport, signal transduction, and plant defense/stress response [22]. To this date, there are no transcriptome studies on mature dwarfed citrus field trees infected with CDVd [9].
To further explore the global effect of CDVd-infection on citrus host mRNAs and gain insight into the symptom development mechanism leading to the dwarfed phenotype observed in field plantings, we performed a transcriptome analysis of CDVd-infected dwarf citrus trees using a high throughput sequencing (HTS) approach. There has been limited transcriptomic research on citrus host response to viroid infection using the bioindicator 'Etrog' citron [23], while other viroid studies have used model plant systems such as Solanum lycopersicum and Nicotiana tabacum [11,24,25].
The improved availability of citrus genome and annotation has enabled the widespread use of HTS studies to develop a global understanding of the molecular machinery behind the response to pathogen infection [26]. Furthermore, HTS allows the discovery of rare transcripts, which would be difficult to identify using traditional methods such as microarrays, and is more time efficient than an RT-qPCR approach. More specifically, in citrus, HTS analysis of mRNA has been widely used to understand host responses to pathogens such as 'Candidatus' Liberibacter asiaticus [27][28][29][30] and citrus tristeza virus [28,31].
In this study, we performed a comparative analysis of differentially expressed genes in the stems and roots of CDVd-infected and non-infected sweet orange scion on trifoliate orange rootstock. This analysis provides valuable molecular information to understand the mechanisms responsible for the citrus dwarf phenotype in the field.

Plant Materials and RNA Isolation
Plant materials were collected in April 2016 from 18-year-old 'Parent Washington' navel orange on 'Rich 16-6 trifoliate orange trees. All trees were planted in the same east-west running orchard located at the University of California, Agriculture and Natural Resources, Lindcove Research and Extension Center (LREC) in Exeter, CA, USA. CDVdinfected trees were planted at high density (3 × 6.7 m), whereas non-infected control trees were spaced at standard density (6.1 × 6.7 m).
Stem and root samples were collected from the south side of the dwarfed CDVdinfected (n = 3) and the north and south sides of the full-size non-infected trees (n = 3). This type of sampling was necessary to represent the different sizes of trees (i.e., dwarfed trees 61.2% smaller than full size, [21]) and the different light distribution within their canopies due to the size and east-west orientation of the tree rows, in the transcriptome analysis.
After 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 with liquid nitrogen in the field. Feeder roots were sampled approximately 1 m away from the trunk and 20 cm deep, near the irrigation emitters, using a 10 cm diameter corer. The feeder root 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 with liquid nitrogen in the field. Between each sample collection and processing, cutting tools and working surfaces were sanitized with 10% bleach solution (0.5-1% sodium hypochlorite) and then rinsed with water. New sterile disposable plasticware and razor blades were used to chop each sample. The frozen samples were transported to the Citrus Clonal Protection Program (CCPP), Citrus Diagnostic Therapy and Research Laboratory at the UC Riverside (Riverside, CA, USA) on dry ice and stored at −80 • C for downstream analysis (CDFA permit 3477).
Total RNA was isolated using the TRIzol ® reagent (Thermo Fisher Scientific, Waltham, MA, USA). For each sample, 300 mg of frozen tissue were ground in liquid nitrogen with sterilized mortar and pestle. The ground material was transferred to a 5 mL Eppendorf tube, and 3 mL of TRIzol ® 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 were assessed with a spectrophotometer and the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) using the Plant RNA Nano assay. The presence or absence of CDVd, as well as that of other graft-transmissible pathogens of citrus endemic to California, was confirmed in each sample by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) as previously described [21,32].

Library Preparation and High Throughput Sequencing
Eighteen cDNA libraries were prepared for each sample using the Illumina TruSeq Stranded mRNA Kit (San Diego, CA, USA), following the manufacturer's recommended protocol for Low Sample (LS) throughput. The libraries were quantified with the Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) using the High Sensitivity DNA Kit. The libraries were quantified using the in-house SeqMatic HT1 quantitative polymerase chain reaction (qPCR) assay (SeqMatic, Fremont, CA, USA) and pooled equimolar. The libraries were sequenced using an Illumina HiSeq TM 4000 (San Diego, CA, USA) instrument with paired-end 100 bp reads (SeqMatic, Fremont, CA, USA). Raw reads were trimmed and demultiplexed for subsequent bioinformatic analysis.

Bioinformatic Analysis
The raw sequencing data of the north and south samples from each of the three non-infected full-size trees were concatenated into one file for analysis. Bioinformatic analyses were performed using the OmicsBox software suite version 2.0.36 (Cambridge, MA, USA) using the reference Valencia orange (version 1.0) (GCF_000317415.1) and trifoliate (v. 1.3.1) [33] genomes. RNA-seq alignment was performed using Spliced Transcript Alignment to a Reference (STAR v. 2.7.8a) [34] with default parameters. The quantification of transcript expression levels was performed using the HTSeq (v. 0.90) package [35], while the differential expression analysis was performed edgeR (v. 3.28.0) [36] with default parameters. Functional analysis was performed using the Blast2GO [37] tool within OmicsBox Suite. BLAST [38] was performed on the DEG sequences, and the gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) functions [39] were assigned. Protein functionality was confirmed with InterPro [40]. The dataset was uploaded into NCBI Sequence Read Archive under the accession numbers SAMN26677719 to SAMN26677722. All figures were created using GraphPad Prism version 9.3.0 (San Diego, CA, USA).

Expression Analysis of Citrus mRNA Target Genes Using RT-qPCR
RT-qPCR assays were designed to verify the expression levels of the predicted target genes (Supplementary Table S1). Actin was used as an internal control gene to determine the relative abundance of the target mRNA expression levels by the comparative Cq method. Reverse transcription was performed using the SuperScript TM II Reverse Transcriptase (RT) (Thermo Fisher Scientific). The reactions were performed using the manufacturer's recommended protocol as follows: 2 µL of oligo (dT) (500 µg/mL), 2 µL of dNTP (10 mM), 4 µL of total RNA (diluted to 100 µg/µL), and 16 µL of nuclease-free water with a final volume of 24 µL. All samples were standardized to the same concentration to ensure equal representation, and incubation steps were performed using the ProFlex thermal cycler (Thermo Fisher Scientific). The mixture was incubated for 5 min at 65 • C and subsequently chilled on ice. The first strand synthesis was prepared with 8 µL of 5x First-Strand Buffer, 4 µL of 0.1 M DTT, and 2 µL of RNaseOUT TM (40 units/µL) (Thermo Fisher Scientific) and then incubated for 2 min at 42 • C. Finally, 2 µL of SuperScript TM II RT (200 units) were added, and the reaction was incubated at 42 • C for 50 min followed by 15 min at 70 • C. Downstream qPCR was performed in triplicates, according to the MIQE guidelines [41], using the TaqMan TM Fast Advanced Master Mix Kit (Thermo Fisher Scientific) as follows: 10µL of master mix, 8 µL of nuclease-free water, 1 µL of primer and probe mixture, and 1 µL of cDNA. FAM fluorophore was used for all qPCR probes. The qPCR was performed using the QuantStudio 12K Flex Real-Time PCR (Thermo Fisher Scientific) with the following conditions: 50 • C for 2 min, 95 • C for 2 min, and 95 • C for 1 s, 60 • C for 20 s for 40 cycles.

Transcriptome Assembly and Annotation
To understand the citrus host transcriptome response to CDVd infection, we prepared and analyzed RNA-seq libraries from the stem and root samples of non-infected and CDVd-infected navel orange citrus trees on trifoliate orange rootstock. Overall, stems generated more reads than roots for non-infected and CDVd-infected trees (Table 1). A high percentage of reads (72-88%) from both the non-infected and CDVd-infected trees were mapped to their relative genomes using the STAR alignment software ( Table 1). The average length of mapped reads in all samples tested was close to the expected size of 200 nt (Table 1).

Differentially Expressed Genes (DEGs) Analysis and Identification
The multidimensional scaling (MDS) plot showed distinct differences between the non-infected and CDVd-infected trees in stem and root tissues (Figure 2A,B). Non-infected stems and roots were similar and clustered closer together than the CDVd-infected tissues (Figure 2A,B).  Table 2). Heat maps of the top 50 DEGs identified from the stem and root tissues of non-infected and CDVd-infected trees were generated based on FDR ( Figure 4). Results indicated a relatively even clustering of stem upregulated and downregulated DEGs in the non-infected and CDVd-infected samples. In contrast, the roots showed clustering skewed towards the DEGs that were significantly downregulated in the CDVd-infected and upregulated in the non-infected samples ( Figure 4A,B).  DEG analysis identified up-and downregulated targets. In the stem, the probable 28S rRNA (cytosine-C(5))-methyltransferase (FDR = 5.05 × 10 −11 ; logFC = 8.83) and the pentatricopeptide repeat-containing protein (FDR = 4.29 × 10 −7 ; logFC = −2.95) had the most statistically significant fold change (Table 3).

Confirmation of Candidate DEGs by RT-qPCR Analysis
To confirm the accuracy of Illumina RNA-seq results, RT-qPCR was performed on selected candidate DEGs. The logFC values of the RNA-seq and RT-qPCR were consistent and showed similar trends for both up-and downregulated stem and root DEGs ( Figure 5).

Functional Classification of DEGs
The RNA-seq data were annotated based on the gene ontology (GO) terms, which categorized the DEGs in response to CDVd infection for both stems and roots into various biological, molecular, and cellular processes ( Figure 6). Among the upregulated DEGs in the stems, there were 143 associated with biological process (BP), 82 14). Only two GO terms were associated with CC (i) cellar anatomical entity (up 26 and down 110) and (ii) protein-containing complex (up 2 and down 8) ( Figure 6B). A list of all GO term annotations for the root DEGs can be found in Supplementary Tables S8 and S9. Kyoto Encyclopedia of Genes and Genome (KEGG) analysis was performed on upand downregulated DEGs from non-infected and CDVd-infected trees from the stem and root tissues (Figure 7). The KEGG results for stem upregulated DEGs identified 82 pathways and 25 sequences linked to specific pathways. The pathways with the most sequences from upregulated stem DEGs were starch and sucrose metabolism (3) (ko00500), NOD-like receptor signaling pathway (2) (ko04621), tryptophan metabolism (2) (ko00980), and protein processing in the endoplasmic reticulum (2) (ko04141). KEGG analysis also identified 96 pathways and 25 sequences in stem downregulated DEGs, with the highest downregulated pathway being the plant circadian rhythm (3) (ko04712) ( Figure 7A). KEGG analysis for roots upregulated DEGs found 73 pathways and 26 sequences and 118 pathways and 53 sequences from downregulated roots. The highest number of sequences associated with upregulated roots were protein processing in the endoplasmic reticulum (10 sequences) (ko04141). In contrast to the downregulated roots, the NOD-like receptor signaling pathway (ko04621) (7), MAPK signaling pathway (7) (ko04010), and plant hormone signal transduction (5) (ko04075) were the highest ( Figure 7B). A complete list of all KEGG results can be found in Supplementary Tables S10-S13.

Discussion
Citrus is among the most economically valuable fruit crops. In 2019, citrus production exceeded 157.9 million tons in over 9.8 million hectares worldwide [1]. The ever-reducing global natural resources (i.e., decreased arable land and water availability) combined with climate change, agricultural labor shortages, and the increased disease pressure of deadly citrus diseases such as HLB demand the adoption of novel citriculture practices. Highdensity citrus plantings, the possible implementation of mechanized citrus harvest, and citrus production under protective structures are some of the most promising trends for sustainable citriculture. The ability to produce dwarf citrus trees is critical for implementing these innovations and maximizing returns on investments [19,[42][43][44][45][46].
Dwarfed sweet orange trees propagated on trifoliate orange rootstock, where reduced vegetative growth is achieved without compromising fruit yield or quality, have been produced using CDVd. CDVd is part of a group of so-called "graft-transmissible dwarfing agents," the use of which has been explored for almost fifty years since it was initially proposed [16,18,[47][48][49][50]. Previously, we observed CDVd significantly reduced sweet orange canopy volume on trifoliate rootstock by reducing vegetative growth [19][20][21]. Long-term field observations indicated that CDVd might be used as a possible tool for high-density plantings of citrus. In addition, key findings on the possible biological mechanism through which CDVd affects specific rootstock scion combinations to reduce tree canopy volume were obtained [14]. Understanding the detailed molecular mechanisms leading to reduced canopy volume of commercial citrus tree varieties in response to CDVd infection is critical. Indeed, it would provide the information to produce reduced-size citrus trees without a graft-transmissible viroid agent.
Transcriptomic profiling provides a rapid and cost-effective approach to identifying target genes responsible for the observed reduced vegetative growth of CDVd-infected sweet orange trees. Commercial field citrus trees consist of two different species: the scion, producing the canopy of the tree (i.e., branches, stems, and leaves), and the rootstock, from which the root system of the tree (i.e., taproot, lateral, and feeder roots) develops [51]; therefore, any observed tree phenotype is a result of gene expression not in one but two species that closely interact with each other. In this study, we performed a transcriptome analysis of CDVd-infected stems and roots from sweet orange on trifoliate rootstock trees grown in high-density plantings and compared them with non-infected controls. Here, we showed that CDVd infection affects a wide range of biological functions via different stem and root mRNA transcripts. Furthermore, it complements recent molecular studies targeting CDVd-derived small RNAs on dwarfed citrus [32]. As a result of differential expression analysis, a total of 409 CDVd-dependent differentially expressed genes (DEGs) were deregulated (131 upregulated and 278 downregulated), with the majority of upregulated genes in the stems (83/131), and most downregulated genes in the roots (186/278) ( Table 2).
The five most significant upregulated genes in CDVd-infected stems include: LOC102607140, a probable 28S rRNA (cytosine-C(5))-methyltransferase, an ortholog of the Arabidopsis thaliana (AT5G26180), which plays a role in regulating translation in response to stress [52]; LOC102611719, a germin-like protein subfamily T member. Germin-like proteins (GLPs) are encoded by multigene families in several plant species, and some subfamily members play a role in defense against pathogen attack [53,54]; LOC102610118, a pentatricopeptide repeatcontaining protein (ortholog of AT2g17525). Pentatricopeptide repeat (PPR) proteins belong to a large gene family in plants. Some PPR proteins play important roles in organellar RNA metabolism and organ development in Arabidopsis and rice, while their functions remain unknown in woody species [55]; LOC102616051, cytosolic sulfotransferase 15-like, an ortholog of the Arabidopsis cytosolic sulfotransferase genes. In Arabidopsis, a hydroxyjasmonate sulfotransferase was reported to participate in the hydroxylation and sulfonation reactions that might be components of a pathway that inactivates excess jasmonic acid in plants [56]. The LOC102614284 is associated with aspartic proteinase Asp1-like, and plant aspartic proteinases have been implicated in diverse cellular functions, including protein processing and/or degradation, plant senescence, stress responses, programmed cell death, and reproduction [57].
The five most significantly downregulated genes in CDVd-infected stems include: LOC102627502, a pentatricopeptide repeat-containing protein, an ortholog of At4g21300; LOC102620512, SKP1-like protein 21, which is involved in ubiquitination and subsequent proteasomal degradation of target proteins [58]; LOC102610900, SHORT-ROOT-like protein, which controls radial patterning of the Arabidopsis root, and therefore root growth [59], which is consistent with reduced root mass in viroid infected plants [9]; LOC102620992, an amino-acid permease BAT1-like, which plays a role in primary carbon metabolism and plant growth by mediating the transport of gamma amino butyric acid from the cytosol to mitochondria in Arabidopsis [60]; LOC102628393, serine/threonine-protein kinase ATG1a, involved in autophagy in a nutritional condition-dependent manner [61].
In CDVd-infected roots, the five most significantly upregulated genes include: Ptrif.0009s1914.1 (ortholog of AT5G62350/LOC_Os06g49760), a plant invertase/pectin methylesterase inhibitor superfamily protein, which inhibits pectin methylesterases (PMEs) and invertases (PMIs) and has been implicated in the regulation of fruit development, carbohydrate metabolism, and cell wall extension as well as inhibition of microbial pathogen PMEs. The interplay between PME and PMEI is a determinant of cell adhesion, cell wall porosity, and elasticity, as well as a source of signaling molecules released upon cell wall stress [62]. The Ptrif.0004s1154.1 (ortholog of AT5G07310/Cre14.g620500/LOC_Os04g32620) is an integrase-type DNA-binding superfamily protein involved in the regulation of transcription and recently reported to integrate the jasmonate and cytokinin signaling pathways to repress adventitious rooting in Arabidopsis [63]; Ptrif.0005s1389.1, an ortholog of Cre05.g240400/LOC_Os12g25450, S-adenosyl-L-methionine-dependent O-methyltransferase/ ethylene-responsive transcription factor, which methylates proteins, small molecules, lipids, and nucleic acids [64]; Ptrif.0005s2465.1 (ortholog of AT5G52010/LOC_Os06g07020) C2H2like zinc finger protein, which contains one of the most common domains found in the transcription factors of higher eukaryotes [65]; and Ptrif.0006s1595.1 (ortholog of AT5G03800) another pentatricopeptide repeat (PPR) superfamily protein.
Among the additional interesting CDVd-responsive stem genes whose differential expression was verified by RT-qPCR ( Figure 5A), several were found including the plant defense receptor-like protein 9DC3 (LOC102630240) [70]; the MADS-box transcription factor 23-like (LOC102630220) [71], which is part of the transcription factor (TF) families (WRKY, MADS-box and MYB) that activate unique abiotic and biotic stress-responsive strategies considered as key determinants for defense and developmental processes in most eukaryotic plants [72]; the phospholipase A1-Ibeta2 (LOC102609966) [73], the chloroplast import protein TIC 55 (LOC102614849), the chloroplast import protein [74]; the chloroplast pentatricopeptide repeat-containing protein ortholog of At1g08070 (LOC102619843), which is involved in RNA editing events within the chloroplast [75]; the transcription factor MYB13 (LOC102616893), which regulates meristem function by being a component of a regulatory network controlling the establishment and/or development of the shoot system in Arabidopsis [76]; protein IQ-DOMAIN 1, transcript variant X1 (LOC102628194) which plays a role with calmodulins or calmodulin-like proteins as well as involved in scaffolding in cellular signaling and trafficking [77,78]; hydroxyproline O-galactosyltransferase GALT3, transcript variant X3 (LOC102626136), which commits arabinogalactan proteins to the first step in arabinogalactan polysaccharide addition. AGP glycans play essential roles in both vegetative and reproductive plant growth [79]; cytochrome P450 82C4-like, which is involved in the early Fe deficiency response in Arabidopsis [80]; LOC102622795; thaumatinlike protein 1 (LOC102619709), involved in local responses of roots to colonization by non-pathogenic plant growth-promoting rhizobacteria [81].
Gene ontology (GO) assignments were used to classify the function of the identified DEGs. GO term enrichment analysis showed that in both stems and roots of CDVd-infected trees, for the biological process category, the GO terms corresponding to metabolic process and cellular process were the most enriched terms in response to CDVd; for molecular function, the catalytic activity and binding are the most enriched, and for the cellular component category, the cellular anatomical entity term was most enriched in the stems ( Figure 6A,B).
In order to gain further insight into the function of DEGs correlating to a reduced vegetative growth as observed in CDVd-infected trees, the KEGG pathway classification was performed. Overall, circadian rhythms (ko04712), starch and sucrose metabolism (ko00500), and protein processing in the endoplasmic reticulum (ko04141) were most altered in response to CDVd infection. Circadian rhythms (ko04712), downregulated in CDVdinfected stems, are known to cross-talk with defense signaling against bacteria in plants [85] and to be important determinants in the outcome of plant-pathogen interactions [86]. Starch and sucrose metabolism (ko00500) were upregulated in both CDVd-infected stems and roots. Sucrose, the main form of assimilated carbon produced during photosynthesis, has important roles as a signaling molecule, and it is involved in many metabolic processes in plants; it is essential for plant growth and development and plays a role in plant defense by activating plant immune responses against pathogens. Upon infection, pathogens hijack the plant metabolism to access the plant sugars and, in doing so, trigger plant defense responses. Invertases appear to be involved in establishing plant defense responses [87].
Protein processing in the endoplasmic reticulum (ko04141) was upregulated in roots and stems, consistent with higher metabolic activity levels. Tryptophan metabolism (ko00380) was up-and downregulated in stems. Tyrosine and tryptophan are precursors for the plant defense compounds dhurrin and indole glucosinolates, respectively. In addition, tryptophan is a precursor for the essential phytohormone indole-3-acetic acid [88]. Other KEGG pathway classification terms that were altered in response to CDVd infection include: Lipid metabolism (glycerophospholipid (ko00564), which was downregulated in roots and stems; ether lipid (ko00565) was downregulated in stems; ABC transporters (ko 02010), Ras signaling pathway, mTOR signaling pathway (ko04150), which were downregulated in stems.
To date, only one previous citrus gene expression profiling study in response to CDVd infection has been reported [22], and it differs from ours in many important aspects. Indeed, it predates the widespread use of HTS applications and the availability of reference genomes for C. sinensis and P. trifoliata. Instead, the authors profiled the transcriptome using the differential display technique (DDRT-PCR). Most importantly, however, the authors studied the sensitive viroid bioindicator 'Etrog' citron, Arizona 861-S1, newly inoculated with CDVd (18 months old) grown under growth chamber conditions and expressing symptoms of leaf drooping and petiole necrosis. Provided these experimental conditions, as expected, among the upregulated genes, a suppressor of RNA silencing was identified [22]. In contrast, our study used HTS technologies to analyze 18-year-old, field-grown commercial citrus tree varieties that do not express typical viroid symptoms of leaf epinasty and midvein and petiole necrosis [89]. Therefore, it is not surprising that the transcriptome profiling was different.
When taken together, our results indicate that CDVd modulates the expression profile of important citrus growth and developmental processes that may participate in the cellular changes leading to the observed phenotype of reduced vegetative growth and overall smaller tree size. In this study, most of the transcriptome expression reprogramming appeared to occur in the stems (C. sinensis) rather than in the roots (C. trifoliata), which is consistent with the phenotypic observation of the reduced stem vegetative growth of the stems and striking dwarfed citrus tree phenotype in CDVd-infected trees [21]. It is also in agreement with the observation that the trifoliate orange rootstock does not display major symptoms in response to CDVd infection [90][91][92]. Even though our study demonstrates the potential use of modern molecular technologies to decipher complex plant responses to biotic factors using plants growing in agricultural production systems, future studies employing transgenics and model plant systems could provide additional evidence to dissect further the genes and pathways participating in the development of the dwarf phenotype identified in this study.
Lastly, the lack of major alterations in the pathogen defense transcriptome profile in our study is in contrast to studies with other viroids that demonstrated altered expression of plant defense genes [11,24]. This finding supports the idea that when viroid RNAs such as CDVd are used for horticultural purposes, such as reducing tree canopy volume, perhaps they are better described as "Transmissible small nuclear ribonucleic acids" (TsnRNAs). This term provides a more detailed description of the TsnRNA's hallmark properties (i.e., transmissibility, small genome size, site of replication, and RNA nature) as well as dissociation from a generic virology term, such as "viroid", which implies disease and crop loss. Therefore TsnRNAs that do not express a disease syndrome but rather act as modifying agents of tree performance resulting in desirable agronomic traits, such as reduced tree size for high-density plantings with potential economic and sustainability advantages, can be distinguished from the pathogenic viroids [18,20].