Potato Virus Y Infection Alters Small RNA Metabolism and Immune Response in Tomato

Potato virus Y (PVY) isolate PVYC-to induces growth reduction and foliar symptoms in tomato, but new vegetation displays symptom recovery at a later stage. In order to investigate the role of micro(mi)RNA and secondary small(s)RNA-regulated mechanisms in tomato defenses against PVY, we performed sRNA sequencing from healthy and PVYC-to infected tomato plants at 21 and 30 days post-inoculation (dpi). A total of 792 miRNA sequences were obtained, among which were 123 canonical miRNA sequences, many isomiR variants, and 30 novel miRNAs. MiRNAs were mostly overexpressed in infected vs. healthy plants, whereas only a few miRNAs were underexpressed. Increased accumulation of isomiRs was correlated with viral infection. Among miRNA targets, enriched functional categories included resistance (R) gene families, transcription and hormone factors, and RNA silencing genes. Several 22-nt miRNAs were shown to target R genes and trigger the production of 21-nt phased sRNAs (phasiRNAs). Next, 500 phasiRNA-generating loci were identified, and were shown to be mostly active in PVY-infected tissues and at 21 dpi. These data demonstrate that sRNA-regulated host responses, encompassing miRNA alteration, diversification within miRNA families, and phasiRNA accumulation, regulate R and disease-responsive genes. The dynamic regulation of miRNAs and secondary sRNAs over time suggests a functional role of sRNA-mediated defenses in the recovery phenotype.


Introduction
In the last two decades, thanks to the recognition of noncoding small RNAs (sRNAs) as leading regulatory factors, studies on long-established mechanisms that intervene in the evolution, development, and control of many biological processes have seen an impressive degree of progress. Different classes of sRNAs have been identified to date as important regulators of gene expression and adaptive responses to adverse conditions [1,2]. sRNAs regulate different plant biological processes, such as physiology, alters the expression of sRNAs in its solanaceous hosts depending on the plant sensitivity. Certain sRNAs have been shown to be involved in the establishment of tolerance and resistance responses to PVY infection in potato and tobacco plants [13,28]. Additionally, several studies have reported that, similar to other plant hosts, tomato responds to viral pathogens by generating infection-responsive sRNAs of different classes including miRNAs and phasiRNAs [29][30][31][32][33]. In this study, we examined PVY-tomato interactions, looking into the complexity and diversity of defense pathways involving sRNAs. We used a combination of experimental and computational approaches in order to identify and characterize PVY-responsive sRNAs in infected tomato plants. We generated four sRNA libraries from PVY C -to-infected and mock-inoculated tomato leaves at two time points in order to capture the process of recovery in tomato. Using high-throughput sequencing (sRNA-seq), we performed a genome-wide identification and expression profiling of miRNAs and secondary sRNAs. The data herein offer new insights into sRNA-regulated host responses to PVY infection, and indicate a functional role of sRNAs in the molecular mechanisms through which disease tolerance in tomato is controlled.

Plant and Virus Material
A PVY isolate belonging to the strain PVY C , found on cultivated tomato in Italy and named PVY C -to (accession no. EU482153), was maintained on and purified from tobacco (Nicotiana tabacum cv. Samsun) plants as described [34]. Two-week-old seedlings of tomato cv. UC82 were used for plant inoculation. For PVY inoculation, a 10 µL aliquot of a purified virus preparation at 50 ng/µL in 30 mM Na 2 HPO 4 , or the same buffer alone as the negative mock-inoculated control, was rubbed on the cotyledons of twelve plants per inoculum.

RNA Extraction
Symptomatic and recovered (symptom-free) leaves from six PVY C -to-inoculated plants were collected at 21 and 30 dpi, respectively. Leaves from six mock-inoculated plants were sampled at the same time points as the negative controls. The leaves from each treatment were pooled together into four samples.
Total RNA was extracted from about 100 mg of leaf tissue using the TRIzol reagent (Thermo Fisher Scientific, Waltham, MA USA) according to the manufacturer's instructions. RNA preparations were subjected to DNase digestion using TURBO DNA-free™ Kit (Thermo Fisher Scientific) and resuspended in nuclease-free water. Total RNA was used for the detection of viral RNA and sRNA sequencing. The leaves from an additional set of six mock-and six PVY C -to-inoculated tomato plants at 21 dpi were collected and pooled into three biological replicates per treatment for separate extraction of high and low molecular weight RNA fractions (HMW and LMW, respectively) using the mirVana miRNA isolation kit™ (Ambion, Austin TX, USA), according to the manufacturer's instructions. These RNA preparations, after DNase treatment, were used for gene and miRNA quantification by reverse transcription-quantitative polymerase chain reaction (RT-qPCR). RNA integrity was estimated by agarose gel electrophoresis on a 1.2% agarose gel. RNA concentration and purity were measured with a Nanodrop ND-1000 Spectrophotometer (Nanodrop Technologies, Thermo Fisher Scientific).

Primer Design and RT-qPCR Analysis for Detection of Viral RNA, mRNAs and miRNAs
Six RS-related and four disease-related genes were selected among those described as miRNA targets in tomato. Primer pairs used for the tomato gene quantification or PVY genomic RNA detection were designed with the support of the IDT's PrimerQuest software (http://eu.idtdna.com/Scitools/ Applications/Primerquest/) (Table S11). For viral RNA and tomato mRNA quantitative analysis, RT-qPCR was performed as previously described [35]. Briefly, RNA samples were reverse-transcribed using the High Capacity cDNA Reverse Transcription Reagents (Thermo Fisher Scientific), according to the manufacturer's instructions. qPCR reactions were performed in a 96-well CFX96 RealTime PCR System (Bio-Rad, Hercules, CA, USA) using 2X PowerUp™ SYBR Green PCR Master Mix (Thermo Fisher Scientific), 400 nM forward and reverse primers, and 1 µL of a 1:10 dilution (approx. 5 ng) of reverse-transcribed RNA in a total volume of 10 µL.
The miRNA 1st-Strand cDNA Synthesis Kit and miRNA QPCR Master Mix (Agilent Technologies) were used for RT-qPCR analysis of 21 miRNAs in order to validate sRNA-seq experiments. cDNA was amplified using the provided universal reverse primer and the miRNA-specific forward primers (Table S11) according to the manufacturer's instructions. Forward primers were homologous to the selected known and novel miRNAs, and included an additional 6-nt sequence at the 5' terminus for the amplicon stabilization [22].
Relative viral RNA, mRNA, and miRNA expression levels were calculated according to the (2 -∆∆Ct ) method using tomato UBI3 (for viral RNA, mRNAs) or both tomato snoU6 and miR167a (for miRNAs) as the endogenous controls and corrected for PCR amplification efficiencies as previously described [22,35]. Data were statistically analyzed by Student's t-test.
For the validation of sRNA-seq data with relative miRNA expression levels estimated by RT-qPCR, sRNA-seq relative expression values were calculated by summing the normalized RPM counts of mature miRNAs that perfectly matched with the specific forward primers used for RT-qPCR, and RPM counts of all isomiRs that matched at least three nucleotides at the 3'-end of the primers.

Sequencing and sRNA Data Analysis
Four high quality sRNA libraries were constructed using the TruSeq Small RNA Library Prep Kit (Illumina) and sequenced in a single-read mode (50 bp) on a HiSeq 2000 (Illumina) by IGA Technology Services (Udine, Italy).
To identify novel unannotated miRNAs and their loci of origin (MIR loci), reads were submitted to the two plant miRNA prediction tools ShortStack 3.8.3 [40] and miR-PREFeR [41]. Predictions were performed using default parameters, except that no mismatches were allowed during mapping on the reference tomato genome (Genome build SL3.0, solgenomics.net). As novel miRNAs were considered only with more than ten raw reads in at least two sRNA libraries, the miRNA sequence and corresponding miRNA* and MIR locus should have been predicted with both miRNA prediction tools. Within the prediction analyses, the reads that could be mapped to more than 30 locations in the tomato genome were also discarded as being too repetitive to be miRNAs. The output of miRNA prediction tools also contained the predictions of already annotated tomato MIR loci; therefore, to separate them from potential novel MIR loci candidates, annotated pre-miRNA precursors from the miRBase database release 22 were mapped to the reference tomato genome (Genome build SL3.0, solgenomics.net) using bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml) [42]. Next, genome locations were extracted and compared with predicted MIR loci locations using our internally developed script. If no overlap was detected, the predicted MIR loci were regarded as novel. Novel tomato miRNAs were further classified into known or novel miRNA families by clustering their predicted pre-miRNA sequences with sequences of known plant pre-miRNAs from miRBase using CD-HIT-EST with an identity threshold of 0.8 (http://weizhong-lab.ucsd.edu/cdhit-web-server/cgi-bin/index.cgi?cmd=cd-hit-est) [43]. The sequences showing similarities with annotated pre-miRNAs were grouped into corresponding known miRNA families, and sequences that did not show similarity with known plant miRNAs were classified as novel.
Additionally, sequence miRNA variants (isomiRs) of known and novel miRNAs were identified using computational pipeline isomiRID [44]. Only sRNAs perfectly matching known or novel pre-miRNA sequences, known as templated isomiRs [45], were considered. In addition, to ensure robustness of isomiR predictions, only sequences present at a minimum of 5 raw reads in at least two sRNA libraries were regarded as isomiRs.
Predictions of PHAS loci was performed using ShortStack 3.8.3. PHAS loci were detected by mapping preprocessed sRNA reads to tomato transcriptome sequences (ITAG release 3.2, solgenomics.net). Analyses of phasing were performed in 21-and 24-nt intervals using the default settings.

Identification of Differentially Expressed miRNAs
Preprocessed reads from sRNA-seq samples were mapped with no mismatches to all identified known, novel miRNAs and isomiRs, and counted. sRNA counts were further normalized to total library sizes. To identify PVY-responsive miRNAs, normalized miRNA counts (given in RPM) were compared between PVY and mock samples. miRNAs accumulating beyond a threshold of ±2 log 2 fold change (log 2 fc) in infected vs. healthy plants and at least 1 RPM (approx. 20 raw reads in one of the data set) were considered differentially expressed. miRNA expression profiles were subjected to hierarchical clustering data analysis using the online ClustVis tool (https://biit.cs.ut.ee/clustvis/) [46] and applying the following parameters to the analysis: Clustering distance for rows-Euclidean; Clustering method for rows: complete; Tree ordering for rows: tightest cluster first; and Number of clusters in rows: 30. A heat map was generated to show the changes of each miRNA expression according to time course after PVY C -to infection, and miRNAs showing similar expression pattern were grouped.
Results of miRNA-target (PHAS loci) interactions were used to reveal miRNA triggers of the phasiRNA production. Only 22-nt miRNAs were kept as potential triggers [48,49].
A phylogenetic analysis was performed for the 19 disease-related NLRs and RLP/RLKs genes, selected by generating phasiRNAs and representing targets of 22-nt miRNAs. Only genes with a phasescore higher than two were included in the analysis. Selected NLR/RLK/RLP genes and additional tomato genes active in resistance against fungi, nematodes, and viruses, retrieved from the literature (Table S12) were used to infer a phylogenetic tree. All sequences were aligned using ClustalW within the BioEdit program (http://www.mbio.ncsu.edu/BioEdit/bioedit.html) [50], considering only the C-terminal LRR domain superfamily. The neighbor-joining method [51] was used to build up the phylogenetic tree using the Molecular Evolutionary Genetics Analysis (MEGA) software version X (https://www.megasoftware.net/home) [52] with a bootstrap of 1000 replicates. The tomato Ubiquitin gene sequence (Solyc11g012950) was used as the out-group.
All DEMs were analyzed for functional over-representation in biological pathways with the MapMan software (https://mapman.gabipd.org/) [53]. miRNAs were grouped into MapMan s according to their predicted targets determined in silico and by Degradome-Seq analysis, and using tomato ontology information (for version ITAG3.2) (http://www.gomapman.org/export/current/mapman). Moreover, a Wilcoxon rank sum test was carried out to identify the pathways in which miRNAs underwent coordinated changes [53].

Degradome-Seq Target Validation
Six degradome datasets produced from virus-infected and healthy tomato leaves [30,54] were analyzed with CleaveLand4 (http://sites.psu.edu/axtell/software/cleaveland4/) [55] using all our experimentally-identified miRNAs and the tomato transcriptome sequences (ITAG release 3.2). All identified degradation targets were classified into five categories (0 to IV) as previously described [56]. Only categories with high confidence of cleavage (0, I, II, III) were considered for biological interpretation. The results of miRNA-target (PHAS loci) interactions were also used to confirm miRNA triggers of the phasiRNA production determined in silico.

Validation of miRNA Targets with RLM 5´-RACE
An RNA Ligase-Mediated Rapid Amplification of cDNA Ends (RLM 5'-RACE) protocol [57] was used to validate miRNA cleavage events of selected NLR genes in PVY-infected tomato leaf tissues. Four gene-specific reverse primer pairs and a common forward primer for semi-nested PCR were designed based on the downstream sequence of the predicted miRNA-guided cleavage site of NLR transcripts Cf-9/Solyc01g005870 and R1A-4/Solyc05g008070 using the IDT's PrimerQuest software (Table S11). Four µg of HMW RNA were ligated to a synthetic 5' HD adapter (Table S11) using a T4 RNA ligase (New England Biolabs, Beverly, MA, USA) as described [58]. After purification with the RNA Clean & Concentrator kit™ (Zymo Research, Orange, CA, USA), and the RT reaction was set up using gene-specific reverse primers 5RACE-2_REV (Table S11) and SuperScript ® IV Reverse Transcriptase (Thermo Fisher Scientific). A semi-nested PCR assay was then performed to amplify the cleaved 5 -end using 5' HD adapter-specific forward primer RP-1_FOR in both amplification rounds, while reverse primers were gene-specific 5RACE-2_REV and 5RACE-1_REV in the first and the second amplification round, respectively (Table S11). The first and second PCR rounds were performed in a final volume of 50 µL containing 1.25 U DreamTaq DNA Polymerase, 1X DreamTaq Green Buffer (Thermo Fisher Scientific), 400 nM of each forward and reverse primer, 200 nM of each dNTP, and 2 µL of cDNA (or 1 µL of 1:100 diluted templates from the first-round PCR). The amplification conditions of the first-round PCR consisted of 60 s at 95 • C, followed by 10 cycles of 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 30 s, followed by a final step of 72 • C for 5 min. For the second-round PCR, conditions were 60 s at 95 • C, followed by 35 cycles of 95 • C for 30 s, 60 • C for 30 s, and 72 • C for 20 s, and a final step of 72 • C for 5 min. All PCRs were conducted in a T100™ Thermal Cycler (Bio-Rad).
Amplicons from the second round PCR were separated by 1.2% agarose gel electrophoresis. Bands of the expected size (approx. 180-200 bp) were excised from the agarose gel, purified with QIAEX II ® Gel Extraction Kit (Qiagen, Hilden, Germany) and cloned into chemically-competent E. coli DH5α cells using the pGEM®-T Easy Vector (Promega, Madison, WI, USA), following the manufacturer's protocol. For each cloned amplicon, about 10 clones were selected and sequenced with the M13 reverse primer (Table S12) at Macrogen Europe (Amsterdam, the Netherlands). Sequences were aligned using MultAlin software (http://multalin.toulouse.inra.fr/multalin/) [59] and analyzed to identify specific cleavage sites within homologous groups of sequences.

Data Availability and Retrieval
The sRNA sequencing data herein produced were deposited in NCBI Sequence Read Archive (SRA) database, under the accession number PRJNA563278.

Tomato plants Exhibit Symptom Recovery at Later Stages of PVY C -to Infection
Viral infection was detected at 7 dpi in all PVY C -to inoculated plants by dot blot hybridization (not shown). Starting from 10 to 12 dpi, PVY C -to induced symptoms of mild mosaic, leaf distortion and plant growth reduction. Symptoms reached the highest severity at 21 dpi, when the first set of leaf tissues for RNA extraction was sampled. After nine days (30 dpi), PVY C -to infected plants started showing the signs of disease recovery, and newly-grown leaves showed milder disease symptoms ( Figure S1a). Symptom observations were repeated and confirmed in three independent experiments.

Viral siRNA Accumulation Levels Correlate with Viral RNA Titer
To investigate the correlation between viral RNA titer and production of vsiRNAs in symptomatic leaves (21 dpi) and during the following symptom recovery phase (30 dpi), four sRNA libraries were generated from either healthy (mock-inoculated) or infected (PVY C -to-inoculated) tomato leaf tissues from both time points.
Libraries yielded 20 million reads per sample on average, with 21-24-nt sRNA constituting the predominant size classes (Table 1, in bold). In libraries from mock-inoculated plants, 24-nt sRNAs were largely prevalent (56 and 53% in Mk-21 and Mk-30, respectively), followed by 23-nt (15%) and similar numbers of 22 and 21-nt sRNA (about 10 and 12% in Mk-21 and Mk-30, respectively). The size distribution of host sRNAs was altered in libraries from PVY-inoculated plants, where 21-nt sRNA were predominant (56 and 40% in PVY-21 and PVY-30, respectively) and 24-nt sRNA were reduced to 16 and 27% in PVY-21 and PVY-30, respectively. Also, 21-nt vsiRNAs accounted for 38 and 27% of sRNAs in PVY-infected samples from symptomatic (21 dpi) and recovered leaves (30 dpi), respectively (Table 1). By measuring the accumulation levels of PVY RNA, we found that level of vsiRNAs in PVY-infected plants positively correlated with the viral genomic RNA concentration. PVY-infected symptomatic leaves at 21 dpi exhibited higher vsiRNAs levels and higher virus RNA titers, whereas the vsiRNA abundance was lower as viral RNA load decreased by about 30% in recovered leaves at 30 dpi ( Figure S1b).

Known and Novel miRNAs and Their isomiR Variants Were Identified in Tomato Leaf Tissues
A total of 175 unique miRNAs were identified belonging to 64 known (i.e., reported in miRBase) miRNA families. Out of those 175, 123 matched to previously described tomato miRNAs and 52 to miRNAs found in other plant species. Of the latter 52 known miRNAs, 47 were miRNA-related isomiR forms (i.e. sequence variants), while 5 were mapped to novel tomato MIR genes (Table S1, Figure S2). Using in silico miRNA prediction, 24 novel MIR genes were discovered, coding for 30 novel unannotated canonical miRNAs and 73 previously unannotated novel miRNA variants. Based on the similarity of their predicted precursor sequences with other annotated plant precursor sequences, novel miRNA sequences were assigned to 15 known and 11 novel miRNA families (Tables S1 and S2).
Many studies have shown that isomiRs are functional and can act cooperatively with their canonical miRNAs to target common biological pathways [60]. Hence, we examined the diversity of MIR loci in production of isomiRs and their biological relevance in the tomato-PVY interaction. Our libraries contained a total 587 previously unannotated isomiRs ( Figure S2). Only 19 MIR loci did not generate isomiR variants, while 18 MIR loci produced more than 10 isomiRs, and 4 of them, namely MIR6023, MIR7981f, MIR9474, and MIR10532, were the most diverse, giving rise to more than 20 different miRNA variants. Among novel MIR loci, MIR11 locus generated the highest number of different sequence variants (33) (Table S1). One example of sequence diversity generated from a single MIR locus is MIR6023, encoding canonical sly-miR6023, a well-characterized miRNA which regulates R genes in tomato [4]. Thirty-one sly-miR6023 variants were identified in at least two of the four libraries. By mapping all isomiRs to the pre-miR6023 sequence (miRBase acc. MI0020240), five distinct clusters were identified (Figure 1a). According to the proposed stem-loop structure of the miR6023 precursor, two cluster pairs (Figure 1b, clusters 1-4 and 2-3, respectively) were apparently organized as matching sequences (miRNA/miRNA* duplexes), whereas for cluster 5, a putative miRNA* sequence was not detected (Figure 1b). Five miRNA sequences originating from the miR6023 precursor, one per each isomiR cluster, were quantitatively predominant: the canonical 22-nt sly-miR6023 on the 5p-arm (cluster 1) and its putative miRNA* sequence sly-miR6023.15 (cluster 4); the 22-nt sly-miR6023.28 on the 3p-arm (cluster 3) and its putative miRNA* sequence sly-miR6023.1 (cluster 2); and the 19-nt sly-miR6023.23 on the 3'-terminal region of the 3p-arm (cluster 5) ( Table 2).

Regulation of Known and Novel Tomato mRNA upon PVY C -to Infection
We identified 247 and 219 miRNAs at 21 and 30 dpi respectively that were upregulated by PVY infection. Among those, 192 were upregulated at both analyzed time points during infection. On the other hand, 66 and 11 miRNAs at 21 and 30 dpi respectively were downregulated, with only two miRNAs showing downregulation at both time points ( Figure S3).
As reported in Table 3, 28 out of 192 miRNAs upregulated at both time points were canonical tomato miRNAs already recorded in miRBase, and 13 were canonical, newly-discovered miRNAs from either known or novel miRNA families. Among downregulated canonical miRNAs, seven were downregulated at 21 dpi, whereas their levels remained unchanged at 30 dpi (Table 3). Interestingly, all novel miRNAs retrieved in the present work were regulated in at least one time point: novel-sly-miR4 was downregulated at 21 dpi, novel-sly-miR10 was upregulated at 30 dpi, and all other mature novel miRNAs were upregulated at both time points (Table 3, Table S1). Many isomiRs also showed altered accumulation levels upon PVY infection; among miRNAs altered at 21 dpi, for instance, 81.4% and 89.4% upregulated and downregulated miRNAs, respectively, were isomiRs (Table S1). In some cases, e.g. sly-miR6022, sly-miR6025 and sly-miR10534, the canonical miRNA forms were less abundant than some of their isomiRs in our libraries (Table S1).  Some miRNAs were found to accumulate specifically in virus-infected samples whilst being expressed at null or very low levels in healthy plants. Table S3 lists miRNAs that were found in infected libraries (at both time points) while being not or barely detectable in healthy tomato samples. Fifty-three miRNA sequences belonging to 22 known miRNA families were detected only in PVY-infected libraries, and four miRNAs from novel miRNA families also responded to these criteria. A relevant presence of defense-related miRNA families (e.g. sly-miR396, miR398, miR482, miR6022, miR6023, miR6027) was found in this dataset. In the case of sly-miR6023, significantly, 17 out of 31 sly-miR6023 isomiR forms were specifically associated with PVY infection, indicating that the wide diversification of this miRNA family is, at least in part, a direct consequence of viral infection ( Table 2, Table S3).
In order to identify the miRNAs potentially involved in the regression of viral symptoms and replication at 30 dpi, differences in miRNA expression patterns at the two time points were examined. Hierarchical clustering was applied to the dataset to highlight the changes of miRNA expression profiles according to time course after PVY C -to infection, and miRNAs displaying similar expression patterns were grouped together ( Figure 2, Table S4).
In clusters 3, 5, 9, 11, and 13-30, miRNAs showing similar levels of expression at both time points were grouped. Clear expression differences between 21 and 30 dpi were evident in two series of clusters, with the first series comprising clusters 1, 2, 10, and 12 including miRNAs whose accumulation levels decreased from 21 to 30 dpi. Notably, among the well-characterized miRNAs, sly-miR156 and sly-miR166 were negatively regulated only at 30 dpi, and sly-miR396a, sly-miR396b, and sly-miR397-3p exhibited higher accumulation at 21 dpi and then substantially decreased at 30dpi (Figure 2, Table S4). The second series (clusters 4, 6, 7, and 8) included miRNA with increasing accumulation from 21 to 30 dpi. Defense-related miRNA such as sly-miR482d-3p, sly-miR6023, sly-miR6024, and sly-miR6027-3p were predominantly downregulated at 21 dpi and increased their expression at the following time point, and members of the sly-miR156, sly-miR169, and sly-miR171 families were substantially upregulated at 30 dpi while their levels were not altered at 21 dpi ( Figure 2, Table S4).
Altogether, these results indicated that the effects of PVY C -to infection on miRNA metabolism declined from 21 to 30 dpi, and suggested that a set of miRNAs with changing expression pattern between time points could have an important biological effect in the transition from diseased to recovered conditions.

RT-qPCR Expression Analysis of Known and Novel miRNAs and Their Target mRNAs
To confirm the sRNA-seq results, RT-qPCR expression analysis was performed on a randomly-selected subset of known and novel miRNAs, as well as on their target mRNAs. Thirteen miRNAs were upregulated and two were downregulated at 21 dpi according to both sRNA-seq and RT-qPCR, although fold change figures were different in some cases (Figure 3a). For one out of 16 miRNAs shown in Figure 3a, sly-miR172, RT-qPCR and sRNA-seq results for miRNA expression were discordant.  RT-qPCR also provided additional experimental evidence that 7 novel miRNAs were authentically expressed in tomato tissues. Novel-sly-miR4 was the only novel miRNA detected by sRNA-seq downregulated upon PVY infection, whereas 6 additional novel miRNAs were all significantly upregulated, as observed by sRNA-seq and confirmed by RT-qPCR (Figure 3a).

Target Genes Prediction and Functional Characterization of Differentially Expressed Tomato miRNAs
To better understand the role of miRNA-mediated pathways on host defense responses to PVY infection, we performed a functional analysis based on the known or the putative miRNA-regulated genes function. A total of 5867 unique miRNA/mRNA pairs were predicted by the psRNATarget tool, 1181 of which involved canonical miRNAs and the remaining part isomiRs (Table S5). We further tested whether the accumulation of selected miRNA-regulated transcripts in tomato could be altered as a consequence of PVY C -to ability to interfere with miRNA expression profiles. The comparison of the expression levels of the selected genes from infected and mock-inoculated plants at 21 dpi was made by RT-qPCR using the same leaf samples used for miRNA analysis (Figure 3). The abundance of AGO1a, DCL2a, DCL2d, and NAC1 transcripts was significantly increased upon PVY infection compared with mock-inoculated plants (Figure 3c). AGO1 is a target of sly-miR168a/b, DCL2 family members are targets of sly-miR6026, and NAC1 is a target of sly-miR164a/b. Thus, since the latter miRNAs were all upregulated in virus-infected plants (Table S1), in these four cases, a positive correlation between transcript and cognate miRNA accumulation levels was observed. On the other hand, the expression profiles of DCL1 and Homeobox leucine-zipper protein (HD-Zip III), and of their cognate miRNAs, sly-miR162 and sly-miR166a/b respectively, did not differ significantly upon viral infection, indicating that their gene expression homeostasis was not altered by PVY C -to (Figure 3c, Table S1). An unvaried level of sly-miR6026, as inferred by sRNA-seq data, corresponded to an increased accumulation of its targets DCL2a and DCL2d, whereas DCL2b did not show significant modulation (Figure 3c).
We identified a list of predicted mRNA targets that were also for novel miRNAs (Table S7). We obtained a list of 148 putative targets since, as in the case of known miRNAs, for each individual novel miRNA, more than one (up to 27 for novel-sly-miR7) putative target were bioinformatically retrieved. Interestingly, three mRNAs coding for a Glycine-rich protein, a Disulfide-isomerase-like protein and a Protein phosphatase 2C family protein were identified as being post-transcriptionally regulated by novel-sly-miR1, novel-sly-miR6, and novel-sly-miR10 respectively, an outcome that was additionally supported by the results of the Degradome-seq analysis performed in silico on publicly-available datasets (Tables S7 and S9). Eight out of 15 putative targets of novel_sly-miR10 were NLRs, mitogen-activated protein kinases, or other resistance-related genes. Disease resistance genes were also predicted targets of other novel miRNAs (novel-sly-miR1, -miR4, -miR5, miR7, and -miR8). Finally, novel-sly-miR4, -miR5, -miR7, and -miR9 putatively targeted mRNAs classified in different transcription factor families (e.g. WRKY, MYB, ERF) (Table S7).
To investigate the main biological pathways affected by the miRNA-mediated regulation, a MapMan analysis was performed based upon DEMs detected at 21 and 30 dpi and their predicted targets. The results revealed that miRNA which targets WRKYs and genes involved in jasmonic acid, ethylene, and brassinosteroid biosynthesis and signaling were concordantly upregulated at both time points, whereas miRNAs targeting nucleotide metabolism components were mostly downregulated (Table S6). When we focused on processes involved in plant response to biotic stress [61] (Figure 4), we discovered that several DEMs were targeting mRNAs coding for different immune receptors (NLRs, RLKs), RS proteins such as AGOs and DCLs, proteins maintaining redox homeostasis, proteins involved in biosynthesis and signaling of different phytohormones, calcium signaling, as well as transcription factors belonging to AP2/ERF, MYB, ARF, and bZIP family proteins were upregulated at 21 dpi and 30 dpi (Figure 4, Table S6). By closely inspecting the expression level of miRNAs targeting NLRs and RLKs, we observed that many miRNAs were commonly upregulated at both time points; however, several NLR and RLKs-targeting miRNAs exhibiting downregulation at 21 dpi were weakly downregulated or their level was unaltered at 30 dpi (Figure 4). The same trend was also observed for certain miRNAs targeting ARF and MYB transcription factors, auxin, ethylene, and gibberellin signaling components, as well as mRNAs encoding proteins involved in secondary metabolism, cell wall synthesis and degradation, and calcium signaling, suggesting that their diminished regulation could be important in the recovery process.
Therefore, our data indicate that PVY-infected tomato plants could activate and/or regulate via miRNA pathways several different defense responses mediated by NLRs and RLP/RLK receptors, transcription factors, RNA silencing, and other defense-related metabolic functions. The overall impact of PVY infection on the expression profiles of defense-related miRNAs was higher at 21 dpi and decreased at 30 dpi.

PVY C -to Infection Induces Secondary Phased siRNA Accumulation
miRNAs can potentiate their silencing capability by triggering the production of phasiRNAs from their targets [62]. Over 900 PHAS loci generating 21-or 24-nt phasiRNAs were identified (Table S8). The highest number (535) of PHAS loci was detected in PVY-21 library, and 180 out of these 535 loci were not detected in any of other libraries, whereas 63 PHAS loci were identified only in PVY libraries (Table S8). When inspecting the PHAS loci with the highest degree of phasing, we observed that the majority of them encoded immune receptor proteins harboring NB-ARC and/or LRR domains (Table S8). Next, we integrated 22-nt miRNAs, in silico predictions, and the Degradome-Seq data to identify the triggers for the predicted PHAS loci. MiRNA triggers for 114 and 38 PHAS loci were identified in silico and Degradome-Seq data, respectively (Table S8, Table S9), with a large number of predicted triggers belonging to the miR482, miR6023, miR6024, miR6025, and miR6027 families. For 19 interactions predicted in silico, the cleavage events were also detected by Degradome-Seq analysis (Table S8, Table  S9). Different classes of immune receptors, including a majority of CC-NBS-LRRs (CNLs), but also TIR-NBS-LRRs (TNLs) and RLPs, were targeted for phasiRNA biogenesis in PVY-infected tomato plants ( Figure S4, Table S12).
It is noteworthy that a restricted set of transcripts was predicted to be specifically targeted by isomiRs only. By limiting the search to disease-responsive genes (e.g. NLRs, RLKs, RS-related genes), 9 isomiRs, including 2 isomiRs from novel miRNAs, were specifically associated with a unique degraded transcript (Table S9). Three isomiRs' predicted targets, i.e., Solyc12g100030 (RLP), Solyc01g008800 (TNL), and Solyc11g008540 (DCL2b), were defense-related genes also found as phasiRNA-generating loci (Table S8), which strongly suggests that those genes could be miRNA/isomiR targets in PVY-infected tomato plants. Additionally, an isomiR, sly-miR6023.7, was the only 22-nt long isomiR in that set. Sly-miR6023.7 mapped with a +4 nt shift on the pre-miRNA sequence compared to mature sly-miR6023 and was found exclusively in PVY-infected sRNA libraries ( Figure 1, Table 2). The degradome-seq analysis indicated four potential targets for sly-miR6023.7, including the RLK Solyc09g064670 (Table S9).
To assess whether the virus infection can induce or suppress the phasiRNA production over the time of infection, the accumulation level of sRNAs mapping to PHAS loci was compared between PVY and mock sRNA samples at the two time points. By normalizing the four libraries to RPM and selecting genomic loci showing in-phase accumulation of secondary sRNAs (phase score >1), 241 and 140 loci gave origin to higher accumulation of sRNAs in PVY-infected samples than in mock samples at 21 and 30 dpi, respectively, displaying at least 2-fold higher sRNA accumulation. Of those PHAS loci, 54 were found to be upregulated at both time points, but 10 loci displayed an opposite regulation profile (i.e. phasiRNAs over-accumulated at 21 dpi and decreased accumulation at 30 dpi in PVY vs. healthy samples) ( Figure S5). Eighty-nine and 84 loci showed higher phasiRNA accumulation in mock versus infected samples at 21 and 30 dpi respectively, with 22 loci being downregulated by PVY infection at both time points ( Figure S5, Table S8).
Overall, these data indicate that cleavage by 22-nt miRNAs may trigger the production of phased secondary 21-nt siRNAs as a mechanism of fine regulation of defense-related gene expression that is dynamic over time, and that isomiRs are likely involved in the same mechanism in tomato plants infected by PVY C -to.

miRNAs Mediate Cleavage of Transcripts Encoding NLRs and RLPs
To analyze in more detail the activity of miRNAs, the phasiRNAs biogenesis and the post-transcriptional regulation of RLP and NLR gene expression during the PVY infection, we identified two PHAS loci, i.e., Solyc01g005870 and Solyc05g008070.
According to the NCBI Solanum lycopersicum Annotation Release 103, Solyc01g005870 codes for the receptor-like protein Cf-9 [63] (NCBI Ref. Seq.: XM_004228576.4), and Solyc05g008070 codes for a CNL, the putative late blight resistance protein homolog R1A-4 [64] (NCBI Ref. Seq.: XM_010322406). Cf-9 was chosen for further analysis since, among sRNA reads mapping to its sequence, we found a very abundant (>3000 raw reads) 21-nt antisense sRNA in the PVY-21 library (Figure 5c). Searches in the Tomato Functional Genomics Database for small RNAs [65] identified a homologous miRNA named M00093 ( Figure S6a). Although this miRNA family was neither retrieved in the miRBase nor identified as a novel miRNA by our own analysis, it was found in the literature, described as sly-miR00093 [66] or as tomato novel miRNA pc-27 [67]. Sly-miR00093 was predicted to target at least other seven RLPs which were closely related to Cf-9 ( Figure S6b). The same tomato transcript Cf-9 was predicted also as the target of 22-nt-long sly-miR6023 and 21-nt-long sly-miR6022 by in silico and degradome-seq prediction analyses, respectively (Figure 5a, Table S8, Table S9). The second transcript, i.e., R1A-4, was the NLR gene with the third highest secondary sRNA accumulation level, and was predicted as the target of 22-nt sly-miR482b/c, sly-miR6024, sly-miR6026, and sly-miR6027-3p (Figure 5b, Table S8). At 21 dpi, a redundant set of 275 phasiRNA (81 RPM, phasescore = 4.4) mapped to Cf-9, while 1548 phasiRNA (456 RPM, phasescore = 25.2) mapped to R1A-4 (Table S8). miRNAs that were predicted to regulate NLR genes post-transcriptionally were quantified by RT-qPCR to further confirm the effects of viral infection on their accumulation levels. A general trend of upregulation was found upon PVY C -to infection at 21 dpi, as demonstrated for sly-miR6023, that was significantly overexpressed by 5.32 log 2 fc in infected vs. healthy tomato plants, for sly-miR6024 (log 2 fc = 2.18), for sly-miR6027-3p (log 2 fc = 2), and for sly-miR00093 (log 2 fc = 2.9) (Figure 3b). Divergence between sRNA-seq and RT-qPCR quantification could depend upon the fact that template RNA was extracted from two separate experiments, and could also be explained by the fact that these specific miRNAs, forward qPCR primers, can potentially detect and quantify more than one isomiR variant along with the canonical miRNA sequence (Table S10).
To further verify the possibility that tomato miR00093, miR6022, and miR6023 concurred in the post-transcriptional regulation of Cf-9, we investigated the presence of specific cleavage products in PVY-infected tomato leaf tissues using RLM 5'-RACE. The predicted miRNA-guided cleavage products, corresponding to the position between nucleotides 10th and 11th of the miRNA sequences, were correctly amplified and sequenced in 8/8 PCR fragments for 22-nt sly-miR6023 and 7/11 PCR fragments for 21-nt sly-miR00093, and in the latter case, an additional cleavage site between nucleotides 11th and 12th of sly-miR00093 was identified in 3/11 PCR products (Figure 5e), indicating that both one-hit (22-nt miRNA) and two-hit (two 21-nt miRNA triggers) mechanisms were employed for the generation of phasiRNAs from Cf-9. In addition, the degradome-Seq data and results by Bai and colleagues [54] showed that Cf-9 can be also targeted by miR6022 (Figure 5a). Since our RLM 5'-RACE experiment failed to detect any cleaved mRNA products at this specific site in PVY-infected tomato tissues (Figure 5e), our results suggest that the regulation of Cf-9 by individual miRNAs might be disease-specific. , as correctly predicted for the sly-miR6024-guided activity. In 12/12 cases, PCR fragments from RLM 5'-RACE evidenced a cleaved 5'-end which was compatible with this prediction (Figure 5f). Despite the demonstrated miRNA-guided cleavage activity on Cf-9 and R1A-4, RT-qPCR analysis showed that both genes were overexpressed upon PVY C -to infection (Figure 3c). Therefore, the post-transcriptional expression regulation of these genes by miRNAs did not result in the straightforward suppression of gene accumulation in PVY-infected vs. healthy plants, but rather, in a sophisticated miRNA-triggered fine-tuning of gene expression, also involving the induction of phasiRNAs production from their targets.

Discussion
In this study, we used a high-throughput sRNA-Seq approach to investigate the effects of PVY infection on sRNA metabolism and on regulation of immune responses in tomato. PVY infection in tomato induced a huge reprogramming of miRNA accumulation levels. The general effect of PVY C -to in infected tomato tissues was the upregulation of miRNAs. Similarly, a large predominance of upregulated miRNAs was also observed upon PVY infection in other hosts [13,28]. In this work, more than 200 miRNA sequences were shown to be consistently overexpressed at 21 and 30 dpi. Those can be considered as the "core" set of PVY-responsive miRNAs. Induction of miRNA expression is frequently regulated by transcription factors acting in trans on MIR gene-promoter sequences [68]. In virus-infected plants, an additional mechanism of miRNA accumulation depends on the presence of VSR proteins and their ability to bind sRNAs. HC-Pro, a VSR of PVY and other potyviruses, has been proven to bind miRNA/miRNA*, preventing miRNA loading into RISC and determining the abnormal accumulation of inactive miRNA and miRNA* forms [69][70][71]. The large number of over-accumulating miRNA* sequences found in our dataset indicates that the HC-Pro activity may be a major cause of alteration of miRNA metabolism also in PVY C -to-infected tomato plants. In fact, this observation implies that a conspicuous part of detected miRNA sequences is inactive as miRNA/miRNA* duplexes, and therefore, impaired in their gene expression regulatory function [69,70,72]. Accordingly, our work shows the occurrence in PVY-infected tomato plants of cases of: i) anomalous accumulation of miRNA* sequences, not detected in healthy plants (e.g. novel-sly-miR11*, miR164a-3p, miR168a-3p, miR171b-3p, miR396a-3p, etc.) and ii) positive correlation between accumulation levels of miRNAs and their target mRNAs.
Besides the quantitative effects on miRNA accumulation levels, PVY infection correlated with a remarkable increase of sequence variation within individual miRNA families. Although frequently underestimated, some authors warned about the importance of miRNA diversification and isomiR production in plant biology, since new isomiR biogenesis is proposed to strengthen or to fine-tune RS-based antiviral defenses [73,74]. Several enzymatic mechanisms have been shown to contribute to isomiR biogenesis (reviewed in [74]), some of which could be activated by virus infections possibly interfering with host immune responses. In rice for instance, severe infection caused by Rice stripe virus, but not mild infection caused by Rice dwarf virus, induced the accumulation of novel isomiR sequences from conserved pre-miRNAs. Those novel isomiRs were shown to regulate the expression of their predicted target mRNAs [75,76]. Novel miRNA variants were also observed by Hu and colleagues [77] in Arabidopsis plants infected with Oilseed rape mosaic virus. Our data indicate that many isomiRs are produced as a consequence of PVY C -to infection, and suggest that PVY-associated isomiRs such as sly-miR396a.2, sly-miR397.5, sly-miR6022.15, and sly-miR6023.7 have the potential to regulate the expression of RLPs, NLRs and RS-encoding genes by inducing mRNA cleavage and phasiRNA accumulation. To the best of our knowledge, our results represent the first example of the production of several different isomiRs strictly associated with the infection of a virus in a solanaceous host. Future research on this intriguing subject will expand our knowledge of the multifaceted role of miRNAs in plant defenses.
MapMan analysis indicated that the altered expression of many miRNAs potentially affects the signaling and defense pathways involved in the plant response to biotic stresses. We showed that at 21 dpi during the acute phase of PVY infection, many defense gene-targeting miRNAs are downregulated, presumably allowing the increased expression of their respective targets that are deployed by host cells to counteract the pathogen. At 30 dpi, when decreased virus replication and alleviated symptoms were observed, the expression of defense genes could be relieved through the reduced alteration of miRNA levels ( Figure 4). It is well known that many miRNA targets are transcription factors and genes regulating phytohormones biosynthesis, implying that plant miRNAs are master regulators of many physiological functions [78]. Following viral infection in plants, the regulation of transcription and hormone synthesis/signaling are often found among altered functional categories for miRNA targets [13,79]. Since many growth and development processes are controlled by miRNA-regulated factors, it has been proposed that the alteration of miRNA metabolism may represent one of the major mechanisms through which pathogenic viruses induce symptoms on infected plants [22,30,33,80].
Data from the present work also show significant changes in expression profiles of miRNAs controlling various classes of R genes including NLRs (TIR-NBS-LRR and CC-NBS-LRR), RLPs, and RLKs. Upon PVY infection, miR482, miR6022, miR6023, miR6024, and miR6027 families, known for targeting R genes in tomato [81][82][83], were found to be either up-or down-regulated in our libraries. Sly-miR6023 and sly-miR6024 were here validated as post-transcriptional regulators of Cf-9 and R1A-4, respectively. An additional miRNA species, referred to as sly-miR00093, was also experimentally validated as another miRNA targeting Cf-9 for post-transcriptional silencing. We provide evidence that miRNA-guided cleavage of R genes is a specific event which strictly depends on the host-pathogen combination: in fact, whereas R1A-4 was proved to be a target for miR6024 in both PVY-infected (this work) and Tomato yellow leaf curl Sardinia virus (TYLCSV)-infected tomato cells [84], we failed to confirm the Cf-9 cleavage operated by sly-miR6022, in contrast to what observed in Tomato yellow leaf curl virus-infected plants [54]. The low levels of mature sly-miR6022 in our libraries might explain the lack of cleavage, although sequence variations at miRNA recognition sites among host genotypes should also be taken into account as a possible reason for impaired miRNA activity [84].
The consequence of miRNA activity on target R genes is not only limited to mRNA cleavage, but can also lead to biogenesis of secondary, RDR6-dependent phasiRNAs [81,85]. PHAS/R loci discovery prompted scientists to suggest a mechanism by which R genes accumulation could be controlled by miRNAs by triggering secondary siRNA production finalized to potentiate post-transcriptional repression of resistance gene expression [86,87]. miRNA-mediated down-regulation of NLRs would allow the plant to balance in a dynamic process the metabolic costs of accumulating constitutively expressed NLR proteins and the benefits of effective protection against pests and pathogens [88,89]. The mechanism generates a "cascade" effect, since phasiRNAs would, in turn, amplify the defense response either by acting post-transcriptionally in cis, by silencing more R genes in trans [81], or by affecting the disease phenotype targeting other non-R genes [32,84]. Our data further corroborate the notion of phasiRNA-mediated defense mechanisms. In fact, the NLR Solyc02g036270, the PHAS/R locus with the highest degree of phasing in our study, and other loci such as Solyc04g005540 and R1A-4, were also reported as secondary siRNA-producing loci in previous studies, [81,84,90]. However, the present work also provides new evidence that sheds light on this mechanism. Importantly, we show that a highly-significant phasiRNA biogenesis activity occurred at 21 dpi upon PVY infection, and therefore, correlated with the infection peak, whereas a more limited phasiRNA accumulation was observed at 30 dpi and in libraries from healthy plants. The additional observation that, at least for a limited number of PHAS loci, phasiRNA biogenesis can significantly vary over time from overexpression to suppression at an individual locus, suggests the existence of a highly-dynamic mechanism oriented towards the tight control of gene expression, which is particularly important for NLR/RLP/RLK genes. It is noteworthy that a reduced phasiRNA accumulation was observed at 30 dpi in comparison to 21 dpi, when infected leaf tissues exhibited a regression of symptoms. According to previous works, upon pathogen attack, plants respond with an enhanced R gene expression driven by suppressing the miRNA-mediated silencing pathway. The PHAS/R regulatory mechanism would be then activated for controlling excessive R gene expression [12,81]. More research is needed to refine a mechanistic model to explain phasiRNAs role in host recovery from PVY-induced diseases.
Tomato genes encoding enzymes of the RS pathway were also found as phasiRNA-producing loci. It has been demonstrated that sly-miR6026 targets tomato DCL2a, DCL2b, and DCL2d in vivo, and that its cleavage is responsible for phasiRNA production at these loci [91]. In PVY-infected tomato, we could detect phasiRNA mapping on DCL2a and DCL2b, but not on DCL2d. Conversely, DCL2a and DCL2b did not show upregulation upon PVY infection, whereas DCL2c and DCL2d significantly increased. This observation suggests that DCL2 family members have redundant roles, and although Wang and colleagues (2018) indicated DCL2a/b as the major antiviral DCL2 forms, their post-transcriptional silencing could be replaced in infected tomato cells by increased levels of DCL2c/d. In fact, in Tobacco mosaic virus-infected dcl2ab tomato double mutants, the 22-nt viral sRNAs were less abundant but not absent, indicating residual DCL2 activity [91]. Apparently, HC-Pro, a VSR of PVY, suppresses RS in infected tomato plants without inhibiting the following process of secondary siRNA biogenesis where RDR6/SGS3 and DCL4 guide dsRNA synthesis and degradation, respectively. In contrast, TYLCSV RS suppressor, the V2 protein, was found to bind SGS3 and block antiviral silencing at the stage of single-to dsRNA polymerization operated by the RDR6/SGS3 complex. Probably for this reason, the production of phasiRNAs was abolished upon TYLCSV infection, and was only reported to occur in healthy tomato plants [84].

Conclusions
The work here presented confirms the existence of a wide range of endogenous sRNA-based responses encompassing quantitative miRNA alteration, qualitative diversification within miRNA families, and miRNA-driven secondary sRNA accumulation. From our investigation on host recovery from PVY-induced disease, we were able to identify a set of tomato sRNAs showing changes in accumulation levels over time, between 21 and 30 dpi. Although the molecular mechanisms leading to disease symptom suppression is still unknown, the intriguing possibility that recovery from viral infections may be supported by the direct involvement of miRNAs and other endogenous sRNA has been already suggested by other authors [92] and needs to be further explored. Suggestions of a putative role of these mechanisms in plant defenses and immunity are valuable, since they will provide the basis to reveal further plant-virus interactions, and will contribute to the development of control measures against viral diseases.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4915/11/12/1100/s1. Figure S1. PVYC-to-infected tomato plants exhibit disease recovery at 30 dpi. Figure S2. Classification of miRNA and isomiR sequences found in libraries from mock-inoculated and PVYC-to-infected tomato leaf tissues. Figure  S3. Venn diagrams showing the overlap of DEMs between 21 and 30 dpi. Figure S4. Phylogenetic tree of LRR-domain-containing tomato NLR or RLP coding sequences. Figure S5. Venn diagrams showing the overlap of PHAS loci showing increase or suppression of phasiRNAs accumulation in PVY-infected vs. mock-inoculated plants at 21 and 30 dpi. Figure S6. Characterization of sly-miR00093. Table S1. miRNAs identified in PVY-and mock-inoculated tomato samples together with their expression levels. Table S2. Novel MIR loci and miRNAs identified in PVY-and mock-inoculated tomato plants. Table S3. List of miRNAs and isomiRs specifically associated with PVY-infected plants. Table S4. Clustering results based on the expression profiles of mature miRNAs and isomiRs upon PVYC-to infection at 21 and 30 dpi. Table S5. Predicted targets of miRNAs by in silico approach. Table S6. MapMan analysis results. Table S7. Predicted targets of novel miRNAs by in silico approach. Table S8. Identified PHAS loci. Table S9. Bioinformatic validation of mRNA targets of tomato miRNAs by Degradome-Seq. Table S10. Quantification of miR6023, miR6024 and miR6027 by sRNA-Seq and RT-qPCR. Table S11. List of primers used in RT-qPCR and 5 RACE experiments in this work. Table S12. List of the 20 disease-related genes encoding NLRs and transmembrane RLP/RLK receptors, selected by the capability to generate phasiRNAs.

Conflicts of Interest:
The authors declare no conflict of interest.