Skip to main content

In silico analysis of the grapefruit sRNAome, transcriptome and gene regulation in response to CTV-CDVd co-infection



Small RNA (sRNA) associated gene regulation has been shown to play a significant role during plant-pathogen interaction. In commercial citrus orchards co-infection of Citrus tristeza virus (CTV) and viroids occur naturally.


A next-generation sequencing-based approach was used to study the sRNA and transcriptional response in grapefruit to the co-infection of CTV and Citrus dwarfing viroid.


The co-infection resulted in a difference in the expression of a number of sRNA species when comparing healthy and infected plants; the majority of these were derived from transcripts processed in a phased manner. Several RNA transcripts were also differentially expressed, including transcripts derived from two genes, predicted to be under the regulation of sRNAs. These genes are involved in plant hormone systems; one in the abscisic acid, and the other in the cytokinin regulatory pathway. Additional analysis of virus- and viroid-derived small-interfering RNAs (siRNAs) showed areas on the pathogen genomes associated with increased siRNA synthesis. Most interestingly, the starting position of the p23 silencing suppressor’s sub-genomic RNA generated a siRNA hotspot on the CTV genome.


This study showed the involvement of various genes, as well as endogenous and exogenous RNA-derived sRNA species in the plant-defence response. The results highlighted the role of sRNA-directed plant hormone regulation during biotic stress, as well as a counter-response of plants to virus suppressors of RNA-silencing.


Plants respond to pathogen infection through a number of gene regulatory pathways. RNA-silencing is a form of regulation where double-stranded or hairpin-structured RNA precursors give rise to small RNAs (sRNAs), which become control elements for the expression of target genes [1,2,3]. Several types of sRNA species have been identified and characterised in plants, and their involvement in biotic stress responses have been suggested. These include microRNAs (miRNAs) [4,5,6,7,8], phased-siRNAs (phasiRNAs) [9,10,11,12], natural-antisense transcript siRNAs (nat-siRNAs) [12,13,14], repeat-associated siRNAs (rasiRNAs) [7] and tRNA-derived RNA fragments (tRFs) [15, 16].

Stem pitting is a destructive symptom in fruit crops caused by various virus species. Citrus tristeza virus (CTV) is a highly destructive, phloem limited, pathogen of citrus species, causing three different disease syndromes [17] of which stem pitting currently is considered the highest threat to the industry [18]. CTV belongs to the genus Closterovirus in the family Closteroviridae [19]. The ~19,300 nt genome of CTV is organised into 12 open reading frames and represents the largest known plant virus genome [20].

Due to the complex disease aetiology of CTV, which is strongly influenced by the combination of host factors and virus genotypes, the mechanism(s) behind symptom expression is poorly understood [21, 22]. A recent study, however, has suggested the involvement of different combinations of p33, p18 and p13 expression during stem pitting symptom development [18]. Three CTV RNA-silencing suppressors have been identified, namely p20, p23 and p25 (coat protein), all of which can play a role in symptom expression [23].

In addition to viruses, citrus species are also affected by viroid-infections. Citrus dwarfing viroid (CDVd), a member of the genus Apscaviroid (family Pospiviroidae), has been suggested for use in high-density orchards since it causes dwarfing of citrus varieties grafted onto Poncirus trifoliata (P. trifoliata) and its hybrids, without reducing production [24,25,26]. The fact that citrus species are often co-infected with CTV and viroids prompted a recent study which investigated the co-infection of CTV and CDVd in their respective indicator plants, Mexican lime and etrog citron [27]. A host-specific increase in the accumulation of CDVd was observed in the presence of CTV, along with the synthesis of CDVd-associated sRNAs. These observations were ascribed mainly to the involvement of the CTV-encoded silencing suppressor, p23. It is also interesting to note that the co-infection did not affect symptom expression under their experimental conditions [27].

Understanding the mechanisms involved in pathogen infection and symptom expression provides the information required to study and potentially engineer disease resistance. In this study, we used a next-generation sequencing (NGS) approach to investigate the plant responses to the co-infection of CTV and CDVd in two commercial grapefruit (Citrus paradisi) cultivars, on both the sRNA and transcriptome level. Our results highlighted the association of CTV-CDVd co-infection with the expression of various genes and sRNA species, these include sRNAs involved in the regulation of plant hormones.


Sample preparation

Grapefruit plants (cultivars ‘Marsh’ and ‘Star Ruby’) on ‘Carrizo’ citrange rootstocks were bark-inoculated from an asymptomatic Citrus sinensis (sweet orange) plant that was confirmed to be infected with CTV (genotype T3) and CDVd by RT-PCR using previously described protocols [28, 29]. Briefly, bark-inoculation was performed by patch-grafting two bark chips of the source plant to each grapefruit scion, once the scion was approximately 7 mm thick. All plants were inoculated at the same height. After inoculation, the scions were cut back approximately 10 cm above the inoculation point. One shoot of the new growth was allowed to grow from the top bud. Un-inoculated plants served as healthy controls.

Total RNA was extracted from the phloem material of three replicates of healthy and infected plants of each cultivar following an adapted CTAB method from Carra et al. [30]. Virus and viroid status was confirmed using the above-mentioned RT-PCR assays.

Next-generation sequencing and data preparation

Total RNA extracted from each sample was sent for sequencing on an Illumina HiSeq instrument (Fasteris, Geneva, Switzerland). Two libraries per sample were prepared and sequenced. An sRNA library was generated from 18 to 30 nt size-selected RNA and sequenced in a 1 × 50 nt run, as well as a transcriptome library generated from ribo-depleted total RNA and sequenced in a 2 × 125 nt run. Adapter sequences were trimmed from the data using cutadapt [31]. Fastx-toolkit [32] was used to remove all low quality reads from the sRNA data, while Trimmomatic [33] was used to filter and trim the transcriptome data for quality. sRNA reads, 18–26 nts in length and transcriptome reads, 20 nts and longer were retained for further analyses.

Virus and viroid infection status of samples were confirmed bioinformatically with the mapping of virus-derived siRNA (vsiRNA) and viroid-derived siRNA (vd-siRNA) derived NGS data against the respective genomes as described below. BLASTn [34] analysis of assembled contigs (described below) against NCBI’s nt database was used to exclude the possible presence of any other viruses or viroids from the data. The viral status of samples were further verified using CTV-specific e-probes [35, 36].

Grapefruit transcriptome-assembly, differential expression analysis and natural-antisense transcript (NAT) identification

Trinity [37, 38] was used to assemble the transcriptome data into contigs, applying default parameters. Transcript differential expression analysis was performed using the DESeq2 [39] method in Trinity. Gene ontology analyses were performed using Trinity and Blast2GO [40].

Assembled contigs were used to identify NATs by aligning contigs to each other, using BLASTn [34] to identify overlapping regions of 50 nts and longer with 100% identity. Duplex formation of the overlaps was validated with UNAfold [41].

Grapefruit sRNA identification, differential expression analysis and target prediction

To identify novel miRNAs and phased (PHAS) transcripts within grapefruit, all sRNA datasets were simultaneously submitted to ShortStack [42, 43], with default parameters. The assembled contigs served as template for precursor identification. DESeq2 was used to analyse the differential expression of all unique sRNA sequences. Differentially expressed sRNAs were characterised based on their comparison to different sRNA species as described below.

The predicted novel miRNAs, along with the plant entries in miRBase 21 [44,45,46,47], as well as the predicted phasiRNAs (sRNAs that fell into a dominant phasing register), served as database for the identification of differentially expressed sRNAs, which were miRNAs or phasiRNAs, respectively. Differentially expressed sRNAs were also mapped, using Bowtie [48], onto the identified PHAS transcripts, the overlapping regions of the NATs, plant repetitive sequences in RepBase [49, 50], and the sequences of plant tRNAs in the PlantRNA database [51], to identify PHAS-associated sRNAs (not in phase with the dominant phasing register), nat-siRNAs, rasiRNAs and tRFs, respectively.

psRNATarget [52] was used, applying default settings, to predict targets for the differentially expressed sRNAs using the assembled transcriptome as a list of potential targets.

Pathogen-derived sRNA analysis

vsiRNAs (associated with CTV) and vd-siRNAs (associated with CDVd) were identified by mapping the sRNA reads, using Bowtie, onto the CTV-T3 (Accession No. KC525952) and CDVd (Accession No. AF184149) genomes allowing a single, or no mismatches, respectively.


Symptom expression in grapefruit

Healthy ‘Marsh’ and ‘Star Ruby’ grapefruit plants were co-infected with CTV and CDVd, using an asymptomatic (CTV and CDVd infection confirmed) sweet orange plant as source. The co-infection was confirmed with RT-PCRs and supported through NGS read mapping analysis (shown below). No additional viruses or viriods were identified in the NGS data. ‘Star Ruby’ plants showed more distinct leaf cupping and stem pitting symptoms than ‘Marsh’ plants (Fig. 1).

Fig. 1
figure 1

CTV-CDVd co-infected plants. Three healthy followed by three CTV-CDVd co-infected (a) ‘Marsh’ and (b) ‘Star Ruby’ plants. A representative stem sample (after bark removal) is given below each plant (c and d)

NGS data preparation

sRNA and transcriptome NGS datasets were generated for each RNA sample extracted from phloem tissue. The raw data ranged from 14,740,885 to 22,862,616 sRNA reads and 12,726,094 to 16,410,600 transcriptome read-pairs per sample, while the high-quality datasets ranged from 8,550,133 to 12,715,167 sRNA reads (after quality filtering) and 12,012,340 to 15,458,590 transcriptome read-pairs (after quality filtering and trimming) per sample.

Grapefruit transcriptome assembly and differential expression

High-quality reads from all the transcriptome datasets were combined and de novo assembled into transcripts (contigs). Altogether 214,371 transcripts were generated, which could be grouped into 120,991, Trinity-defined, “genes”.

Differential expression analysis was subsequently performed to identify transcripts involved in CTV-CDVd co-infection. In ‘Marsh’ and ‘Star Ruby’ 675 and 1204 transcripts, respectively, showed altered expression (Additional file 1: Tables S1 and S2). The results also identified 154 “genes” for which at least one transcript were differentially expressed between healthy and infected plants, across both grapefruit varieties (Additional file 1: Table S3). According to similarity searches, these included 21 potential disease response genes, as well as 60 membrane and 10 photosystem-associated genes, highlighted through gene ontology (GO) analysis (Fig. 2, Additional file 1: Tables S4-S6 and Additional file 2: Figures S1-S3). Many transcripts were however homologous to hypothetical or uncharacterised citrus proteins, resulting in 33% of the differentially expressed “genes” remaining unidentified.

Fig. 2
figure 2

Gene ontology classification of differentially expressed genes. Barr-graph illustrating the number of differentially expressed genes assigned to biological process, molecular function and cellular component gene ontology terms (level 3)

Endogenous sRNA identification and regulation

Combined analysis of the transcript and sRNA data predicted miRNAs from 60 grapefruit miRNA genes (MIRs) that were expressed in at least one of the samples (Additional file 3: Table S7). For 38 of these MIRs, neither of the mature miRNAs predicted, were represented by any mature plant-derived miRNA in miRBase. In addition to the predicted miRNAs, reads with sequences identical to 216 known plant miRNAs were also present in the data (Additional file 3: Table S8). Significant phasing was also seen in 7268 transcripts (called PHAS transcripts), producing 63,943 phasiRNAs in total, which were in phase with the dominant phasing register. To facilitate the identification of nat-siRNAs, transcripts were subjected to NAT identification. The duplex formation of 25,378 transcripts, predicted to be part of one or more NAT pair, were validated in silico. The overlapping regions were extracted for subsequent nat-siRNA analysis.

Differential expression analysis revealed 761 sRNAs with altered expression levels resulting from pathogen infection (Additional file 4: Tables S9 and S10). Of these, 577 were variety-specific, while 184 showed differential expression across varieties (Additional file 4: Table S11). These sRNAs were characterised based on sequence homology to either the predicted, or other plant miRNAs, predicted phasiRNAs, PHAS transcripts, the overlaps of NATs, as well as plant repetitive DNA-regions and tRNAs. While a number of sRNAs could be classified as nat-siRNAs (22), miRNAs (five), rasiRNAs (17) or tRFs, the majority (59) of differentially expressed sRNAs were phasiRNAs. Some could potentially be classified into more than one sRNAs species, for example two sRNAs derived from a repetitive region and seven sRNAs derived from the overlapping region of a NAT that were all processed in a phased manner (Fig. 3).

Fig. 3
figure 3

Species identification of differentially expressed sRNAs. Venn diagram illustrating the overlapping sRNA species (miRNA, phasiRNA, nat-siRNA, rasiRNA, tRF) identities of differentially expressed sRNAs

To determine the biological role that sRNAs play during CTV-CDVd co-infection, all differentially expressed sRNAs were subjected to in silico target prediction. Only three sRNAs (one uncharacterised and two derived from PHAS transcripts) showed significant inverse-regulation with respect to their predicted target transcripts (Additional file 5: Table S12), one of which was across the two varieties, and the other two specific to ‘Star Ruby’. Homology searches identified the across-variety sRNA target as a chloroplastic Magnesium-chelatase subunit ChlH (CHLH) and the ‘Star Ruby’-specific sRNA targets as Cytokinin dehydrogenase 6 (CKX) and a hypothetical protein.

Pathogen-derived sRNAs

Virus-derived siRNAs (vsiRNAs) were found associated with 97% of the CTV-T3 reference genome. The majority of the vsiRNAs were 21 or 22 nts in length (Fig. 4). The distribution of the sRNA reads on the genome showed an increase in vsiRNAs mapping towards the 3′ end of the virus (Fig. 5). A prominent hotspot for sRNA synthesis was observed on the negative-strand at the sub-genomic RNA initiation site of p23.

Fig. 4
figure 4

Size-distribution of vsiRNA and vd-siRNA reads. Histogram illustrating the number of vsiRNA and vd-siRNA reads, 18 nt to 26 nt in length, from the infected samples, all (redundant) as opposed to unique (non-redundant, NR), as a percentage of the vsiRNA and vd-siRNA reads in this size-range respectively

Fig. 5
figure 5

Distribution of vsiRNA reads along the CTV genome. vsiRNA profile generated from sRNA reads depicted as heat maps showing the reads that mapped onto the positive (+) or negative (−) strand of CTV. A schematic representation of the genome above the heat maps illustrates the genomic position of the vsiRNA reads. The start of the p23 subgenomic RNA, which forms an vsiRNA hotspot, is indicated with an arrow

In the case of CDVd, viroid-derived siRNAs (vd-siRNAs) were mostly 22 nts in length, followed by reads 21 and 24 nts in length (Fig. 4), and covered the complete viroid genome (Fig. 6). A specific area on the negative-strand of the viroid, overlapping the central and variable regions, gave rise to an abundance of sRNAs, indicating a potential target area.

Fig. 6
figure 6

Distribution of vd-siRNA reads along the CDVd genome. vd-siRNA profile generated from sRNA reads depicted as heat maps showing the reads that mapped onto the positive (+) or negative (−) strand of CDVd. A vd-siRNA hotspot was formed on the negative-strand of the viroid, overlapping the central and variable regions. C, central domain; P, pathogenic domain; TCR, terminal conserved region; CCR, central conserved region; TL, terminal left domain; TR, terminal right domain; V, variable domain


In the field, citrus plants are often subjected not only to CTV infection but also to the co-infection with viroids. It is therefore necessary to understand this combined pathogen interaction with the plant host in order to improve disease resistance strategy design. In this study, healthy ‘Marsh’ and ‘Star Ruby’ grapefruit plants developed leaf cupping and stem pitting symptoms after being co-infected with CTV (genotype T3) and CDVd. These symptoms were not present in the co-infected sweet orange plants, which served as inoculation source. Although the T3 genotype is associated with increased stem pitting, previous studies have shown that CTV isolate is not the only determinant of symptom expression, but that host species also plays a role [18, 21, 22]. The mechanism(s) that drive the severity and host-specific symptom expression remains to be elucidated.

An sRNA and transcriptome next-generation sequencing approach was followed to study the gene-regulatory pathways involved in the CTV-CDVd co-infection of grapefruit. To compensate for the limited genomic information available for grapefruit, the transcriptome data were de novo assembled to generate a case-specific grapefruit transcriptome. Since the assembled transcripts represented both coding and non-coding transcripts, they were used to identify both the precursors and targets of sRNAs.

As potential precursor source, the transcripts were first used for miRNA discovery. The miRNA registry, miRBase, currently holds no entries for grapefruit. Here we report on the identification of 60 grapefruit MIR genes along with their mature miRNAs, based on in silico prediction analysis. Many of the other sRNA reads represented homologous plant miRNAs. Once more grapefruit genome information becomes available, these homologous sequences may still prove to be true miRNAs, expressed from grapefruit MIR genes. In addition to miRNAs, phasiRNAs and PHAS transcripts were also identified, based on the assembled transcriptome, along with NATs that form the precursors of nat-siRNAs.

Many diverse “genes” were found differentially regulated in response CTV-CDVd co-infection across both grapefruit varieties. Similar to the results of previous studies on the infection of either CTV [53,54,55,56,57] or CDVd [58] in different citrus species and in P. trifoliata, grapefruit metabolic, disease response, structural and phytohormone-related pathways seem to be affected by the co-infection. Genes with known involvement in plant disease responses were mostly grouped into leucine-rich repeat (which include tobacco mosaic virus resistance protein N-like) and ankyrin repeat-containing protein coding genes. The number of chloroplast-associated genes with altered expression supports the hypothesis that the CTV-CDVd co-infection influences the photosynthetic pathways of grapefruit, which was previously shown for CTV infection in sweet orange [56, 57] and Mexican lime [54, 59].

Many members of different endogenous sRNA species, such as miRNAs, phasiRNAs, nat-siRNAs, rasiRNAs and tRFs also showed variation in expression resulting from the co-infection. Target prediction was performed to determine the biological role of the differentially expressed sRNAs. Despite many transcripts and sRNAs showing differential expression resulting from co-infection, an inverse-regulation was only seen for three sRNA-transcript pairs. The apparent disconnect between sRNAs and their predicted targets could be due to a number of factors. First, the target prediction software used, was designed specifically for miRNA and phasiRNA target prediction and may therefore not be as effective for other sRNA species. Target prediction models remain to be developed for the other species, following the characterisation of their mechanism(s) of action. Second, although sRNA action may lead to the indirect inhibition of protein expression, the presence of any (cleaved or uncleaved) target transcript-related reads in the transcriptome dataset will count towards transcript-associated reads during differential expression analysis of the genes. Last, the relatively low read counts associated with many of the predicted target transcripts could have influenced the statistical significance of the differential expression analysis.

The combined results from this study suggested the sRNA-directed gene regulation of plant hormone pathways during CTV-CDVd co-infection in grapefruit. The expression of chloroplastic Magnesium-chelatase subunit ChlH (CHLH) showed inverse-regulation with respect to that of its regulating sRNA, across both ‘Marsh’ and ‘Star Ruby’ plants. In Arabidopsis, CHLH plays a role in chlorophyll biosynthesis, the expression of photosynthesis-related proteins, as well as abscisic acid (ABA) signal regulation [60]. CHLH was shown to repress the expression of the disease response WRKY40 gene in the ABA signalling pathway [61]. Therefore, unsurprisingly, WRKY40 showed inverse-regulation to that of CHLH resulting from the co-infection. The involvement of the ABA pathway in plant-virus response has previously been described [62,63,64], and includes the restriction, to some extent, of virus movement through callose deposition [62]. In addition, ABA was shown to contribute to virus-resistance through the regulation of Argonauts [64]. The sRNA-directed regulation of CHLH during the CTV-CDVd co-infection could therefore potentially have a down-stream effect on virus resistance through ABA regulation. Cytokinin dehydrogenase 6 (CKX), on the other hand, is involved in the irreversible down-regulation of cytokinins [65,66,67], and showed inverse-regulation with respect to that of its regulating sRNA in ‘Star Ruby’ plants. Cytokinins can stimulate plant defence response upon pathogen infection [68], and may lead to either an enhancement in resistance [69,70,71], or susceptibility to viral infections [72]. The observed down-regulation of CKX may lead to an increase in cytokinin levels, contributing to the grapefruit defence response. A tissue-specific up- or down-regulation of CKX resulting from virus infection was also recently observed in Arabidopsis roots and shoots, respectively [73].

Pathogen-derived sRNA profiles have been found to vary upon different infection. The majority of the vsiRNAs, derived from CTV, were 21 or 22 nts in length, which is common for vsiRNAs produced by many plants, including citrus [8, 74]. The vsiRNAs also seem to favour the 3′ end of the virus. While this increase towards the 3′ end of CTV was observed before, not all genotype-host combinations showed the same pattern [8, 74]. Factors that could influence patterns of vsiRNA synthesis include; a host-specific response, the virus genome secondary structure, virus gene expression and host-specific suppression of virus proteins.

Interestingly, when looking at the vsiRNA profile of CTV, a highly prominent hotspot was observed at the sub-genomic RNA initiation site of p23, associated with the negative strand, which indicates an area targeted by the host-response. p23 is known to play a major role in CTV-host interaction, especially as a suppressor to counter the host’s RNA-silencing response [23]. The targeting of p23 by the host therefore adds another layer to the CTV-host interaction.

Recent studies have also investigated the plant siRNA response to viroid infection [75,76,77,78,79]. While the majority of CDVd-derived vd-siRNAs were 22, 21 or 24 nts in length respectively, the observed size-distribution may be tissue dependent [76, 77]. As was seen for CTV, an sRNA hotspot was observed associated with an area on the negative-strand of CDVd. A similar hotspot was previously observed for another member of the same genus, Potato spindle tuber viroid, in tomato [75, 78]. The implication of the sRNA targeting of this specific area of the viroid remains to be elucidated.


In silico analysis of the transcriptome and sRNAs generated in response to CTV and CDVd co-infection in grapefruit, suggested the involvement of sRNAs in the regulation of plant hormone pathways in infected plants. Further analysis revealed regions within both the CTV and CDVd genomes that form hotspots for vsiRNA and vd-siRNA synthesis. The specific vsiRNA hotspot associated with p23, suggests a first report of potential vsiRNA counter-response by the plant to the p23-silencing suppressor of CTV. An exciting future prospect for these pathogen-derived sRNAs could be their application in disease resistance strategies.



abscisic acid


Citrus dwarfing viroid


chloroplastic magnesium-chelatase subunit ChlH


cytokinin dehydrogenase 6


Citrus tristeza virus


miRNA gene




natural-antisense transcript


natural-antisense transcript siRNA



P. trifoliata :

Poncirus trifoliata


phased transcript




repeat-associated siRNA


small-interfering RNA


small RNA


tRNA-derived RNA fragment


viroid-derived siRNA


virus-derived siRNA


  1. Hamilton AJ, Baulcombe DC. A species of small antisense RNA in posttranscriptional gene silencing in plants. Science. 1999;286:950–2.

    Article  CAS  PubMed  Google Scholar 

  2. Mette MF, Aufsatz W, van der Winden J, Matzke MA, Matzke AJ. Transcriptional silencing and promoter methylation triggered by double-stranded RNA. EMBO J. 2000;19:5194–201.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Brodersen P, Sakvarelidze-Achard L, Bruun-Rasmussen M, Dunoyer P, Yamamoto YY, Sieburth L, et al. Widespread translational inhibition by plant miRNAs and siRNAs. Science. 2008;320:1185–90.

    Article  CAS  PubMed  Google Scholar 

  4. Bester R, Burger JT, Maree HJ. Differential expression of miRNAs and associated gene targets in grapevine leafroll-associated virus 3-infected plants. Arch Virol. 2017;162:987–96.

    Article  CAS  PubMed  Google Scholar 

  5. Gao R, Wan ZY, Wong S-M. Plant growth retardation and conserved miRNAs are correlated to hibiscus chlorotic ringspot virus infection. PLoS One. 2013;8:e85476.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Naqvi AR, Haq QM, Mukherjee SK. MicroRNA profiling of tomato leaf curl new delhi virus (ToLCNDV) infected tomato leaves indicates that deregulation of mir159/319 and mir172 might be linked with leaf curl disease. Virol J. 2010;7:281.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Romanel E, Silva TF, Corrêa RL, Farinelli L, Hawkins JS, Schrago CEG, et al. Global alteration of microRNAs and transposon-derived small RNAs in cotton (Gossypium Hirsutum) during cotton leafroll dwarf polerovirus (CLRDV) infection. Plant Mol Biol. 2012;80:443–60.

    Article  CAS  PubMed  Google Scholar 

  8. Ruiz-Ruiz S, Navarro B, Gisel A, Peña L, Navarro L, Moreno P, et al. Citrus tristeza virus infection induces the accumulation of viral small RNAs (21–24-nt) mapping preferentially at the 3′-terminal region of the genomic RNA and affects the host small RNA profile. Plant Mol Biol. 2011;75:607–19.

    Article  CAS  PubMed  Google Scholar 

  9. Du P, Wu J, Zhang J, Zhao S, Zheng H, Gao G, et al. Viral infection induces expression of novel phased microRNAs from conserved cellular microRNA precursors. PLoS Pathog. 2011;7:e1002176.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Li F, Pignatta D, Bendix C, Brunkard JO, Cohn MM, Tung J, et al. MicroRNA regulation of plant innate immune receptors. Proc Natl Acad Sci. 2012;109:1790–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Shivaprasad PV, Chen H-M, Patel K, Bond DM, Santos BACM, Baulcombe DC. A microRNA superfamily regulates nucleotide binding site-leucine-rich repeats and other mRNAs. Plant Cell. 2012;24:859–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Quintero A, Pérez-Quintero AL, López C. Identification of ta-siRNAs and cis-nat-siRNAs in cassava and their roles in response to cassava bacterial blight. Genomics Proteomics Bioinformatics. 2013;11:172–81.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Katiyar-Agarwal S, Morgan R, Dahlbeck D, Borsani O, Villegas A, Zhu J-K, et al. A pathogen-inducible endogenous siRNA in plant immunity. Proc Natl Acad Sci. 2006;103:18002–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Katiyar-Agarwal S, Gao S, Vivian-Smith A, Jin H. A novel class of bacteria-induced small RNAs in Arabidopsis. Genes Dev. 2007;21:3123–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Loss-Morais G, Waterhouse PM, Margis R. Description of plant tRNA-derived RNA fragments (tRFs) associated with argonaute and identification of their putative targets. Biol Direct. 2013;8:6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Visser M, Maree HJ, Rees DJ, Burger JT. High-throughput sequencing reveals small RNAs involved in ASGV infection. BMC Genomics. 2014;15:568.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Dawson WO, Garnsey SM, Tatineni S, Folimonova SY, Harper SJ, Gowda S. Citrus tristeza virus-host interactions. Front Microbiol. 2013;4 doi:10.3389/fmicb.2013.00088.

  18. Tatineni S, Dawson WO. Enhancement or attenuation of disease by deletion of genes from citrus tristeza virus. J Virol. 2012;86:7850–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Bar-Joseph M, Garnsey SM, Gonsalves D. The Closteroviruses: a distinct group of elongated plant viruses. Adv Virus Res. 1979;25:93–168.

    Article  CAS  PubMed  Google Scholar 

  20. Karasev AV, Boyko VP, Gowda S, Nikolaeva OV, Hilf ME, Koonin EV, et al. Complete sequence of the citrus tristeza virus RNA genome. Virology. 1995;208:511–20.

    Article  CAS  PubMed  Google Scholar 

  21. Garnsey SM, Civerolo EL, Gumpf DJ, Paul C, Hilf ME, Lee RF, et al. Biological characterization of an international collection of citrus tristeza virus (CTV) isolates. In: Proceedings of th 16th conference of International Organization of Citrus Virologists; 2005. p. 75–93. Accessed 30 Mar 2017.

    Google Scholar 

  22. Hilf ME, Mavrodieva VA, Garnsey SM. Genetic marker analysis of a global collection of isolates of citrus tristeza virus: characterization and distribution of CTV genotypes and association with symptoms. Phytopathology. 2005;95:909–17.

    Article  CAS  PubMed  Google Scholar 

  23. Lu R, Folimonov A, Shintaku M, Li W-X, Falk BW, Dawson WO, et al. Three distinct suppressors of RNA silencing encoded by a 20-kb viral RNA genome. Proc Natl Acad Sci U S A. 2004;101:15742–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Hutton RJ, Broadbent P, Bevington KB. Viroid dwarfing for high density citrus plantings. In: Janick J, editor. Horticultural Reviews. Oxford: John Wiley & Sons, Inc; 1999. p. 277–317. doi:10.1002/9780470650776.ch6.

  25. Hardy S, Sanderson G, Barkley P, Donovan N. Dwarfing citrus trees using viroids. 704th ed. New South Wales, Australia: NSW Department of Primary Industries; 2007. Accessed 31 Mar 2017

    Google Scholar 

  26. Vidalakis G, Pagliaccia D, Bash JA, Afunian M, Semancik JS. Citrus dwarfing viroid: effects on tree size and scion performance specific to Poncirus trifoliata rootstock for high-density planting. Ann Appl Biol. 2011;158:204–17.

    Article  Google Scholar 

  27. Serra P, Bani Hashemian SM, Fagoaga C, Romero J, Ruiz-Ruiz S, Gorris MT, et al. Virus-viroid interactions: citrus tristeza virus enhances the accumulation of citrus dwarfing viroid in Mexican lime via virus-encoded silencing suppressors. J Virol. 2014;88:1394–7.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Cook G, van Vuuren SP, Breytenbach JHJ, Burger JT, Maree HJ. Expanded strain-specific RT-PCR assay for differential detection of currently known citrus tristeza virus strains: a useful screening tool. J Phytopathol. 2016;164:847–51.

    Article  CAS  Google Scholar 

  29. Cook G, van Vuuren SP, Breytenbach JHJ, Steyn C, Burger JT, Maree HJ. Characterization of citrus tristeza virus single-variant sources in grapefruit in greenhouse and field trials. Plant Dis. 2016;100:2251–6.

    Article  Google Scholar 

  30. Carra A, Gambino G, Schubert A. A cetyltrimethylammonium bromide-based method to extract low-molecular-weight RNA from polysaccharide-rich plant tissues. Anal Biochem. 2007;360:318–20.

    Article  CAS  PubMed  Google Scholar 

  31. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjournal. 2011;17:10–2.

    Google Scholar 

  32. FASTX-Toolkit. Accessed 28 Mar 2017.

  33. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinforma Oxf Engl. 2014;30:2114–20.

    Article  CAS  Google Scholar 

  34. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.

    Article  CAS  PubMed  Google Scholar 

  35. Visser M, Burger JT, Maree HJ. Targeted virus detection in next-generation sequencing data using an automated e-probe based approach. Virology. 2016;495:122–8.

    Article  CAS  PubMed  Google Scholar 

  36. Jooste TL, Visser M, Cook G, Burger JT, Maree HJ. In silico probe-based detection of citrus viruses in NGS data. Phytopathol. 2017;107:988–93.

  37. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29:644–52.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with trinity. Nat Protoc. 2013;8:1494–512.

    Article  CAS  PubMed  Google Scholar 

  39. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Conesa A, Gotz S. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genomics. 2008;2008 doi:10.1155/2008/619832.

  41. Markham NR, Zuker M. UNAFold: software for nucleic acid folding and hybridization. Methods Mol Biol Clifton NJ. 2008;453:3–31.

    Article  CAS  Google Scholar 

  42. Axtell MJ. ShortStack: comprehensive annotation and quantification of small RNA genes. RNA. 2013; doi:10.1261/rna.035279.112.

  43. Johnson NR, Yeoh JM, Coruh C, Axtell MJ. Improved placement of multi-mapping small RNAs. G3 Genes Genomes Genet. 2016;6:2103–11.

    Google Scholar 

  44. Griffiths-Jones S, Grocock RJ, Dongen S, van Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(suppl 1):D140–4.

    Article  CAS  PubMed  Google Scholar 

  45. Griffiths-Jones S, Saini HK, Dongen S, van Enright AJ. miRBase: tools for microRNA genomics. Nucleic Acids Res. 2008;36(suppl 1):D154–8.

    CAS  PubMed  Google Scholar 

  46. Kozomara A, Griffiths-Jones S. miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic Acids Res. 2011;39(suppl 1):D152–7.

    Article  CAS  PubMed  Google Scholar 

  47. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res, 2013. 42(Database issue):D68–73.

  48. Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10:R25.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J. Repbase update, a database of eukaryotic repetitive elements. Cytogenet Genome Res. 2005;110:462–7.

    Article  CAS  PubMed  Google Scholar 

  50. Kapitonov VV, Jurka J. A universal classification of eukaryotic transposable elements implemented in Repbase. Nat Rev Genet. 2008;9:411–2.

    Article  PubMed  Google Scholar 

  51. Cognat V, Pawlak G, Duchene A-M, Daujat M, Gigant A, Salinas T, et al. PlantRNA, a database for tRNAs of photosynthetic eukaryotes. Nucleic Acids Res. 2012;41:D273–9.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Dai X, Zhao PX. psRNATarget: a plant small RNA target analysis server. Nucleic Acids Res. 2011;39(suppl 2):W155–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Cristofani-Yaly M, Berger IJ, Targon MLPN, Takita MA, de Dorta SO, Freitas-Astua J, et al. Differential expression of genes identified from Poncirus Trifoliata tissue inoculated with CTV through EST analysis and in silico hybridization. Genet Mol Biol. 2007;30:972–9.

    Article  CAS  Google Scholar 

  54. Yang F, Wang G, Jiang B, Liu Y, Liu Y, Wu G, et al. Differentially expressed genes and temporal and spatial expression of genes during interactions between Mexican lime (Citrus aurantifolia) and a severe citrus tristeza virus isolate. Physiol Mol Plant Pathol. 2013;83(Supplement C):17–24.

    Article  CAS  Google Scholar 

  55. Cheng C, Zhang Y, Zhong Y, Yang J, Yan S. Gene expression changes in leaves of Citrus sinensis (L.) Osbeck infected by citrus tristeza virus. J Hortic Sci Biotechnol. 2016;91:466–75.

    Article  CAS  Google Scholar 

  56. Laino P, Russo MP, Guardo M, Reforgiato-Recupero G, Valè G, Cattivelli L, et al. Rootstock-scion interaction affecting citrus response to CTV infection: a proteomic view. Physiol Plant. 2016;156:444–67.

    Article  CAS  PubMed  Google Scholar 

  57. Fu S, Shao J, Zhou C, Hartung JS. Co-infection of sweet orange with severe and mild strains of citrus tristeza virus is overwhelmingly dominated by the severe strain on both the transcriptional and biological levels. Front Plant Sci. 2017;8:1419.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Tessitori M, Maria G, Capasso C, Catara G, Rizza S, De Luca V, et al. Differential display analysis of gene expression in Etrog citron leaves infected by citrus viroid III. Biochim Biophys Acta. 2007;1769:228–35.

    Article  CAS  PubMed  Google Scholar 

  59. Pérez-Clemente RM, Montoliu A, Vives V, López-Climent MF, Gómez-Cadenas A. Photosynthetic and antioxidant responses of Mexican lime (Citrus aurantifolia) plants to citrus tristeza virus infection. Plant Pathol. 2015;64:16–24.

    Article  Google Scholar 

  60. Du S-Y, Zhang X-F, Lu Z, Xin Q, Wu Z, Jiang T, et al. Roles of the different components of magnesium chelatase in abscisic acid signal transduction. Plant Mol Biol. 2012;80:519–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Shang Y, Yan L, Liu Z-Q, Cao Z, Mei C, Xin Q, et al. The mg-chelatase H subunit of Arabidopsis antagonizes a group of WRKY transcription repressors to relieve ABA-responsive genes of inhibition. Plant Cell. 2010;22:1909–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Iriti M, Faoro F. Abscisic acid is involved in chitosan-induced resistance to tobacco necrosis virus (TNV). Plant Physiol Biochem PPB. 2008;46:1106–11.

    Article  CAS  PubMed  Google Scholar 

  63. Alazem M, Lin K-Y, Lin N-S. The abscisic acid pathway has multifaceted effects on the accumulation of bamboo mosaic virus. Mol Plant-Microbe Interact MPMI. 2014;27:177–89.

    Article  CAS  PubMed  Google Scholar 

  64. Alazem M, He M-H, Moffett P, Lin N-S. Abscisic acid induces resistance against bamboo mosaic virus through Argonaute 2 and 3. Plant Physiol. 2017;174:339–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Whitty CD, Hall RH. A cytokinin oxidase in Zea mays. Can J Biochem. 1974;52:789–99.

    Article  CAS  PubMed  Google Scholar 

  66. Galuszka P, Frébort I, Sebela M, Sauer P, Jacobsen S, Pec P. Cytokinin oxidase or dehydrogenase? Mechanism of cytokinin degradation in cereals. Eur J Biochem. 2001;268:450–61.

    Article  CAS  PubMed  Google Scholar 

  67. Ashikari M, Sakakibara H, Lin S, Yamamoto T, Takashi T, Nishimura A, et al. Cytokinin oxidase regulates rice grain production. Science. 2005;309:741–5.

    Article  CAS  PubMed  Google Scholar 

  68. Novák J, Pavlů J, Novák O, Nožková-Hlaváčková V, Špundová M, Hlavinka J, et al. High cytokinin levels induce a hypersensitive-like response in tobacco. Ann Bot. 2013;112:41–55.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Masuta C, Tanaka H, Uehara K, Kuwata S, Koiwai A, Noma M. Broad resistance to plant viruses in transgenic plants conferred by antisense inhibition of a host gene essential in S-adenosylmethionine-dependent transmethylation reactions. Proc Natl Acad Sci U S A. 1995;92:6117–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Clarke SF, Burritt DJ, Jameson PE, Guy PL. Influence of plant hormones on virus replication and pathogenesis-related proteins in Phaseolus vulgaris L. infected with white clover mosaic potexvirus. Physiol Mol Plant Pathol U K. 1998;53:195–207.

  71. Pogány M, Koehl J, Heiser I, Elstner EF, Barna B. Juvenility of tobacco induced by cytokinin gene introduction decreases susceptibility to tobacco necrosis virus and confers tolerance to oxidative stress. Physiol Mol Plant Pathol. 2004;65:39–47.

    Article  Google Scholar 

  72. Baliji S, Lacatus G, Sunter G. The interaction between geminivirus pathogenicity proteins and adenosine kinase leads to increased expression of primary cytokinin responsive genes. Virology. 2010;402:238–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Vitti A, Nuzzaci M, Scopa A, Tataranni G, Remans T, Vangronsveld J, et al. Auxin and cytokinin metabolism and root morphological modifications in Arabidopsis thaliana seedlings infected with cucumber mosaic virus (CMV) or exposed to cadmium. Int J Mol Sci. 2013;14:6889–902.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Licciardello G, Scuderi G, Ferraro R, Giampetruzzi A, Russo M, Lombardo A, et al. Deep sequencing and analysis of small RNAs in sweet orange grafted on sour orange infected with two citrus tristeza virus isolates prevalent in Sicily. Arch Virol. 2015;160:2583–9.

    Article  CAS  PubMed  Google Scholar 

  75. Machida S, Yamahata N, Watanuki H, Owens RA, Sano T. Successive accumulation of two size classes of viroid-specific small RNA in potato spindle tuber viroid-infected tomato plants. J Gen Virol. 2007;88:3452–7.

    Article  CAS  PubMed  Google Scholar 

  76. Navarro B, Pantaleo V, Gisel A, Moxon S, Dalmay T, Bisztray G, et al. Deep sequencing of viroid-derived small RNAs from grapevine provides new insights on the role of RNA silencing in plant-viroid interaction. PLoS One. 2009;4:e7686.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Martinez G, Donaire L, Llave C, Pallas V, Gomez G. High-throughput sequencing of hop stunt viroid-derived small RNAs from cucumber leaves and phloem. Mol Plant Pathol. 2010;11:347–59.

    Article  CAS  PubMed  Google Scholar 

  78. Li R, Gao S, Hernandez AG, Wechter WP, Fei Z, Ling K-S. Deep sequencing of small RNAs in tomato for virus and viroid identification and strain differentiation. PLoS One. 2012;7:e37127.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Su X, Fu S, Qian Y, Xu Y, Zhou X. Identification of hop stunt viroid infecting Citrus limon in China using small RNAs deep sequencing approach. Virol J. 2015;12:1–5.

    Article  Google Scholar 

Download references


The authors would like to acknowledge Citrus Research International (CRI), the Technology and Human Resources for Industry Programme (THRIP) and the National Research Foundation (NRF) for their financial assistance towards this research. Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF.


This research project was funded by Citrus Research International (grant number 1100) as well as the Technology and Human Resources for Industry Programme (grant number TP13081327563).

Availability of data and materials

The datasets supporting the results of this article are available in the BioProject repository of the National Centre for Biotechnology Information, BioProject: PRJNA384115 in

Author information

Authors and Affiliations



MV participated in the design of the study, RNA samples preparation, performed data analysis and drafted the manuscript. GC participated in the design of the study, preparation of the healthy and inoculated plant material and contributed to drafting the manuscript. JTB participated in the design of the study and contributed to drafting the manuscript. HJM participated in the design of the study and contributed to drafting the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Hans J. Maree.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Table S1.

Results for the transcript differential expression analysis of ‘Marsh’ plants. Transcripts were considered to be differentially regulated as a result of infection if a |log2 fold change| > =1 and padj < =0.05 were observed. Table S2. Results for the transcript differential expression analysis of ‘Star Ruby’ plants. Transcripts were considered to be differentially regulated as a result of infection if a |log2 fold change| > =1 and padj < =0.05 were observed. Table S3. Results for the gene differential expression analysis across both ‘Marsh’ and ‘Star Ruby’ plants. Table S4. Biological process based characterisation of genes differentially expressed across both grapefruit varieties. Table S5. Molecular function based characterisation of genes differentially expressed across both grapefruit varieties. Table S6. Cellular component based characterisation of genes differentially expressed across both grapefruit varieties. (XLSX 392 kb)

Additional file 2: Figure S1.

Biological process based network of genes differentially expressed across both grapefruit varieties. Figure S2. Molecular function based network of genes differentially expressed across both grapefruit varieties. Figure S3. Cellular component based network of genes differentially expressed across both grapefruit varieties. (PDF 2159 kb)

Additional file 3: Table S7.

miRNA prediction results. Table S8. Homologous plant miRNA results. (XLSX 26 kb)

Additional file 4: Table S9.

Results for the sRNA differential expression analysis of ‘Marsh’ plants. sRNAs were considered to be differentially regulated as a result of infection if a |log2 fold change| > =1 and padj < =0.05 were observed. Table S10. Results for the sRNA differential expression analysis of ‘Star Ruby’ plants. sRNAs were considered to be differentially regulated as a result of infection if a |log2 fold change| > =1 and padj < =0.05 were observed. Table S11. Results for the sRNA differentially expressed across both ‘Marsh’ and ‘Star Ruby’ plants. sRNAs were considered to be differentially regulated as a result of infection if a |log2 fold change| > =1 and padj < =0.05 were observed. (XLSX 106 kb)

Additional file 5: Table S12.

Differential expression analysis results showing sRNAs with anti-correlated expression to their targets. (XLSX 53 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Visser, M., Cook, G., Burger, J.T. et al. In silico analysis of the grapefruit sRNAome, transcriptome and gene regulation in response to CTV-CDVd co-infection. Virol J 14, 200 (2017).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: