Tolerance to PVY infection in potato is conditioned by perturbation of small RNA– gibberellin regulatory network

Potato virus Y (PVY) is the most economically important potato viral pathogen worldwide. We aimed at unraveling the roles of small RNAs (sRNAs) in the complex immune signaling network controlling the establishment of tolerant response of potato cv. Désirée to the virus. We identified and quantified the endogenous miRNAs and phasiRNAs as well as virus-derived sRNAs. We constructed a sRNA regulatory network connecting sRNAs and their targets identified to link sRNA level responses to physiological processes. We discovered an interesting novel sRNAs – gibberellin (GA) regulatory circuit being activated as early as 3 days post inoculation, before viral multiplication can be detected. Increased levels of miR167 and phasiRNA931 were reflected in decreased levels of transcripts involved in GA biosynthesis. Moreover, decreased concentration of GA confirmed this regulation. We have additionally showed that this regulation is salicylic acid dependent as the response of sRNA network was attenuated in salicylic acid-depleted transgenic counterpart NahG-Désirée expressing severe disease symptoms. Moreover, we found that differentially expressed miR6022, which targets leucine-rich-repeat receptor-like kinases, likely links GA signaling with immune responses. Besides downregulation of GA signaling, regulation of several other parts of sRNA network in tolerant Désirée revealed striking similarities to responses observed in mutualistic symbiotic interactions. The intertwining of different regulatory networks revealed here shows how developmental signaling, symptomology and stress signaling are balanced.

Recent findings revealed that endogenous RNA silencing mediated by microRNAs (miRNAs) and small interfering RNAs (siRNAs) could play important roles in plant immunity [13]. These 18-24-nt long noncoding sRNAs are able to negatively regulate gene expression by binding to the specific mRNA targets which leads to either promoting their degradation, inhibiting their translation, or suppressing transcription by epigenetic modification [8,14]. The endogenous RNA silencing can be amplified by the production of secondary phased siRNAs (phasiRNAs), triggered by 22-nt miRNAs/siRNAs [15,16]. phasiRNAs are generated in phase relative to positions of the miRNA cleavage site, can be produced from both coding and non-coding transcript (PHAS loci) and are able to target transcripts not only in trans but also their PHAS loci of origin in cis and thus additionally contribute to the autoregulation [17].
miRNAs have been associated with defense responses against several pathogens [18]. Arabidopsis miR393 was the first plant miRNA reported to play a key role in antibacterial immunity by repressing auxin signaling [19]. Recently, several studies have uncovered the miRNA-mediated silencing of immune receptor gene (R-gene) transcripts. Infection by pathogens e.g. viruses or bacteria, relieves the silencing, leading to the accumulation of R proteins and activation of immune responses [20][21][22].
This growing body of evidence suggests that sRNAs are integral components of plant immunity.
However, none of the studies performed so far investigated the sRNA regulatory network in potatovirus interaction at the systems level linking it to transcriptional regulation. The aim of this study was to investigate sRNAs' role in establishment of the tolerant response of potato to PVY NTN , hence we have studied response in the early stage of viral infection, before the viral multiplication can be detected.
Employing high-throughput sequencing technology, we characterized and compared the sRNA expression patterns between PVY-infected and healthy tolerant potato plants. In addition, this information was linked with expression profiles of their target transcripts identified by in silico and Degradome-Seq and used for sRNA regulatory network construction. Besides the already described regulation of R-gene transcripts we have discovered a previously undescribed connection between sRNAs and gibberellin (GA) biosynthesis representing an important link between defensive and developmental signaling pathways which was confirmed by hormonal content measurements.
Additionally we showed that response of the discovered sRNA network was attenuated in susceptible, SA deficient NahG-Désirée, confirming the crucial role of SA in mediating responses leading towards establishment of the tolerance against PVY infection. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Novel endogenous sRNAs identified in potato leaves
We identified 245 different previously described miRNAs (including 38 miRNA variants; isomiRs), belonging to 95 miRNA families in control and PVY NTN -infected leaves of cv. Désirée using sRNA-Seq (S1 Fig, S1 Dataset). In addition, 141 novel miRNAs were detected, of those 12 were coded by multiple MIR loci. Novel miRNA sequences were assigned to 123 novel miRNA families (S1 and S2 Datasets).
When assessing the miRNA regulatory network, the amplification of silencing through phasiRNA biogenesis was also considered. In total, more than 4000 PHAS loci were predicted, coding for 2558 phasiRNAs. 124 loci were located on protein-coding regions of genes, with the majority encoding NBS-LRRs (nucleotide binding site-leucine rich repeat proteins) and LRR-RLKs (leucine rich repeatreceptor-like kinases) (S3 and S4 Datasets). Interestingly, several novel PHAS loci were discovered in coding regions of genes associated with stress signaling, such as HSP70, superoxide dismutase and auxin signaling F-box proteins (see S1 Table for the full list of genes used in this study together with their descriptions and corresponding IDs). We also confirmed that StAGO1 is a PHAS locus, as previously described in Arabidopsis [23] (S3 Dataset).
We further compared miRNA expression profiles of PVY NTN -infected versus mock inoculated leaves in early stages of virus infection (3 dpi, 1 day before detectable viral multiplication; [5,7]). In total, 57 unique miRNAs were found to be significantly differentially expressed in early stages of PVY NTN infection (3 dpi) of Désirée plants. Virus infection predominantly caused an increase in miRNAs levels (Fig 1, S1 Dataset). Additionally, we identified 35 phasiRNAs as differentially regulated 3 dpi, mainly originating from noncoding PHAS and NBS-LRR loci (Fig 1, S4 Dataset).
To validate the obtained sRNA-Seq results, abundance of six differentially expressed miRNAs was analyzed by stem-loop RT-qPCR. As shown in S2 Fig, all sRNA-Seq differential expression results were confirmed. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.  In previous studies, a plethora of potato miRNA/phasiRNAs have been shown to be differentially expressed following pathogen infection. However, the biological relevance of these differences remains largely unknown. To translate the data obtained on sRNA level into changes in physiological processes, we performed in silico sRNA target prediction, both at the levels of translational inhibition and target cleavage (S5 Dataset). Additionally, the predictions of target cleavage were experimentally validated by Degradome-Seq (S6 Dataset). Based on this information we constructed a potato sRNA regulatory network connecting miRNAs with phasiRNAs and their targets (Supplemental Online Files 1 and 2).
This revealed several already known and many novel connections related to the plant immune signaling (see example in S3 Fig; S5 and S6 Datasets). Interesting targets of differentially expressed miRNAs not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Differentially expressed miRNAs regulate various immune receptor genes in tolerant response
Analysis of the differentially expressed miRNAs and phasiRNAs together with the levels of the target transcripts (data published in [7]) revealed the presence of many known, as well as novel, regulatory cascades involving NBS-LRRs. Several NBS-LRRs were predicted to be targeted by miR482, miR6024 and miR6027 family members (S5 and S6 Datasets). Moreover, NBS-LRRs are regulated also by phasiRNAs, where most phasiRNAs have multiple NBS-LRR targets (S4 Fig, S5 and S6 Datasets) due to the shared conserved P-loop or Walker A motif [21]. In all of the previous studies miR482 family members were downregulated following pathogen infection [21,26,27]. In our study, however, the one regulated member of the miR482 family (miR482e-5p) targeting NBS-LRR transcripts was upregulated following PVY NTN infection (S1 Dataset), similarly as observed in establishment of mutualistic symbiosis in soybean roots [28]. Moreover, several miRNAs that were upregulated in response to PVY NTN in cv. Désirée, such as miR164, miR167, miR171, miR319, miR390, miR393, miR397 and miR398 have also been reported to regulate nodulation in different legume species (S1 Dataset; [29][30][31][32]. In addition to NBS-LRR proteins, LRR-RLKs are also important mediators of immune signaling cascades. We have identified a novel miRNA-LRR-RLK interaction in which miR6022 levels decrease in response to PVY NTN infection in cv. Désirée, which is further linked to upregulation of its predicted target genes encoding LRR-RLKs (S4A Fig).

Several PVY NTN -derived siRNAs trigger degradation of host transcripts
The primary plant defense mechanism against invading viruses is RNA silencing involving the production of vsiRNAs. The population of vsiRNAs detected in the infected samples consisted of more than 46 000 unique sequences of 20-24 nt in length (S7 Dataset). In order to take into account the unlikely possibility that PVY NTN produces its own miRNAs, we first ran the miRNA prediction pipeline on vsiRNAs and the viral genome. However, we found no sequence that would fullfill the criteria for a viral miRNA. Subsequently, we searched for potential host transcripts targeted by vsiRNAs in our experimentally validated target degradation dataset (S8 Dataset). We found that vsiRNAs are indeed not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint able to target multiple potato transcripts, among them mRNAs coding for immune receptor proteins, various transcription factors and proteins involved in hormonal signaling pathways (S8 Dataset). For example, several vsiRNAs were detected with the confirmed ability to guide the cleavage of transcripts involved in auxin signaling, transcripts encoding IAA-amino acid hydrolases (StILR1 and StIAR3), Aux/IAA transcriptional repressors (StIAAs) and the transcription factor StARF2.

sRNA-mediated downregulation of GA biosynthesis genes is reflected in lower GA3 levels
Interestingly, we found that GA biosynthesis and downstream signaling are targeted by a sRNAmediated regulatory network and that the changes in sRNA levels following PVY NTN infection are reflected also in the changes of their target transcripts levels (Fig 4A, S5 Dataset). GA20-oxidase (GA20ox) and GA3-oxidase (GA3ox) are enzymes that catalyze the last steps in the formation of bioactive GAs [33]. We found that in Désirée plants the StGA20ox transcript is regulated by miR167e-3p (Fig 4A). An additional layer of GA biosynthesis regulation is represented by increased production of phasiRNA931, which promotes cleavage of the StGA3ox transcript. The transcriptomics results support these interactions as the targeted transcripts are significantly downregulated in Désirée upon PVY NTN infection (Fig 4A, S5 Dataset). Additionally, vsiRNAs were found to target transcripts encoding two enzymes involved in GA biosynthesis StGA1 and StGA20ox (Fig 4A, S8 Dataset). One also has to note that all of the miRNAs/phasiRNAs discovered to be involved in sRNA-GA biosynthesis regulation have so far not been identified in Arabidopsis and among them only miR167 was also discovered in tomato [34,35].
Downstream GA signaling is also targeted by sRNAs in the potato -PVY interaction on multiple levels.
The four miR319 family members were predicted to cleave the transcript encoding StMYB33, a GAMYB-like transcription factor orthologue [36], whereas phasiRNA931 targets a potato orthologue of a DELLA protein, a GA-signaling repressor (Fig 4A).
Such interconnectedness between plant defense related miRNA/phasiRNA signaling and GA signaling has not been previously identified and may represent a link between defensive and developmental signaling. Thus, we decided to evaluate functionally these results by measuring the concentrations of a set of plant hormones. As predicted by the sRNA-target transcript analyses, we detected a reduced level of GA3 in PVY NTN infected Désirée plants (Fig 4B). The levels of SA, jasmonic acid (JA), the JA precursor 12-oxophytodienoic acid (OPDA), abscisic acid (ABA) and indole-3-acetic acid (IAA) remained unchanged at 3 dpi (Fig 4B, S9 Dataset). To inspect if GA deficiency has any impact on plant growth, plants height was monitored till 21 dpi. No differences were observed between all four studied groups of plants.

The miRNA regulatory network and transcriptional regulation are tightly interconnected
Given the critical role of miRNAs in gene regulation, cis-regulatory elements of differentially expressed MIR genes involved in R-gene regulation and GA signaling were investigated. Interestingly, GAMYB binding sites were detected in the promoters of the MIR6022 and MIR482f genes (S10 Dataset).
Moreover, these genes harbor WRKY8/28/48 binding sites, while MIR319a harbors a general WRKY binding W-box regulatory element. Additionally, cis-acting elements involved in SA and JA responsiveness were identified in the promoters of the MIR482f gene. In all three analyzed promoter sequences for NAC transcription factor binding sites, potentially implicated in defense signaling, were detected (S10 Dataset). The promoter of MIR167e is only partially assembled in the current version of the potato genome; thus, the promoter analysis was performed only for the first 80 nt upstream of the predicted hairpin precursor sequence (S10 Dataset).

sRNA response is attenuated in susceptible SA depleted plants following PVY NTN infection
As the link between repression of GA signaling and disease symptoms severity was already established in Arabidopsis [38] we have further investigated whether the activity of discovered sRNA -GA circuit is acting similarly in susceptible potato -PVY NTN interaction. To keep the same genetic background we have used a transgenic counterpart of Désirée that is depleted in SA accumulation, NahG-Désirée. We have shown previously that in NahG-Désirée plants the pronounced disease symptoms appeared both not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint on the inoculated and systemic leaves [5]. Furthermore, the viral multiplication was detected 1 day earlier than in non-transgenic Désirée plants (at 4 dpi), while the final concentrations of the virus were not significantly higher [5]. Here, we performed the analysis of sRNA response as well as the measurements of hormonal concentrations in interaction of this susceptible genotype with the virus.
Interestingly, we found that the overall response of sRNAs was attenuated in NahG-Désirée at 3 dpi. In NahG-Désirée only 20 miRNAs were differentially expressed, with the majority showing a lower degree of induction than in Désirée plants (Figs 1 and 2, S1 Dataset).
Inspecting specifically the discovered link between sRNAs regulation, GA signaling and immune signaling in the sRNA and transcriptional datasets we observed that the responses of miR167e and miR319a are diminished in NahG-Désirée plants (S1 Dataset). Interestingly, although the phasiRNA931 was also upregulated in NahG-Désirée plants, albeit to a lower extent, this was not translated into downregulation at the target transcript level (Fig 4C). Adding to the significance of this finding, in NahG-Désirée plants the level of bioactive GA was not significantly different in infected leaves (Fig 4B). We also compared sRNA-mediated responses in susceptible NahG-Désirée with those observed in tolerant plants of Désirée, which resembled responses in establishment of mutualistic symbioses. Upregulated miR482e targeting NBS-LRR transcripts was observed in both genotypes indicating that this response is not relevant for symptoms development nor is SA-dependent. On the other hand, the relieved silencing of LRR-RLKs by miR6022 that is predicted to potentiate the defence responses and which is also linked with GAMYBs, being involved in symptoms development, was however absent in NahG-Désirée plants (S4B Fig). Similarly, the responses of other miRNAs involved in mutualistic symbiosis, such as miR164, miR167, miR171, miR319, miR390, miR393, miR397, miR398, were diminished in NahG plants (S1 Dataset).
To evaluate the direct effect of the SA deficiency in NahG plants, independently of the viral infection, we also compared the sRNA expression profiles in mock-inoculated leaves of non-transgenic and SAdepleted Désirée (Fig 1). We found 37 miRNAs regulated by SA. It seems that SA in the normal growing conditions generally causes the reduction in the level of miRNAs as majority (28) of miRNAs were detected at significantly higher levels in NahG-Désirée plants (S1 Dataset). When we similarly compared the expression profile of transcripts between the two genotypes, we detected that the most strongly induced by SA signaling are notably different WRKY transcription factors ( Table 1, S11 and S12 Datasets), among them orthologues of Arabidopsis WRKY70, which was already shown to be positively regulated by SA [39]. As miR319a and miR482f promotors are regulated by WRKY transcription factors, we discovered a direct link between SA signaling and miRNA regulatory network in potato. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Discussion
We hypothesized that fine-tuned regulation of subsets of genes involved in defensive signaling can interfere with developmental signaling, which could explain decreased symptom severity in plants expressing tolerance to virus infection. The sRNA's have proven to be an important level for precise regulation of several developmental processes. We here show that the integration of sRNA network and transcriptional regulation is also crucial in entanglement of stress responses and developmental processes.
Integration of the sRNA regulatory network with sRNA expression and expression profiles of their targets confirmed many known, but also revealed several novel regulatory circuits associated with immunity regulation and hormone signaling (Figs 3 and 4, S3 and S4 Figs, S5 and S6 Datasets). When plants are exposed to pathogens, NBS-LRRs and LRR-RLKs are the key players in sensing and transducing stress signals [40]. Viral suppressors of silencing can release the tight control of R-gene silencing by sRNAs and activate immune responses in plants [20,21]. Our study investigated regulatory processes occurring early after infection, before virus concentration significantly increased thus the effects we detected were not caused by extensive HC-Pro or any other viral protein accumulation. Even so, we have detected diverse regulation of NBS-LRRs and LRR-RLKs and their targeting miRNAs as expected according to their specialized roles (S4 Fig). Upregulation of miR482e (S1 Dataset) and resulting downregulation of its NBS-LRR target transcript was shown to enhance the sensitivity of potato to Verticillium dahliae infection [27]. On the other hand, upregulation of miR482e was also observed to promote soybean nodulation [28] and could in our system also be part of the network allowing the establishment of commensalistic symbiosis (in plant-virus interaction termed tolerance) between potato not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint and PVY. Similarity of response between mutualistic symbiosis in soybean and tolerance in potato was also shown by the profile of several other miRNAs (S1 Dataset; [29][30][31][32]). This suggests a similar modulation of immune response and physiology occurs in both mutualistic and disease tolerant (commensalistic) interactions. In tolerance, plants may adopt some processes resembling mutualistic symbiosis to control plant response and minimize severe disease symptoms allowing non-hindered development of the plant and at the same time multiplication of the virus.
Phytohormones modulate plant defense responses against various biotic and abiotic stresses as well as plant growth and development [41]. Till now, several miRNAs were confirmed to participate in this complex network, mainly in connection to repression of auxin signaling [19]. A similarly complex miR393/miR396/phasiTIR auxin signaling network was identified in this study, yet showing links also to other hormonal signaling pathways (Fig 3, S3 Fig).
Most notable is, however, the novel link between sRNA regulatory network and GA signaling. Biotic stress was shown to repress GA signaling pathways [42]. Here, we show that GA biosynthesis and signaling are post-transcriptionally regulated via multiple miRNAs, phasiRNAs as well as vsiRNAs in potato leaves following infection with PVY NTN (Fig 4A). The effect of this regulatory circuit was confirmed by reduced bioactive GA level in the tolerant Désirée plants (Fig 4B). This reduction was however not reflected in decreased plant growth and was thus most probably transient and localized in nature. sRNAs were shown to regulate GA biosynthesis in potato leaves following PVY NTN infection in three crucial steps, the first step after branching point of GA biosynthesis and the general carotenoid pathway and the last two steps in production of bioactive GA (Fig 4A). In other plant species, GA biosynthesis was shown to be indirectly targeted by miRNAs regulating the activity of the corresponding transcription factors (miR319-TCP14-GA2ox/GA20ox; miR393-GRF2-GA3ox/GA20ox) [25], while direct interactions were to our knowledge not yet reported. Also of note, the sRNAs regulating GA biosynthesis identified here were not identified in Arabidopsis and seem to be Solanaceae specific.
Moreover, downstream GA signaling components are affected by the sRNA regulatory network as well. StMYB33, an orthologue of Arabidopsis GA-induced MYB-like proteins (GAMYBs) MYB33, is regulated by miR319a, a close relative of miR159 [43]. The two Arabidopsis GAMYBs were shown to be regulated by miR159 [44]. The functional relation between lower activity of GA signaling was already directly confirmed to be related to disease severity in three different experimental systems.
Arabidopsis GAMYB double knockout showed ameliorated symptoms when infected with Cucumber mosaic virus (CMV) [38], similarly was shown in rice in interaction with bacteria (Xanthomonas oryzae pv. oyrzae) and fungi (Magnaporthe oryzae) using knockout in GA deactivating enzyme [45].
Also in line with this, increase in DELLA protein concentrations was shown to trigger components not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint of rhizobial and mycorrhizal signaling [46,47], showing yet another similarity between tolerant response of potato to viral infection and response of plants in mutualistic symbiosis. With the discovery of GAMYB binding sites in the MIR482f and MIR6022 promoters (S10 Dataset), we found the circuit in sRNA-GA signaling and additionally a link between GA signaling and R-gene expression. The complexity of regulatory responses observed in this study (Fig 5) is in line with the systems biology paradigm that interaction of multiple components and not a single component within a cell leads to much of biological function [48]. Although the reductionist approach is powerful in building logically simple hypotheses and devising ways to test them, it is very difficult to reconstitute the function for a whole biological system based solely on that as the behavior of the system may depend heavily on complex interactions within the system [49]. Thus, we have adopted a systems level confirmation of our hypothesis that sRNA-GA-immune signaling circuit is important for establishment of tolerant interaction. Previously, we had demonstrated that SA regulates plant responses to virus infection; not only by delaying viral multiplication, but also by controlling disease symptom severity, most probably via its effects on host primary metabolism [5]. In this study we have confirmed that response of sRNA regulatory networks controlling potential immune receptors and hormonal signaling is strongly attenuated in the NahG-transgenic plants in the early stage of viral infection (3 dpi ; Figs 1 and 2) showing the role of sRNAs in linking immune signaling and defense. The molecular mechanisms of the link between SA signaling and sRNA network are also complex. SA has been shown to induce RNAdependent RNA polymerase 1 expression, which is crucial for the maintenance of basal resistance to several RNA viruses [50,51] but none of the silencing mechanism related enzymes are regulated in SA deficient NahG-Désirée plants (S11 Dataset; [7]). We have here predicted SA-directed transcriptional regulation of MIR482f, the miRNA linking the GA signaling circuit and R-gene expression (Fig 5, S10 Dataset), which could be additional link between SA signaling and symptoms development in potato-PVY interaction. Regulation of miR482f expression by SA was indeed experimentally confirmed in comparison between non-transgenic and NahG-transgenic plants (S1 Dataset). An additional link are the WRKY transcription factors that are under positive control of SA ( Table 1, S11 and S12 Datasets) and were predicted to regulate promotors of all three miRNAs involved in sRNA -GA circuit (Fig 5).
The outcome of plant-pathogen interactions depends on the delicate balance between plant immune signaling network and its interaction with pathogen. Here, we focused on the roles of sRNA networks in establishment of the tolerance to PVY infection. We showed that the sRNA regulatory network links immune and developmental signaling in potato through newly discovered sRNA -GA circuit. Tolerance of potato to virus infection perturbs sRNA network resulting in downregulation of GA-mediated signaling, as well as modulation of R-gene transcript levels; this results in amelioration of disease symptoms. Supporting this, the responses observed for individual miRNAs, R-genes and GA-mediated not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint signaling bear a striking similarity to effects seen during mutualistic symbioses. Modulation of disease severity through GA-mediated signaling was also observed in Arabidopsis infected with mild and severe strains of CMV virus. It is thus plausible that a similar modulation of plant responses occurs in both mutualistic symbiosis and tolerance. This is in line with growing evidence showing that viruses, like other symbionts, lie on a continuum between antagonistic and mutualistic relationships [52,53]. not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Plant material
Potato leaves of cv. Désirée were mock or PVY NTN -inoculated (isolate NIB-NTN, GenBank accession no. AJ585342) as described previously [7]. Plant material of the inoculated leaves was collected 3 days post inoculation (dpi, corresponding to early stages of viral multiplication in both genotypes and before symptoms development in NahG-Désirée plants) for each treatment. Three and four biological replicates (individual leaves from different plants per group) were analyzed for RNA analysis and for hormonal measurements, respectively. Three plants from each group were monitored for plant height, till 21 dpi, when they all started to senesce. The same experimental set up was designed also for analysis of transgenic NahG-Désirée plants [54].

RNA extraction, library preparation and sRNA sequencing
Total RNA was extracted from 100 mg of homogenized leaf tissue using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA) and MaXtract High Density tubes (Qiagen, Hilden, Germany) following manufacturers' protocols. RNA concentration, quality and purity were assessed using agarose gel electrophoresis and NanoDrop ND-1000 spectrophotometer (Thermo Scientific). sRNA NGS libraries were generated from total RNA samples using the TailorMix

Identification of endogenous and virus-derived sRNAs
The raw reads were quality filtered using Filter Tool of the UEA sRNA Toolkit [55] by discarding low complexity reads (containing at most two distinct nucleotides), reads shorter than 18 nt and longer than 25 nt, reads matching tRNA/rRNA sequences and reads not mapped to the potato genome (PGSC_DM_v4.3) [56]. To identify known annotated miRNAs, the remaining reads were compared to all plant miRNAs registered in miRBase database (release 21) [57], allowing no mismatches. The sequences that matched mature miRNAs from other plants than potato (miRNA orthologs), were mapped to the potato genome to find corresponding MIR loci able to form hairpin structure [58] and named according to annotation of conserved miRNA [59]. miRNAs that had different 5′ and 3′ ends with respect to the mature miRNA, were annotated as isomiRs. To identify novel unannotated miRNAs, filtered reads were submitted to miRCat tool of the UEA sRNA Toolkit using default parameters for plants, considering only reads of lengths 18-24 nt. Reads were first mapped to the potato genome, then the 100 and 200 nt long windows around the aligned reads were extracted [58]. The predicted secondary structures were trimmed and analyzed to verify the characteristic hairpin pre-miRNA structure according to plant miRNAs annotation criteria [59]. An additional criterion we have imposed was that novel not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint miRNAs should be present at least in two analyzed samples with more than five raw reads. Potential novel miRNAs were mapped against miRBase and sequences that matched known plant miRNAs with up to three mismatches were excluded. The novelty of potato specific miRNAs was verified with the miRPlant version 5 [58] using default parameters and additionally rechecked against the latest releases of Rfam [60]; http://rfam.xfam.org/), tRNA [61]; http://gtrnadb.ucsc.edu/) and snoRNA databases [62]; http://snoopy.med.miyazaki-u.ac.jp/). Families of novel miRNAs were determined by clustering their sequences with sequences of known miRNAs using CD-HIT with identity threshold of 0.9 [63]. To identify vsiRNAs, reads of lengths 20-24 nt from all PVY NTN -infected samples were mapped to the reconstructed consensus PVY NTN genome [64] using CLC Genomics Workbench version 8 (http://www.clcbio.com/) allowing only 100 % identity.

Prediction of novel potato phasiRNAs and PHAS loci
Prediction of phasiRNAs and phasiRNA-producing loci (PHAS loci) was performed using ta-siRNA prediction tool [65] utilizing the potato genome [56] and the merged potato gene and unigene sequences StNIB_v1 [66]. To detect PHAS loci with a significant degree of phasing (p < 0.0001) and to investigate whether phasing also occurs in sRNAs with phase sizes different from 21 nt, different phasing intervals ranging from 21 to 24 nt were analyzed.

sRNA quantification and statistical analysis
Differential expression analysis was performed in R (R development Core Team, 2011; version 3.2.2), using the limma package [67]. In short, sRNA counts with a baseline expression level of at least two RPM (reads per million of mapped reads) in at least three samples were TMM-normalized (edgeR package, [68]) and log-transformed using voom function [69]. To identify differentially expressed sRNAs the empirical Bayes approach was used and the resultant p-values were adjusted using Benjamini and Hochberg's (FDR) method. Adjusted p-values below 0.05 were considered statistically significant.

Stem-loop RT-qPCR
Stem-loop RT-qPCR was used to quantify the expression of six target miRNAs in relation to the endogenous control (stu-mir167a-5p.1), which was determined to be the most robustly expressed in a sRNA sequencing dataset of potato plants that were uninfected or infected with a range of viruses (PVY, PLRV, PVS, PVX, PVT) (SRA accession no. SRP083247). TaqMan MicroRNA Assays (Thermo Fisher Scientific) were ordered according to the sRNA-Seq sequence of the selected miRNAs (S2 Table). Total RNA (1 µg) of the same samples as used for sRNA-Seq was DNase I (Qiagen) treated and reverse transcribed using SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific) and stemloop Megaplex primer pool (Thermo Fisher Scientific) following the manufacturer's protocol and previously optimized cycling parameters [70]. No template control RT reactions were set to assess not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprint (which was this version posted April 26, 2017. . https://doi.org/10.1101/130757 doi: bioRxiv preprint potential Megaplex primer pool background. qPCR reactions were performed in 10 µl volume on the LightCycler480 (Roche Diagnostics Ltd., Rotkreuz, Switzerland) in duplicates and two dilutions (8-and 80-fold) per sample using TaqMan Universal Master Mix II, no UNG (Thermo Fisher Scientific) and TaqMan MicroRNA Assays following the manufacturers' protocols. Additionally, for each miRNA assay, a standard curve was constructed from a serial dilution of the pool of all samples. Raw Cq values were calculated using the second derivative maximum method (Roche Diagnostics Ltd.) and miRNA expression was quantified using a relative standard curve method by normalization to the endogenous control [5]. The statistical significance was assessed by Student's t-test.

sRNA target prediction
In silico identification of potato transcripts targeted by identified sRNAs was carried out using the psRNATarget ( [71]; http://plantgrn.noble.org/psRNATarget/) and StNIB_v1 sequences [66], following previously described parameters [72]. http://sites.psu.edu/axtell/software/cleaveland4/) using all our experimentally identified sRNAs and the StNIB_v1 sequences allowing for maximum three mismatches. All identified degradation targets were classified into 5 categories as previously described [73]. Only categories 0-3 with high cleavage signal were considered as reliable cleavage. Results of miRNA-target (PHAS loci) interactions were also used to reveal miRNA triggers of the phasiRNA production. Only 22-nt miRNAs were kept as potential triggers [15,16].

Regulatory network construction
In order to compare the expression of sRNAs with the expression of their target transcripts we used a microarray gene expression dataset generated from the same samples ( [7]; GEO accession no. GSE58593). All differentially expressed miRNAs and phasiRNAs were analyzed for functional overrepresentation in biological pathways with MapMan software [74] using the ontology adapted for potato [66,75]. All sRNAs and their targets, obtained by in silico prediction and Degradome-Seq were integrated with their expression data and used for construction of regulatory networks in Cytoscape 3.4 [76].
not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Identification of cis-regulatory elements in promoter regions of MIR genes
1000 nt sequences upstream of the predicted MIR gene hairpin sequences were extracted as putative miRNA promoter regions [77] and scanned for cis-regulatory elements of plant transcription factors using position weight matrices and transcription binding sites in TRANSFAC [78]; http://www.biobaseinternational.com/product/transcription-factor-binding-sites) and PlantCARE [79]; http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) databases.

Hormonal measurements
Hormone contents were determined by gas chromatography coupled with mass spectrometry (GC-MS).
Approximately 100 mg of plant material was collected for each sample and 1 ml of 100% methanol (HPLC grade) was added. The following stable isotope-labeled compounds were added as internal and letting them rest for 30 min at room temperature. The setting for the GC and the MS were as described previously [80]. hormone concentrations between treatment-genotype groups were determined by ANOVA followed by not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.

Gene set enrichment analysis
To identify SA regulated genes in cv. Désirée the normalized expression values between mock NahG-Désirée vs. Désirée and PVY-infected NahG-Désirée vs. Désirée samples were compared (calculated from the 3 dpi samples; data of Stare et al. [7]. Gene Set Enrichment Analysis (GSEA; [81] was performed (false discovery rate corrected q ≤0.01) comparing expression profiles between both genotypes, using MapMan ontology as the source of the gene sets.

Data deposition and Gene IDs
The sRNA and Degradome-Seq data can be accessed at the NCBI's Gene expression omnibus (GEO) under accession numbers GSE84851 and GSE84967. Full list of gene/protein names used in this manuscript, together with their Gene IDs, short names, Arabidopsis orthologue genes is given in S1 Table.