Potato leaves of cv. Désirée were mock or PVY^NTN-inoculated (isolate NIB-NTN, GenBank accession no. AJ585342) as described previously (Stare et al., 2015). Plant material of the inoculated leaves was collected 3 days post inoculation (dpi), corresponding to early stages of viral multiplication. 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 (Halim et al., 2004).

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 miRNA Sample Preparation Kit (SeqMatic LLC, Fremont, CA) and sequenced on the Illumina HiSeq 2000 Sequencing System at SeqMatic LLC.

The raw reads were quality filtered using Filter Tool of the UEA sRNA Toolkit (Moxon et al., 2008) 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) (Xu et al., 2011). To identify known annotated miRNAs, the remaining reads were compared to all plant miRNAs registered in miRBase database (release 21) (Kozomara and Griffiths-Jones, 2014), 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 (An et al., 2014) and named according to annotation of conserved miRNA (Meyers et al., 2008). 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 (An et al., 2014). The predicted secondary structures were trimmed and analyzed to verify the characteristic hairpin pre-miRNA structure according to plant miRNAs annotation criteria (Meyers et al., 2008). An additional criterion we have imposed was, that novel 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 (An et al., 2014) using default parameters and additionally rechecked against the latest releases of Rfam (Nawrocki et al., 2015; http://rfam.xfam.org/), tRNA (Chan and Lowe, 2016; http://gtrnadb.ucsc.edu/) and snoRNA databases (Yoshihama et al., 2013; 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 an identity threshold of 0.9 (Huang et al., 2010). To identify vsiRNAs, reads of lengths 20–24 nt from all PVY^NTN-infected samples were mapped to the reconstructed consensus PVY^NTN genome (Kutnjak et al., 2015) using CLC Genomics Workbench version 8 (http://www.clcbio.com/) allowing only 100% identity.

Prediction of phasiRNAs and phasiRNA-producing loci (PHAS loci) was performed using ta-siRNA prediction tool (Chen et al., 2007; Moxon et al., 2008) utilizing the potato genome (Xu et al., 2011) and the merged potato gene and unigene sequences StNIB_v1 (Ramšak et al., 2014). Analysis of phasing was performed in 21- and 24-nt intervals. To detect PHAS loci with a significant degree of phasing, very strict criteria were applied to avoid detection of false positives (phasing p < 0.0001; number of unique phasiRNAs detected at specific PHAS locus ≥ 4, also to avoid detection of repeat-associated siRNAs). Additionally, the phasing p-values were corrected with strict Bonferroni correction and only loci with the p < 0.05 were considered as PHAS loci.

Differential expression analysis was performed in R (R Core Team, 2013; version 3.2.2), using the limma package (Ritchie et al., 2015). 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; Robinson et al., 2009) and analyzed using voom function (Law et al., 2014). 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 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 (Table S1). 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 and stem-loop Megaplex primer pool (both Thermo Fisher Scientific) following the manufacturer's protocol and previously optimized cycling parameters (Varkonyi-Gasic et al., 2007). Three different negative controls were included: no template control RT reactions to assess potential Megaplex primer pool background, RT-minus controls to check the presence of the signal that could be the result of contaminating DNA and no template qPCR control reactions for excluding the contamination of the PCR reagents. All controls were negative. 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 using quantGenius (Baebler et al., 2017) The statistical significance was assessed by Student's t-test.

In silico identification of potato transcripts targeted by identified sRNAs was carried out using the psRNATarget (Dai and Zhao, 2011; http://plantgrn.noble.org/psRNATarget/) and StNIB_v1 sequences (Ramšak et al., 2014), following previously proposed stringent parameters (Zhang et al., 2013). Moreover, targets of identified sRNAs were experimentally validated with parallel analysis of RNA ends (PARE) Degradome-Seq. The four degradome libraries (mock Désirée, PVY Désirée, mock NahG-Désirée, PVY NahG-Désirée) were constructed by pooling RNA of the biological replicates and sequenced on the Illumina HiSeq 2500 platform. The data were analyzed at LC Sciences (Houston, TX) with CleaveLand4 (Addo-Quaye et al., 2009; 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 five categories as previously described (Addo-Quaye et al., 2009). Category “0” is defined as >1 raw read at the position, with abundance at a position equal to the maximum on the transcript, and with only one maximum on the transcript. Category “I” is described as >1 raw read at the position, with abundance at the position equal to the maximum on the transcript, and more than one maximum position on the transcript. Category “II” includes >1 raw read at the position and abundance at the position less than the maximum but higher than the median for the transcript. Category “III” comprised the transcripts with >1 raw read at the position, and abundance at the position equal to or less than the median for the transcript. Category “IV” comprised transcripts with one raw read at the cleavage position. Only categories with high confidence of cleavage (0, I, II, III) were considered for biological interpretation. 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 (Chen et al., 2010; Cuperus et al., 2010).

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 (Stare et al., 2015; GEO accession no. GSE58593). All differentially expressed miRNAs and phasiRNAs were analyzed for functional overrepresentation in biological pathways with MapMan software (Usadel et al., 2009) using the ontology adapted for potato (Rotter et al., 2007; Ramšak et al., 2014). All sRNAs and their targets, obtained by in silico prediction and Degradome-Seq were integrated with their expression data and used for the construction of regulatory networks in Cytoscape 3.4 (Shannon et al., 2003).

1000 nt sequences upstream of the predicted MIR gene hairpin sequences were extracted as putative miRNA promoter regions (Megraw and Hatzigeorgiou, 2010) and analyzed for cis-regulatory elements of plant transcription factors using position weight matrices and transcription binding sites in TRANSFAC (Matys et al., 2003; http://www.biobase-international.com/product/transcription-factor-binding-sites) and PlantCARE (Lescot et al., 2002; http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) implemented algorithms.

Hormone contents were determined by gas chromatography coupled with mass spectrometry (GC-MS) from ~100 mg of plant material. One milliliter of 100% methanol (HPLC grade) and a mix of 50 pmol stable isotope-labeled internal standards were added to each sample prior to extraction. The samples were heated (60°C, 5 min) and then incubated (room temperature, 1 h) with occasional vortexing. The methanolic phase was taken to complete dryness in vacuo. The resulting residue was dissolved in methanol (50 μl) to which diethyl ether (200 μl) was added. The samples were sonified (5 min) and centrifuged (5 min, 14,000 g). The particle-free supernatant was loaded to aminopropyl solid-phase extraction cartridges (Chromabond NH₂ shorty 10 mg; Macherey–Nagel GmbH, Düren, Germany). Each cartridge was washed twice with CHCl₃:2-propanol (2:1, v/v, 250 μl) before the hormone-containing fraction was eluted with acidified diethyl ether (2% acetic acid, v/v, 400 μl). The eluates were transferred into 0.8 ml autosampler vials and again taken to dryness in a gentle stream of nitrogen. Prior to MS assessment, the samples were derivatized with a 20 μl of a mix of l of acetone:methanol (9:1, v/v, 220 μ), diethyl ether (27 μl) and (trimethylsilyl)diazomethane solution (2.0 M in diethyl ether, 3 μl) and letting them rest for 30 min at room temperature. The setting for the GC and the MS were as described previously (Sanz et al., 2014). For the determination of endogenous and stable isotope-labeled methylated acidic plant hormones, respectively, the following ion transitions were recorded: MeSA m/z 152 to m/z 120 and m/z 156 to m/z 124 for [²H₄]-MeSA, retention time 6.75 ± 0.4 min; MeOPDA m/z 238 to m/z 163 and m/z 243 to m/z 168 for [²H₅]-MeOPDA, retention time 10.00 ± 0.4 min; MeJA m/z 224 to m/z 151 and m/z 229 to m/z 154 for [²H₅]-MeJA, retention time 11.27 ± 0.5 min; MeIAA m/z 189 to m/z 130 and m/z 191 to m/z 132 for [²H₂]-MeIAA, retention time 13.34 ± 0.4 min; MeABA m/z 162 to m/z 133 and m/z 168 to m/z 139 for [²H₆]-MeABA, retention time 15.78 ± 0.4 min; and MeGA m/z 136 to m/z 120 and m/z 138 to m/z 122 for [²H₂]-MeGA, retention time 21.67 ± 0.6 min. The amounts of endogenous hormone contents were calculated from the signal ratio of the unlabeled over the stable isotope-containing mass fragment observed in the parallel measurements. Significant changes in hormone concentrations between treatment-genotype groups were determined by ANOVA followed by LSD post hoc analysis (Benjamini Hochberg FDR p-value adjustment, alpha = 0.05) using the Agricolae R package.

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. (2015). Gene Set Enrichment Analysis algorithm (GSEA; Subramanian et al., 2005) was run to perform analysis [false discovery rate (FDR) corrected q ≤ 0.01] of expression profiles between both genotypes, using MapMan ontology as the source of the gene sets.

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

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 (Figure S1, Dataset S1). 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 (Datasets S1, S2).

When assessing the miRNA regulatory network, the amplification of silencing through phasiRNA biogenesis was also considered. In total, more than 400 PHAS loci were predicted, coding for 1,513 phasiRNAs. Two hundred and forty-eight PHAS 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 repeat-receptor-like kinases) (Datasets S3, S4). Moreover, we also searched for miRNAs with the ability to trigger phasiRNA production and similarly to previous reports (Li et al., 2011; Shivaprasad et al., 2012), we observed that the vast majority of predicted miRNA triggers belong to miR482, miR6023, miR6024 and miR6027 families (Dataset S5).

We further compared miRNA expression profiles of PVY^NTN-infected vs. mock-inoculated leaves in early stages of virus infection (3 dpi, 1 day before detectable viral multiplication; Baebler et al., 2011; Stare et al., 2015). In total, 61 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 (Dataset S1). Additionally, we identified 36 phasiRNAs as differentially regulated 3 dpi, mainly originating from non-coding PHAS and NBS-LRR loci (Dataset S4). To validate the obtained sRNA-Seq results, abundance of six differentially expressed miRNAs was analyzed by stem-loop RT-qPCR. As shown in Figure S2, all sRNA-Seq differential expression results were confirmed, except in cases were concentration of miRNAs were below the limit of the quantification.

In previous studies, a plethora of potato miRNA/phasiRNAs has 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 (Dataset S6). Additionally, the predictions of target cleavage were experimentally validated by Degradome-Seq (Dataset S7). Based on this information we constructed a potato sRNA regulatory network connecting miRNAs with phasiRNAs and their targets (Supplementary Online Files 1, 2 are available online at http://projects.nib.si/podefsig/datasets). This revealed several already known and many novel connections linking sRNA regulation to the plant immune signaling (see example in Figure S3, Datasets S6, S7). Several miRNA-mRNA pairs conserved across plant species, such as miR156-SPL11, miR160-ARF10, miR172-AP2, or miR396-GRF5 (Curaba et al., 2014), were confirmed also in our system (Datasets S6, S7). Our data also showed the miR393-mediated cleavage of transcripts encoding members of TIR/AFB gene family, receptors implicated in the control of auxin signaling (Figure 1) (Si-Ammour et al., 2011). We also discovered that these transcripts were targets of several TIR1-derived phasiRNAs (phasiTIRs) (Figure 1, Datasets S6, S7). Moreover, our in silico analysis predicted that the miR393- and phasiTIR-network is targeting downstream transcription factor StARF1 and other phytohormone signaling pathways, such as transcripts involved in ethylene signaling, in jasmonate signaling and in brassinosteroid biosynthesis (Figure 1, Datasets S6, S7). Two of these alternative phasiTIR targets, StAP2 and StOPR1, were also confirmed by degradome sequencing. Interestingly, our analysis also showed that the majority of differentially expressed miRNAs and phasiRNAs target genes are coding for defense-related proteins such as pathogenesis-related (PR) proteins and proteins involved in the biosynthesis of secondary metabolites, transcription factors belonging to AP2/ERF, bHLH, MYB, and GRAS family proteins, putative immune receptors (NBS-LRRs, LRR-RLKs) as well as proteins involved in biosynthesis and signaling of different phytohormones (Figure 2A).

Analysis of the differentially expressed miRNAs and phasiRNAs together with the levels of the target transcripts (data published in Stare et al., 2015) 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 (Datasets S6, S7). Moreover, NBS-LRRs are regulated also by phasiRNAs, where most phasiRNAs have multiple NBS-LRR targets (Figure S4, Datasets S6, S7) due to the shared conserved P-loop or Walker A motif (Shivaprasad et al., 2012). In all of the previous studies, miR482 family members were downregulated following pathogen infection (Shivaprasad et al., 2012; Ouyang et al., 2014; Yang et al., 2015). In our study, however, the one regulated member of the miR482 family (miR482e) targeting NBS-LRR transcripts was upregulated following PVY^NTN infection (Dataset S1), similarly as observed in the establishment of mutualistic symbiosis in soybean roots (Li et al., 2010). Moreover, several miRNAs that were upregulated in response to PVY^NTN in cv. Désirée, such as miR164, miR167, miR169, miR171, miR319, miR390, and miR393 have also been reported to regulate nodulation and arbuscular mycorrhizal symbiosis in different legume species (Dataset S1; Lelandais-Brière et al., 2009; Devers et al., 2011; Mao et al., 2013; Yan et al., 2015). In addition to NBS-LRR proteins, LRR-RLKs are also important mediators of immune as well as important triggers of mutualistic symbiosis signaling cascades (Antolin-Llovera et al., 2014; Hohmann et al., 2017). We have predicted a novel miRNA-LRR-RLKs 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 (Figure S4A, Dataset S6).

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 (Dataset S8). 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 fulfill the criteria for a viral miRNA. Subsequently, we searched for potential host transcripts targeted by vsiRNAs in our experimentally validated target degradation dataset (Dataset S9). We found that vsiRNAs are indeed able to target multiple potato transcripts, among them mRNAs coding for immune receptor proteins, various transcription factors and proteins involved in hormonal signaling pathways (Dataset S9). For example, several vsiRNAs were detected with the confirmed ability to guide the cleavage of transcripts encoding Heat shock proteins and RLKs implicated in stress signaling responses, as well as GRAS transcription factors, ARFs and Auxin induced-like proteins involved in regulation of growth and developmental processes (Figure S5).

Interestingly, we found that GA biosynthesis and downstream signaling are targeted by a sRNA-mediated regulatory network and that the changes in sRNA levels following PVY^NTN infection are reflected also in the changes of their target transcripts levels (Figure 2B, Table 1, Dataset S6). GA20-oxidase (GA20ox) and GA3-oxidase (GA3ox) are enzymes that catalyze the last steps in the formation of bioactive GAs (Yamaguchi, 2008). In Désirée plants miR167e was predicted to cause translational repression of the StGA20ox transcript (Figure 2B, Dataset S6). An additional layer of GA biosynthesis regulation is represented by the increased production of phasiRNA931, which promotes cleavage of the StGA3ox transcript (Figure 2B, Table 1). The transcriptomics results support these interactions as the targeted transcripts are significantly downregulated in Désirée upon PVY^NTN infection (Figure 2B, Table 1, Dataset S6). Additionally, vsiRNAs were found to target transcripts encoding two enzymes involved in GA biosynthesis StGA1 and StGA20ox (Figure 2B, Dataset S9). 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 miR167e was also discovered in tomato (Griffiths-Jones et al., 2006; Zhang et al., 2014). sRNAs can target downstream GA signaling in the potato-PVY interaction on multiple levels. Four miR319 family members, close relatives of miR159 family (Palatnik et al., 2007), were predicted to cleave the transcript encoding StMYB33, an ortholog of GAMYB transcription factor involved in GA signal transduction (Millar and Gubler, 2005). Moreover, phasiRNA931 can cleave a potato ortholog of a DELLA protein, a GA-signaling repressor (Figure 2B, Dataset S6).

Such interconnectedness between plant defense-related miRNA/phasiRNA network and GA biosynthesis/signaling has not been previously identified and may represent a link between defensive and developmental signaling. Thus, we decided to functionally evaluate 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 GA₃ in PVY^NTN infected Désirée plants (Figure 2C). 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 (Figure 2C, Dataset S10). To inspect if GA deficiency has any impact on plant growth, plant height was monitored till 21 dpi. No differences were observed between all four studied groups of plants.

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, MIR6022, MIR319a, and MIR167e were investigated together with MIR482f and MIR482c genes encoding miRNA triggers of phasiRNA931 production. Interestingly, GAMYB binding sites were detected in the promoters of the MIR6022 and MIR482f genes (Dataset S11). Moreover, these genes harbor WRKY8/28/48 binding sites, while MIR319a and MIR482c harbor a general WRKY binding W-box regulatory element. Additionally, cis-acting elements involved in SA and JA responsiveness were identified in the promoter of the MIR482f gene (Dataset S11). In all these four analyzed promoter sequences NAC transcription factor binding sites were detected. 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 (Dataset S11).

As the link between repression of GA signaling and disease symptoms severity was already established in Arabidopsis (Du et al., 2014) we have further investigated the activity of discovered sRNA-GA circuit in susceptible potato-PVY^NTN interaction. We have shown previously that SA depletion breaks the equilibrium between disease and tolerance. In NahG-Désirée plants, the pronounced disease symptoms appeared both on the inoculated and systemic leaves (Baebler et al., 2011). 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 (Baebler et al., 2011). Here, we performed the analysis of sRNA response as well as the measurements of hormonal concentrations in the 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 28 miRNAs were differentially expressed, with the majority showing a lower degree of induction than in Désirée plants (Figures 3, 4A, Dataset S1).

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 were diminished in NahG-Désirée plants (Figure 4B, Table 2, Dataset S1). 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 (Figure 4B, Table 2, Dataset S6). 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 (Figure 4C). Furthermore, the relieved silencing of LRR-RLKs by miR6022 that is predicted to modulate the immune response and which is also linked with GAMYBs was absent in NahG-Désirée plants (Figure 4B, Table 2, Figure S4B). We also inspected sRNA-mediated responses which resembled responses in mutualistic symbioses in tolerant plants of Désirée and found that miR482e was also upregulated in NahG plants, while regulation of miR164, miR167, miR169, miR171, miR319, miR390, and miR393, was diminished (Dataset S1).

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 SA-depleted Désirée (Figure 3). We found 36 miRNAs regulated by SA. It seems that SA in the normal growing conditions generally cause the reduction in the levels of miRNAs as the majority (27) of miRNAs were detected at significantly higher levels in NahG-Désirée plants (Dataset S1). When we similarly compared the transcript expression profiles of the two genotypes we noticed that the most strongly induced by SA signaling are notably different WRKY transcription factors (Table 3, Datasets S12, S13), among them orthologs of Arabidopsis WRKY70, which was already shown to be positively regulated by SA (Li et al., 2004). As MIR319a, MIR482c, and MIR482f promoters harbor WRKY transcription factors binding sites (Dataset S11), we discovered a potential direct link between SA signaling and miRNA regulatory network in potato. This link was experimentally confirmed by differential response of these three miRNAs to PVY infection in non-transgenic vs. SA-depleted genotype. miR319a was only induced following PVY infection in non-transgenic plants, supporting the hypothesis that it requires WRKY transcription factor for this response. On the other hand, miR482c and miR482f accumulation only differs in SA-depleted plants indicating more complex regulation indicated also by prediction of additional SA responsive elements in their promoter regions (Figure 5, Datasets S1, S11).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.