Degradome Analysis of Tomato and Nicotiana benthamiana Plants Infected with Potato Spindle Tuber Viroid

Viroids are infectious non-coding RNAs that infect plants. During infection, viroid RNAs are targeted by Dicer-like proteins, generating viroid-derived small RNAs (vd-sRNAs) that can guide the sequence specific cleavage of cognate host mRNAs via an RNA silencing mechanism. To assess the involvement of these pathways in pathogenesis associated with nuclear-replicating viroids, high-throughput sequencing of sRNAs and degradome analysis were carried out on tomato and Nicotiana benthamiana plants infected by potato spindle tuber viroid (PSTVd). Both hosts develop similar stunting and leaf curling symptoms when infected by PSTVd, thus allowing comparative analyses. About one hundred tomato mRNAs potentially targeted for degradation by vd-sRNAs were initially identified. However, data from biological replicates and comparisons between mock and infected samples reduced the number of bona fide targets—i.e., those identified with high confidence in two infected biological replicates but not in the mock controls—to only eight mRNAs that encode proteins involved in development, transcription or defense. Somewhat surprisingly, results of RT-qPCR assays revealed that the accumulation of only four of these mRNAs was inhibited in the PSTVd-infected tomato. When these analyses were extended to mock inoculated and PSTVd-infected N. benthamiana plants, a completely different set of potential mRNA targets was identified. The failure to identify homologous mRNA(s) targeted by PSTVd-sRNA suggests that different pathways could be involved in the elicitation of similar symptoms in these two species. Moreover, no significant modifications in the accumulation of miRNAs and in the cleavage of their targeted mRNAs were detected in the infected tomato plants with respect to the mock controls. Taken together, these data suggest that stunting and leaf curling symptoms induced by PSTVd are elicited by a complex plant response involving multiple mechanisms, with RNA silencing being only one of the possible components.


Introduction
Viroids are non-protein coding RNAs that infect plants [1]. According to their structural, functional and biological properties, viroids have been classified in the families Pospiviroidae (type member, potato spindle tuber viroid, PSTVd) and Avsunviroidae (type member, avocado sunblotch viroid, ASBVd) that group nuclear and chloroplast replicating viroids, respectively [2,3]. Members of both families may cause plant diseases, but how viroid RNAs may elicit pathogenesis is a question only partially answered [4].
Both nuclear and chloroplast replicating viroids are triggers and targets of RNA silencing, a regulatory mechanism of gene expression in most eukaryotes [5]. RNA silencing is mediated by microRNAs (miRNAs) and small interfering RNAs (siRNAs), which are small RNAs (sRNAs) of 21-24 nt originated by type III RNases (Dicer-like, DCLs in plants) PSTVd-infected (P) tomato and N. benthamiana (Nb) plants were separated by denaturing PAGE in 5% (c) or 17% (d) gels and analyzed by Northern-blot hybridization with a DIG-labeled riboprobe for detecting PSTVd (+) strand. Equal loading was assessed by the intensity of the bands corresponding to rRNAs after staining with ethidium bromide (c,d, bottom).
These assays revealed that the circular and linear forms of PSTVd and the (+) PSTVd-sRNAs accumulated at similar levels in both inoculated tomato replicates. The same RNA preparations were then used to generate six libraries of sRNAs and six libraries of degradome RNAs corresponding to the two mock and two infected tomato biological replicates and to the infected and mock N. benthamiana plants. Sequencing of each sRNA library yielded a total of 39-51 million high quality reads 18-26 nt long. Up to 14-17% of the total sRNAs reads from the infected samples were derived from PSTVd, with very similar levels of PSTVd-sRNAs (about 6.3-7.3 million) in the three libraries. In agreement with previous studies [20,26,27], PSTVd-sRNAs of both polarities, with a predominance of those of plus polarity, were identified ( Figure S1). When (+) and (−) PSTVd-sRNAs from the two tomato replicates and one N. benthamiana plant infected by PSTVd were mapped on PSTVd genomic and complementary RNAs, an irregular distribution of the reads was found, with very similar profiles observed for all three libraries ( Figure S2). These findings suggest that the observed read distributions may be host species-independent. These assays revealed that the circular and linear forms of PSTVd and the (+) PSTVd-sRNAs accumulated at similar levels in both inoculated tomato replicates. The same RNA preparations were then used to generate six libraries of sRNAs and six libraries of degradome RNAs corresponding to the two mock and two infected tomato biological replicates and to the infected and mock N. benthamiana plants. Sequencing of each sRNA library yielded a total of 39-51 million high quality reads 18-26 nt long. Up to 14-17% of the total sRNAs reads from the infected samples were derived from PSTVd, with very similar levels of PSTVd-sRNAs (about 6.3-7.3 million) in the three libraries. In agreement with previous studies [20,26,27], PSTVd-sRNAs of both polarities, with a predominance of those of plus polarity, were identified ( Figure S1). When (+) and (−) PSTVd-sRNAs from the two tomato replicates and one N. benthamiana plant infected by PSTVd were mapped on PSTVd genomic and complementary RNAs, an irregular distribution of the reads was found, with very similar profiles observed for all three libraries ( Figure S2). These findings suggest that the observed read distributions may be host species-independent.

Identification of Tomato mRNAs Cleaved by miRNA via RNA Silencing
Between 38-57 million clean tags, representing the 5 ends of uncapped, polyadenylated RNAs, were obtained in the degradome libraries prepared from the two mock-and two PSTVd-inoculated tomato replicates. While tags resulting from the cleavage of mRNAs 4 of 19 directed by host miRNAs and siRNAs were expected in all samples, the PSTVd-infected replicates should also contain tags corresponding to host mRNAs specifically targeted by PSTVd-sRNAs. In addition, because degradome analysis also identifies random degradation products that are not generated by endonucleolytic cleavage mediated by RNA silencing machinery [28], false positives were also expected. To exclude as many as possible of these false positive targets of PSTVd-sRNAs, we adopted several controls and stringent screening conditions. Putative cleaved targets identified by PAREsnip [28] were classified into five categories (i.e., categories 0-4) according to the abundance of degradome tags at the cleaved position. This classification indicates the confidence level of the sRNAmediated cleavages, those in category 0 being the most certain and those in categories 3 and 4 the most uncertain, likely corresponding to background noise generated by random degradation.
As a preliminary test of the quality of the degradome sequencing results, we used the PAREsnip software to identify in the four degradome libraries those mRNAs that were cleaved by the conserved tomato miRNAs. A total of 84 mRNAs were identified as targets of miRNAs belonging to 24 different miRNA families (Table S1), and most (75%) miRNA/mRNA pairs were classified as category 0 or 1 that, according to PAREsnip categorization, provide the strongest empirical evidence for a true cleavage at the specific position. This result was highly reproducible among the four libraries. Most of the miRNA targets identified in this study have been previously reported [23,[29][30][31][32][33]. However, 14 appeared to be new potential targets of previously reported miRNAs (Table S1).
Target plots (t-plots) were generated for several selected tomato mRNA targets of miRNAs, including miR-156, miR-167, miR-169, miR-172, mi-R319, miR-396, miR-398, and miR-408 to verify the quality of the sequenced libraries ( Figure S3). Similar low levels of background degradation were observed for each transcript in libraries from mock-and PSTVd-inoculated samples ( Figure S3), highlighting the reproducibility of the degradome analyses and the consistency of our results with those from other studies [23,[29][30][31][32][33]. When data from the mock and PSTVd-infected samples were compared, changes in the abundance of degradome tags of the targeted mRNAs were not significant (p < 0.05, Table S1) suggesting that miRNA-mediated cleavages remained largely unaffected upon PSTVd infection. Based on data from the corresponding sRNA libraries, levels of all but four miRNAs (i.e., miR169e-3p, miR169-5p, miR-171, and miR-4376) were also not significantly affected upon PSTVd infection ( Figure S4). Notably, changes in the abundance of these four miRNAs were not associated with a significant parallel change in the cleavage of the specific targets (Table S1). Altogether, these data, besides suggesting that PSTVd infection did not significantly affect the cleavage activity mediated by miRNAs in tomato, provided evidence of the high quality and repeatability of the degradome data from tomato plants.

Identification of Tomato mRNAs Cleaved by PSTVd-sRNAs via RNA Silencing
The mRNAs potentially targeted by PSTVd-sRNAs in the viroid-infected tomato plants were identified using the PAREsnip software and a subset of 21-nt PSTVd-sRNAs extracted from the respective sequenced sRNA libraries. Almost 100 potentially targeted mRNAs were identified in each library from infected samples (Tables S2 and S3). These mRNAs contained potential cleavage sites belonging to all categories (from 0 to 4), with category 0 sites being 17-22% of the total identified targets ( Figure 2a). Surprisingly, 60-80 possible mRNA targets were also identified in the two degradome libraries from the mock-inoculated tomato plants, when PAREsnip was applied using all 21-nt PSTVd sRNAs potentially generated by the variant PSTVd-Nb (Tables S4 and S5). In this case, the cleavage sites of the identified mRNAs also belonged to all PAREsnip categories, but targets with cleavage sites of category 0 were significantly less abundant (2-3%) than in the PSTVd-infected samples (Figure 2a). Since mock-inoculated plants did not contain PSTVd-sRNAs, all the potential targets identified in these samples must correspond to degradation products not generated by vd-sRNA-directed endonucleolytic cleavages, and therefore can be considered as non-PSTVd specific (i.e., false positive). When these false targets were removed from the list of potential targets identified in the two PSTVd-infected libraries, a total of 49 and 56 targets of PSTVd-sRNAs belonging to all four categories were left in the T3 and T4 libraries, respectively (Figure 2b), hereafter denoted as PSTVd specific targets.
cleavage sites of the identified mRNAs also belonged to all PAREsnip categories, but targets with cleavage sites of category 0 were significantly less abundant (2-3%) than in the PSTVd-infected samples (Figure 2a). Since mock-inoculated plants did not contain PSTVd-sRNAs, all the potential targets identified in these samples must correspond to degradation products not generated by vd-sRNA-directed endonucleolytic cleavages, and therefore can be considered as non-PSTVd specific (i.e., false positive). When these false targets were removed from the list of potential targets identified in the two PSTVd-infected libraries, a total of 49 and 56 targets of PSTVd-sRNAs belonging to all four categories were left in the T3 and T4 libraries, respectively (Figure 2b), hereafter denoted as PSTVd specific targets. Figure 2. (a) Total number and quality categorization (from 0 to 4, where 0 represents the strongest evidence for true cleavage products) of mRNAs targeted by PSTVd-sRNAs in PSTVd-infected (T3 and T4) and of non-PSTVd specific (false positive) mRNAs targets in mock-inoculated (T1 and T2) tomato biological replicates identified by degradome analysis; (b) number of PSTVd specific and non-specific mRNA targets in each viroid-infected tomato plant (T3 and T4) and degradome category. The PSTVd-specific mRNA targets were those not identified in any of the two mock-inoculated tomato replicates (T1 and T2), which are all considered as false positives.
Taking into consideration the fact that vd-sRNAs targeting host mRNAs are expected to be incorporated into AGO1 and that the nucleotide at the 5′ end of the vd-sRNA contributes to its loading into AGO1 [12], we examined the 5′ terminus of each of the respective PSTVd-sRNA. Interestingly, vd-sRNAs involved in cleavage of PSTVd specific targets falling into category 0 had exclusively U or A at their 5′ terminus, as expected for AGO1 loading (Figure 3a,b). In contrast, a similar bias was not observed at the 5′ nt of the vd-sRNAs from the other categories or for the vd-sRNAs targeting mRNAs classified as probable false positives in the degradome from mock-inoculated controls (see above, Figure  Taking into consideration the fact that vd-sRNAs targeting host mRNAs are expected to be incorporated into AGO1 and that the nucleotide at the 5 end of the vd-sRNA contributes to its loading into AGO1 [12], we examined the 5 terminus of each of the respective PSTVd-sRNA. Interestingly, vd-sRNAs involved in cleavage of PSTVd specific targets falling into category 0 had exclusively U or A at their 5 terminus, as expected for AGO1 loading (Figure 3a,b). In contrast, a similar bias was not observed at the 5 nt of the vd-sRNAs from the other categories or for the vd-sRNAs targeting mRNAs classified as probable false positives in the degradome from mock-inoculated controls (see above, Figure 3c,d).
These data indicate that the chance of identifying false targets of PSTVd-sRNAs in the degradome libraries from infected samples is lower when the cleavage sites fall within the category 0. In addition, since the two biological replicates showed identical symptoms, we reasoned that key mRNA(s) targeted by vd-sRNAs and involved in pathogenesis should be found in both replicates. Consequently, category 0 mRNAs that were present in both replicates after the exclusion of potential false positives were further considered as bona fide targets of PSTVd-sRNAs. These correspond to eight tomato mRNAs targeted by eight different PSTVd-sRNAs derived mainly from the plus polarity strand, with only two of them being of negative polarity ( Table 1). Irrespective of polarity, all PSTVd-sRNAs mediating the cleavage of the bona fide targets originated outside the hot-spots observed in the distribution profile of PSTVd-sRNAs on the viroid RNA ( Figure S2), indicating that sRNA abundance may not be directly correlated with mRNA targeting. Interestingly, these PSTVd-sRNAs are derived from sequence regions that are conserved among several PSTVd strains showing different pathogenicity; e.g., the intermediate, mild, RG1 and lethal variants (Table 1 and Figure S5). Only two of them, with 5 ends mapping at positions 192 and 290, are not conserved in the mild variant ( Table 1).
The bona fide targeted mRNAs encode eight proteins annotated as CWF19 (CWF19, Solyc09g083420.2.1), Little nuclei 1 (LINC1, Solyc03g045050.2.1), subunit of origin recognition complex protein 6 (ORC6, Solyc05g007450.2.1), NAC domain containing protein 38 (NAC038, Solyc04g079940.2.1), Supercentipede 1 (SCN1, Solyc11g012270.1.1), cyclin-dependent serine/threonine kinase F-1 (CDKF-1, Solyc12g007200.1.1), putative serinethreonine kinase (PBS1, Solyc05g024290.2.1) and proton pump interactor 1 (PPI-1, Solyc07g0 63100.2.1) ( Table 1). Four of these proteins (CWF19, NAC38, LINC1 and ORC6) are localized in the nucleus, two of them being transcription factors. CWF19 belongs to the CWFJ-like family protein of transcription factors containing a CCCH-type zinc finger and is likely involved in mRNA splicing, as denoted by its GO classification. NAC38 is a member of the NAC-family of plant transcription factors involved in developmental processes as well as in the responses to abiotic and biotic stresses, plant hormonal control, and defense [34,35]. In Arabidopsis, LINC1 encodes a coiled-coil nuclear protein that acts as a positive regulator of pathogen-associated molecular pattern-triggered immunity (PTI) responses, thereby regulating jasmonic acid signaling and gibberellin biosynthesis [36]. This protein also contributes to the shape and size of the nucleus [37]. ORC6 is one of the components required for assembly of the pre-replication complex necessary to initiate DNA replication in Arabidopsis [38,39], and it is involved in the BRP4-mediated mitotic cell-cycle progression [40]. The other four targeted mRNAs (SCN1, CDKF1, PBS-1, and PPI-1) correspond to genes encoding proteins that are not localized in the nucleus. SCN1 encodes a RhoGTPase GDP dissociation inhibitor (RhoGDI), which regulates root hair growth in Arabidopsis [41]. CDKF-1 encodes a cyclin-dependent serine/threonine protein kinase that is involved in activating phosphorylation of cyclin dependent kinases in Arabidopsis [42] and is a major regulator of cell proliferation and cell expansion [43]. Interestingly, Arabidopsis cdkf;1-1 mutants exhibited defects in cell division and cell elongation, leading to severe growth retardation in both shoots and roots [43]. In Arabidopsis, the serine-threonine kinase PBS1 is involved in plant defense against bacteria [44]. Arabidopsis thaliana PPI1 homolog protein is involved in proton transport through the plasma membrane [45,46]. These data indicate that the chance of identifying false targets of PSTVd-sRNAs in

Evidence for Cleavage of the Bona Fide Targets by PSTVd-sRNAs
Analyses of cleavages throughout the sequence of the eight bona fide target mRNAs (t-plot analyses) showed that, in the PSTVd-infected samples, the frequency of 5 ends consistent with the specific PSTVd-sRNA-mediated cleavage was always higher than the respective degradation background, with similar results obtained for both biological replicates ( Figure 4, Table 1). A comparable degradation background was observed in the t-plots corresponding to the same mRNAs from mock-inoculated samples where, as expected, mRNAs with 5 ends consistent with PSTVd-sRNAs-mediated cleavages were not found ( Figure 4). Therefore, in contrast to previous reports regarding other PSTVd-sRNA-targeted tomato mRNAs coding for CLC-b-like and RPS3a-like proteins [21], cleavage of the bona fide target mRNAs identified here is not associated with a vd-sRNAindependent degradation induced by viroid infection at late stages. To further validate these results, the PSTVd-sRNAs cleavage sites of three of the identified bona fide targets (Solyc09g083420.2.1, Solyc03g045050.2.1 and Solyc11g012270.1.1 encoding CWF19, LINC1 and SCN1, respectively) were determined by 5 RLM-RACE. All cDNA clones derived from Solyc09g083420.2.1 and Solyc11g012270.1.1 mRNA fragments had 5 termini compatible with a cleavage mediated by the respective targeting PSTVd-sRNAs ( Figure 4); In the case of Solyc03g045050.2.1, six out of nine clones had cDNAs with 5 termini corresponding to the predicted cleavage site; While, in three clones, the termini mapped one nucleotide downstream in the mRNA (Figure 4). Taken together, these data confirmed the cleavage of targeted mRNAs according to an RNA silencing mediated mechanism guided by the respective PSTVd-sRNAs.

Correlation between Targeting by PSTVd-sRNAs and Steady State Level of Cleaved Host mRNAs
The effect of the mRNA silencing-mediated cleavage on the steady state levels of the eight bona fide targets of PSTVd-sRNAs was tested by RT-qPCR assays on three biological replicates from three independent inoculation experiments. To amplify cDNAs derived exclusively from the uncleaved mRNAs, specific PCR primers were designed to bind upstream and downstream from the cleavage site of each target. Only the mRNAs encoding CWF19, LINC1, ORC6, and SCN1 showed a significantly lower level of accumulation in the PSTVd-infected than in the mock-inoculated tomato plants (Figure 5), supporting a role for the respective PSTVd-sRNAs in the regulation of the expression of the cognate genes. In contrast, no significant difference was observed in the accumulation of the other four bona fide targets, indicating that cleavage induced by the specific PSTVd-sRNAs is not always followed by significant changes in the steady state level of the cleaved mRNAs ( Figure 5).

Identification of N. benthamiana mRNAs Potentially Cleaved by miRNAs and PSTVd-sRNAs via RNA Silencing
Sequencing of degradome libraries of mock and PSTVd-infected N. benthamiana plants generated 23 and 28 million of clean reads, respectively. PAREsnip analysis identified a total of 79 mRNAs cleaved by conserved miRNAs belonging to 15 different families. Approximately 70% of the cleavage sites in both the mock-inoculated and PSTVd-infected samples were classified as category 0 or 1 according to the PAREsnip categorization (Table S6). This high degree of specificity was consistent with data reported in the literature [47], thus validating the results of this degradome analysis.   Paralleling the search strategy adopted for tomato, a total of 184 or 149 potential PSTVd-sRNA targets belonging to all PAREsnip categories were identified in the RNA samples isolated from PSTVd-infected and mock-inoculated N. benthamiana, respectively (Tables S7 and S8). Nearly half of the targets identified in the PSTVd-infected plant belonged to category 0 ( Figure 6a). Subtraction of false positive targets identified in the mockinoculated plant from the potential targets identified in the PSTVd-infected library yielded a total of 49 targets of category 0 (Figure 6b and Table S9). As observed for tomato, most (a total of 37) mRNAs containing category 0 cleavage sites were targeted by PSTVd-sRNAs with a U at their 5 terminus (Figure 6c), supporting their preferential loading into AGO1. A similar bias was not observed for either vd-sRNAs in the categories 1-4 or those in category 0 identified in the mock-inoculated controls (false positive targets, Figure 6d). t-plots of representative targets of PSTVd-sRNAs in the mock and PSTVd-infected samples showed comparable degradation backgrounds ( Figure S6). Therefore, as observed in tomato, PSTVd-infection in N. benthamiana was not accompanied by a non-specific degradation of mRNAs targeted by PSTVd-sRNA. In the absence of a replicate, it was not possible to determine which were bona fide targets of PSTVd-sRNAs in N. benthamiana; thus, some of the potential targets could be false. It should be noted, however, that none of the 49 mRNAs potentially targeted by PSTVd-sRNAs in N. benthamiana (Table S5) corresponded to bona fide or potential targets identified in tomato infected by the same strain of PSTVd and showing similar symptoms. Moreover, although the PSTVd-sRNAs mediating the cleavage of the bona fide targets in tomato were also identified in N. benthamiana, they are not expected to target the N. benthamiana homologs because nucleotide changes in the mRNA sequence would destabilize the corresponding vd-sRNA/mRNA duplexes and impair their AGO-mediated cleavage (Table S10). sponded to bona fide or potential targets identified in tomato infected by the same strain of PSTVd and showing similar symptoms. Moreover, although the PSTVd-sRNAs mediating the cleavage of the bona fide targets in tomato were also identified in N. benthamiana, they are not expected to target the N. benthamiana homologs because nucleotide changes in the mRNA sequence would destabilize the corresponding vd-sRNA/mRNA duplexes and impair their AGO-mediated cleavage (Table S10). (d) When the same analysis was performed considering the PSTVd-sRNAs targeting the false (non-PSTVd-specific) targets identified in the degradome sequencing from mock-inoculated samples, no Figure 6. (a) Total number and quality categorization (from 0 to 4, where 0 represents the strongest evidence for true cleavage products) of mRNA targets of PSTVd-sRNA and non-PSTVd specific (false positive) mRNAs targets identified by degradome analysis in PSTVd-infected and mock-inoculated N. benthamiana, respectively; (b) number of specific and non-PSTVd specific mRNA targets in PSTVdinfected N. benthamiana plant in each degradome category. The PSTVd-specific mRNA targets were those not identified in any of the two mock-inoculated tomato replicates; (c) analysis of the 5 terminal nucleotide of PSTVd-sRNAs involved in the cleavage of PSTVd-specific targets according to the degradome category in PSTVd-infected N. benthamiana, the 5 terminal nucleotide of PSTVd-sRNAs involved in host mRNA cleavages classified into the category 0 were mainly U. (d) When the same analysis was performed considering the PSTVd-sRNAs targeting the false (non-PSTVd-specific) targets identified in the degradome sequencing from mock-inoculated samples, no bias was observed in the 5 terminal nucleotide of PSTVd-sRNAs targeting mRNAs of category 0, suggesting that they were not targeted by AGO1.

Discussion
Sequence-specific cleavage of host mRNAs mediated by vd-sRNAs via RNA silencing has been proposed to be one of the primary molecular lesions that could initiate a cascade of molecular events leading to the development of macroscopic symptoms in viroid infected plants. To assess the potential involvement of this mechanism in the pathogenesis of PSTVd, the type member of the family of nuclear-replicating viroids, we performed comparative analyses of sRNAs and degradome libraries from tomato and N. benthamiana plants either mock-inoculated or infected by the same PSTVd-Nb variant. These two solanaceous species were selected because they develop similar stunting and leaf curling symptoms when infected by PSTVd-Nb, thus providing the possibility of determining whether or not homologous genes or genes involved in the same developmental pathway were impaired by RNA silencing activated by PSTVd-sRNAs. A genome-wide analysis strategy based on HTS of sRNAs and degradome libraries was adopted to gain a global view of the host mRNAs targeted by vd-sRNAs in the two hosts.
When the degradome libraries were screened for mRNAs targeted by conserved miRNAs, previously reported cleavages of tomato mRNAs by conserved miRNAs were detected [23,[29][30][31][32][33]. Moreover, several new targets were identified, thereby confirming the depth of the degradome sequencing. Comparisons of sRNA data from the mock-inoculated and PSTVd-infected tomato replicates showed that most miRNAs accumulated to similar levels, with significant alterations observed for only four miRNAs. These changes were not associated with significant parallel quantitative changes in the cleavage of the respective targeted mRNAs, however. These data were consistent with those previously reported by Zheng and colleagues [23], who found that increased levels of two of these miRNAs (miR-171 and miR-4376) following PSTVd-infection were not accompanied by significant changes in the cleavage and accumulation of the respective targeted transcripts. Nevertheless, we did not observe significant alterations in the miR-167-and miR-393-guided cleavage of Solyc02g037530 (encoding auxin response factor 8) and Solyc09g074520 (encoding TIR1), respectively, as reported by the same authors.
Based on HTS data, alterations in the accumulation levels of some miRNAs have previously been reported in tomato plants infected by several PSTVd variants [48,49]. A third study by Owens and colleagues reported that miRNA levels were either unchanged or only moderately affected in PSTVd-infected plants [50]. Interestingly, no significant differences in the accumulation levels of any of the miRNAs examined or their respective mRNA targets were detected in tomato plants infected by citrus exocortis viroid, a viroid closely related to PSTVd [51]. This apparent contradiction calls for additional studies to conclusively assess the involvement of miRNAs in plant-viroid interaction.
To screen out potential false targets in the degradome libraries resulting from random degradation products, our search for tomato mRNAs targeted by PSTVd-sRNAs was carried out adopting stringent criteria that took into consideration both false positives identified in the degradome library from mock-inoculated plants and the statistical robustness of the mRNA cleavage. Moreover, to increase the chances of pinpointing real targets, only those targets identified in both degradome libraries from symptomatic plants were considered. Applying this approach, the almost 100 potential targets predicted based on degradome analysis in the absence of constraints were reduced to eight bona fide mRNA targets of PSTVd-sRNAs. All the PSTVd-sRNAs mediating the cleavage of the bona fide targets have either a U or an A nucleotide at the 5 end, thus satisfying the preferences for loading into AGO1 [12][13][14], the key enzyme operating in RNA silencing against endogenous and foreign RNAs [52][53][54]. Interestingly, these PSTVd-sRNAs accumulate to only relatively low levels (Table 1) and did not correspond to those forming the hot-spots observed in the distribution profiles of both PSTVd plus and minus RNAs ( Figure S2). A previous report noted the relatively low abundance of PSTVd-sRNAs involved in the silencing of endogenous tomato mRNAs [21], paralleling the situation observed for the chloroplast replicating viroid PLMVd [19]. Taken together, these data indicate that the biological significance of vd-sRNAs in silencing host mRNAs cannot be inferred solely on the basis of their relative abundance in the infected tissues.
The eight bona fide mRNAs targeted by PSTVd-sRNAs encode proteins having potential roles in relevant developmental, signaling, or defense pathways. These proteins are localized in either the nucleus or cytoplasm. In our study, the t-plots and RACE analysis of infected samples confirmed that the bona fide targets were predominantly cleaved at the predicted sites. Nevertheless, RT-qPCR assays showed that only four of the eight target mRNAs accumulated at significantly lower levels in the viroid-infected hosts than in the mock-inoculated plants. Therefore, the effective cleavage mediated by PSTVd-sRNAs did not produce significant changes in the steady-state level of at least some targeted mRNAs.
This observation is in agreement with results of a previous degradome analyses of PSTVd-infected tomato plants [23], although the targeted mRNAs in the latter study were different from those reported here. It is possible that a temporary transcriptional induction of the targeted mRNAs may compensate for the reduction due to cleavage. In our study, the targeted mRNAs whose accumulation remained apparently unchanged in the infected tissues encode either signaling proteins such as protein kinases CDKF-1 and PBS1, or proteins such as CDKF1 and NAC38 that are involved in specific developmental or defense pathways. Rather than being constitutively expressed, expression of such proteins would be expected to be induced during the viroid infection. Therefore, we cannot exclude the possibility that the basal PSTVd-sRNA-mediated cleavage of the cognate mRNAs may became physiologically relevant during specific plant developmental stages, especially if the transcription of the targeted mRNAs is simultaneously affected by other regulatory networks in infected tissues.
Among the mRNAs targeted by PSTVd-sRNAs and exhibiting a significantly lower steady state level in the infected tissues were those encoding nuclear proteins CWF19, which is likely involved in mRNA splicing, and the homolog of the Arabidopsis SCN1 protein, involved in root hair growth. Previous reports have documented both alterations in the alternative splicing of host protein-coding genes [23] and downregulation of the SCN1 gene [55] in PSTVd-infected tomato plants.
The fact that all the PSTVd-sRNAs targeting the bona fide mRNAs identified in this study derived from regions of the genome that are conserved among different PSTVd strains suggests that cleavage of some of the targeted mRNAs could be relevant for viroid replication and/or infectivity. Further studies are needed to explore this possibility. However, because different PSTVd strains induce symptoms of differing severity, a direct link between these eight particular PSTVd-sRNAs and viroid pathogenicity seems unlikely.
None of the tomato mRNAs identified as targets for degradation by PSTVd-sRNAs in previous studies [20][21][22][23]56] were detected in our degradome libraries, even when searches were extended to targets belonging to all categories or those not present in both replicates. Differences in PSTVd strains, tomato varieties, and the experimental conditions used in each study could explain these discrepancies. Taken together, these data indicate that (i) PSTVd-sRNAs may simultaneously target different mRNAs for sequence-specific cleavage in infected tomato tissues, (ii) PSTVd-sRNA-mediated cleavage is not necessarily associated with decreased accumulation of the targeted mRNA, and (iii) the mRNAs targeted likely differ depending on the PSTVd strain, tomato cultivar, and/or experimental conditions. In such a complex situation, it was not possible to identify which one, if any, of the targeted mRNAs is actually involved in the elicitation of the typical stunting and leaf curling symptoms observed in PSTVd infected tomato.
In the absence of biological replicates, selection of true targets in the N. benthamiana degradome was necessarily less precise, and, therefore, a higher number of possible PSTVd-sRNA targets was identified. Notably, none of these targets was a gene homologous to any of the eight bona fide tomato targets discussed, mainly because the mRNA sequences targeted by the PSTVd-sRNAs in tomato were not conserved in the homologous N. benthamiana genes (Table S6). The fact that both hosts exhibited similar stunting and leaf curling symptoms further supports the conclusion that the eight bona fide targets identified in tomato are likely not directly involved in the elicitation of these symptoms. However, we cannot exclude that they could play a role in plant-viroid interaction.
Sampling in our study was limited to the late stages of infection; thus, we cannot exclude the possibility that earlier targeting of one or more mRNAs by PSTVd-sRNAs could have triggered a signaling cascade leading to the later appearance of macroscopic symptoms. No evidence for such hypothetical targeted mRNA potentially eliciting pathogenesis was detected by our analyses, however. This is quite different from the situation with the RNA silencing-activated pathogenesis of the chloroplast replicating viroid PLMVd. In this case, vd-sRNA-mediated degradation of a specific host mRNA is always spatially and temporally associated with symptom expression, even at late stages [18,19].
Transcriptomic analyses have thoroughly documented the alteration by PSTVd (and other nuclear-replicating viroids) of the expression of genes involved in defense, hormone signaling, photosynthesis, RNA processing and binding, protein metabolism and modification, RNA silencing and plant innate immune responses [23,55,[57][58][59][60]. The ability of viroids to interfere with the translational machinery [61,62], as well as with RNA-dependent DNA methylation [63][64][65][66][67] has also been reported. Besides the sequence-specific degradation of some host mRNA(s) guided by PSTVd-sRNAs, all these molecular pathways should be taken into consideration when PSTVd pathogenesis is addressed.
In summary, applying degradome analyses, we found several host mRNAs targeted by PSTVd-sRNAs for a specific cleavage mediated by RNA silencing machinery. However, based on comparative analysis in tomato and N. benthamiana plants, the direct involvement of these mRNAs in PSTVd-induced stunting and leaf curling symptoms in these hosts appears unlikely. Additional degradome-based studies performed at several infection stages could further clarify the role of RNA silencing in triggering PSTVd pathogenesis.

Plant Material
Seedlings of tomato (Solanum lycopersicum L.) cv. 'Rutgers' at the cotyledon/first true leaf stage were mechanically inoculated with nucleic acid extracts from tomato plants infected with PSTVd (Nb variant, GenBank accession number AJ634596.1) in 0.2 M boric acid-NaOH, pH 9 or with the inoculation buffer (mock-inoculated plants). Nicotiana benthamiana seedlings were agroinoculated with Agrobacterium tumefaciens C58 carrying either an empty plant binary expression plasmid (pBIN19Ø) or one containing a head-totail dimeric insert of PSTVd-Nb under the control of the 35S promoter of cauliflower mosaic virus [68]. Plants were grown in a controlled environment chamber at 28 • C with fluorescent light for 16 h and at 25 • C in darkness for 8 h. Leaf samples were collected 30 days postinoculation (dpi). Three independent inoculation experiments were performed.

Total RNA Isolation and Northern Blot Hybridization
Total RNA was extracted from apical leaves of mock-inoculated and PSTVd-infected plants with TRizol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. Total RNA preparations were analyzed by polyacrylamide gel electrophoresis under denaturing conditions (8M Urea and 1× TBE buffer) in 5% gels, for detecting genomic PSTVd RNAs, and in 17% gels, for detecting vd-sRNAs. After ethidium bromide staining, the RNAs were electrotransferred to Hybond-N+ membranes (Roche Diagnostics GmbH, Basel, Switzerland) and were hybridized with digoxigenin-labeled fulllength riboprobes specific for detecting PSTVd of plus polarity as described previously [12]. Hybridization signals were revealed with anti-DIG alkaline phosphatase conjugate and the chemiluminescent substrate CSPD (Roche Diagnostics GmbH, Basel, Switzerland) and exposure to X-ray film. Ribosomal RNAs stained with ethidium bromide were used as equal loading controls.

Library Construction and High-Throughput Sequencing
cDNA libraries of sRNAs were generated as reported previously [26] and subjected to high-throughput sequencing on Hi-Seq 2500 (Illumina, San Diego, CA, USA) by Fasteris custom services (Fasteris, Geneva, Switzerland), generating 1 × 50 bp single-end reads.
Degradome libraries (Parallel Analysis of RNA Ends, PARE) were constructed according with the protocol described by German et al. [69]. Briefly, poly(A)+ RNA was isolated from 50 ug of total RNA obtained as indicated before. After ligation of an RNA adapter (Solexa Gex 1 adapter) to the 5 phosphate poly(A) + RNAs, the first cDNA strand was synthesized using an oligo-dT adapter primer, and the resulting cDNA was amplified by PCR. The amplicons were digested with MmeI and purified by PAGE before the Illumina 3 TruSeq adapter was ligated to the MmeI overhang. After PCR amplification using a modified Gex1 adapter, the library products were sequenced (pair-end 2 × 25 bp) on Illumina NextSeq 500 system (Illumina, San Diego, CA, USA) by Vertis Biotechnologie AG (Freising, Germany). Raw sequences of the sRNAs and degradome libraries have been deposited in the European nucleotide archive (ENA) under the project accession number PRJEB42233.

Analysis of sRNA and Degradome Libraries
Raw sRNAs reads generated by HTS were filtered for quality and, after removing the adapters, they were sorted by sequence length. Reads of size between 18 to 26 nt were mapped by Bowtie 1.2.2 [70] on the PSTVd-Nb reference sequence and those of 21-nt length were selected for further degradome analysis. Tomato and N. benthamiana conserved miRNAs were identified by mapping (Bowtie 1.2.2) against the plant miRNA in miRbase (Release 22.1) and the published tomato and N. benthamiana miRNAs.
After quality control, adapter trimming, and filtering for 100% overlapping paired reads, high quality degradome reads were analyzed, using PAREsnip software tool [28] with a cut-off p-value of 0.05, for the identification of PSTVd-sRNA and miRNA targets. The transcripts data set used for mapping degradome reads was the release ITAGv2.3 of Solanum lycopersicum, and the v1.0.1 N. benthamiana GENOME predicted cDNA (www.solgenomics.net, accessed on 30 March 2021). The sRNA data set used in the PAREsnip tool to identify potential cleaved host targets of PSTVd-sRNAs in the degradome data from infected samples was constituted by the 21-nt sRNAs mapped in PSTVd obtained in the corresponding sRNA library (see above). Degradome data from mock controls was analyzed using a sRNA data set generated in silico to contain all possible sRNA derived from PSTVd-Nb. The resulting targets of PSTVd-sRNA in mock samples were considered as false targets.
The putative cleavage sites identified by PAREsnip analysis in each degradome library were assigned to one of five categories (categories 0-4) according to degradome fragment abundance at the cleaved position as described in [28]. This classification indicates the confidence level of target prediction, those in category 0 being the most likely to be true targets. In category 0, maximal abundance of the degradome tag occurs at the predicted cleavage site in the transcript. Category 1 is defined similarly, but there is more than one cleavage site that is associated with degradome tags of similar abundance. Category 2 sites are cleavage positions that are associated with tags whose abundance is lower than the maximum but higher than the median for that transcript. Category 3 sites are defined as for category 2 but where the abundance of tags at the predicted cleavage position is lower than or equal to the median. In category 4, degradome tag abundance may be as low as one raw read at that position. Abundance values for degradome tags in categories 0-3 are always more than one raw read.

RNA Ligase Mediate-Rapid Amplification of CDNA Ends (5 RLM-RACE)
Aliquots (2 µg) of total RNA treated with turbo-DNase (Invitrogen, Carlsbad, CA, USA) were mixed with 0.25 µg of the RNA adapter (5 -CGACUGGAGCACGAGGACACUG ACAUGGACUGAAGGAGUAGAAA-3 ) and incubated at 95 • C for 2 min and snap-cooled on ice. Then, the mixture was incubated at 37 • C for one hour with 5U of T4 RNA ligase (Fermentas, Waltham, MA, USA), 1× buffer T4 RNA ligase (Fermentas, Waltham, MA, USA), 1 mM ATP and 10 U RNase inhibitor. First, cDNA synthesis of adapter-ligated RNAs was carried out with Superscript II (Invitrogen, Carlsbad, CA, USA) following the supplier instructions and using a specific primer complementary to the transcript of interest (Table S7). Then, the cDNA was amplified by PCR using the same primer of the first cDNA synthesis and a primer homologous to the RNA adapter (5 -TGGAGCACGAGGACACTGACATG-3 ) and Go-Taq DNA polymerase (Promega, Corporation, Madison, WI, USA). PCR products of the expected size (Table S11) were eluted from agarose gel and re-amplified in a nested PCR reaction using a primer complementary to the transcript of interest and internal with respect to the primer used in the first PCR and a primer homologous to the RNA adapter and nested within the previous one (5 -GACACTGACATGGACTGAAGGAGTAG-3 ). The PCR products were eluted from agarose gel, cloned in pGEMT-easy (Promega, Madison, WI, USA) and sequenced by Sanger custom sequencing service (Macrogen, Amsterdam, The Netherlands).

Quantitative RT-PCR (RT-qPCR) for Gene Expression Analysis
Total RNA extracted as described above was treated with Turbo DNase free (Invitrogene, Carlsbad, CA, USA) and 1 µg of the treated RNA was used for the first strand cDNA synthesis with Superscript II (Invitrogen, Carlsbad, CA, USA)) and random hexamers according to the manufacturer's protocol. The cDNA reaction was diluted 5-fold with water and used as template for qPCR amplification. The PCR reaction was carried out in a total volume of 20 µL, containing: 4 µL of 5× PyroTaqEvaGreen qPCR Mix Plus (Bioline, London, UK), 0.5 µL of 10 µM of each primer, and 2 µL of the diluted cDNA reaction. Three biological replicates from three independent inoculation experiments were tested, and three technical replicates were performed for each biological replicate. Negative control reactions without cDNA template for each primer pair were included in each amplification experiment. The specificity of amplification was confirmed by melting curve analysis. Gene-specific primer pairs used for qPCR amplification are listed in Table S11. Efficiencies for each primer pair were calculated from standard curves obtained using serial dilutions of mock-inoculated tomato cDNA. To verify the sequence of the amplified cDNA, the amplicon generated in one of the RT-qPCR reactions for each pair of primers was cloned and sequenced. The relative expression level (fold change) of the selected transcripts was calculated by the 2 ∆∆CT method using the ubiquitin conjugating enzyme (UBC) as an internal reference gene for normalization of expression levels. Statistical significance was assessed by a t-test analysis, and the results were considered significant at p < 0.05. Error bars represent the standard error of the mean.