During this study, S. lycopersicum plants were maintained in a greenhouse at 23°C during the day and at 21°C at night (16 h light/8 h dark regime) and relative humidity of 60%. S. lycopersicum cv. Moneymaker (MM) was used as susceptible control and a Ty-1 introgression line was derived from S. chilense LA1969 (Verlaan et al., 2013). Plants were infected with TYLCV via agroinoculation with an infectious clone of the TYLCV Israel strain isolated from Almeria, Spain, as described by Morilla et al. (2005) (GenBank AJ489258.1) or with untransformed agrobacteria as control (strain LBA4404).

The protocol as described in Voorburg et al. (2020) was followed to perform the agrobacterium transient transformation assays, which is a slightly modified protocol of Bucher et al. (2003). In brief, A. tumefaciens was grown overnight at 28°C in 3 ml LB3 (10 gl⁻¹ trypton, 5 gl⁻¹ yeast, 4 gl⁻¹ NaCl, 1 gl⁻¹ KCl, 3 gl⁻¹ MgSO₄·2H₂O) medium containing proper antibiotic selection pressure. From this culture, 600 μl was incubated overnight in 3 ml induction medium (10.5 gl⁻¹ K₂HPO₄, 4.5 gl⁻¹ KH₂PO₄, 1 gl⁻¹ (NH₄)₂SO₄, 0.5 gl⁻¹ sodium citrate. 2H₂O, 1 mM MgSO₄·7H₂O, 0.2% (w/v) glucose, 0.5% (v/v) glycerol, 50 μM acetosyringone, 10 mM 2-(N-Morpholino) ethanesulfonic acid [(MES), pH 5.6)]. The next day bacteria were pelleted by centrifugation (15 min. at 2,670 g) and resuspended in MS MES buffer [Murashige and Skoog medium (Duchefa biochemie) supplemented with 150 μM acetosyringone, 10 mM MES and 87 mM sucrose] at an OD600 of 0.5. The abaxial side of the first two true leaves of the 3-week old tomato seedlings was infiltrated with the agrobacteria by pressure inoculation with a needle-less syringe.

First, systemically infected leaves (top leaves) were snap-frozen in liquid nitrogen and stored at −80°C. RNA was isolated from this tissue using the mirVana miRNA isolation kit (Life Technologies) according to the protocol of manufacturers. In brief, upon grinding the leaf tissue in liquid nitrogen, the cells were lysed by adding lysis/binding buffer in a 1:10 ratio (w/v) and mixing thoroughly. Subsequently, 1/10 volume of miRNA homogenate additive was added, followed by a 10-min incubation on ice. RNA was extracted by adding acid-phenol: chloroform, mixing and spinning the sample. The aqueous phase was transferred to a new tube and 1/3 volume of 100% ethanol was added. Upon mixing thoroughly, the sample was passed through a Filter Cartridge. From this point, both the filtrate, containing the sRNAs, as the filter cartridge, containing the RNA fraction depleted of sRNAs, was further processed. To the filtrate 2/3 volume, 100% ethanol was added and this mixture was passed through a second filter cartridge. This filter cartridge was washed once with miRNA wash solution 1 and subsequently washed two times with wash solution 2/3. The sRNA fraction was eluted with MiliQ. The first filter cartridge containing the RNA fraction depleted of sRNAs was washed once with miRNA wash solution 1 and subsequently washed two times with wash solution 2/3. This RNA fraction was also eluted with MiliQ. The integrity of the RNA fractions was checked on an RNAse free agarose gel. From the samples used for the RNA isolations, also DNA was isolated following the cetyltrimethyl ammonium bromide (CTAB) method of Doyle and Doyle (1987) with slight modifications as described by Doyle and Doyle (1987), Fulton et al. (1995), and Voorburg et al. (2020). In brief, grinded plant material was mixed with CTAB buffer (0.1 M TRIS, 0.7 M NaCl, 0.01 M EDTA, 2% CTAB) and subsequently incubated for 1 h at 65°C. After chloroform isoamyl alcohol extraction, the DNA was precipitated from the aqueous phase by adding isopropanol in a 1:1 ratio. DNA was pelleted, dried, and dissolved in MilliQ. DNA and RNA concentrations were measured with a Nanodrop ND-1000 device.

The sRNA fraction was analyzed by capillary electrophoresis on a Shimadzu MultiNA microchip. The QIAseq miRNA library kit (Qiagen) was used according to the manufacturer's instructions to synthesize cDNA. In short, oligonucleotide adapters were ligated to the 5′ and 3′ ends of the RNA samples. First-strand cDNA was synthesized using M-MLV reverse transcriptase and the 3′ adapter as a primer. Amplification of the resulting cDNA was performed by PCR with a high-fidelity DNA polymerase in 16 PCR cycles. The 3′ sequencing adapters included a sample specific six nucleotide long barcode sequence. The cDNA was purified using the magnetic beads provided in the kit. The cDNA was analyzed by capillary electrophoresis on a Shimadzu MultiNA microchip. The cDNA samples were pooled in approximately equimolar amounts and the size range of 160–200 bp was selected from a polyacrylamide gel. An aliquot of the size fractionated pool was analyzed by capillary electrophoresis on a Shimadzu MultiNA microchip. The cDNA pool was sequenced on an Illumina NextSeq 500 system using a 75 bp read length. The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/sra/PRJNA764425/.

The data analysis was performed on the Galaxy server (Afgan et al., 2018). Quality assessments were performed with the FastQC tool (Galaxy version 0.72+galaxy1). From each read, the adapter sequences were removed using the clip tool (Galaxy Version 1.0.3+galaxy0). Sequence length was determined by using the compute sequence length tool (Galaxy Version 1.0.3) and the count tool (Galaxy Version 1.0.3). Reads were mapped to the tomato genome (solgenomics, SL4.0) using Bowtie2 (Galaxy Version 2.3.4.3+galaxy0, default settings), and read length of tomato aligned reads calculated with the compute sequence length tool (Galaxy Version 1.0.3) and the count tool (Galaxy Version 1.0.3). Reads were mapped to the viral genome (TYLCV, GenBank AJ489258.1) using Bowtie2 (Galaxy Version 2.3.4.3+galaxy0, allowing 1 mismatch). From the bowtie output, genome profiles were generated by using the count (Galaxy Version 1.0.3), filter (Galaxy Version 1.1.1), join two datasets (Galaxy Version 2.1.3), compute (Galaxy Version 1.3.1), and advanced cut (Galaxy Version 1.1.0) tools. To determine the 5′ end nucleotides, the sequence logo tool (Galaxy Version 3.5.0) was used.

Relative virus titers of TYLCV were determined by qPCR using actin (Solyc04g011500) as the internal control. The reaction mixture contained 1× SybrSelect (Applied Biosystems), 300 nM forward primer, 300 nM reverse primer, and 10 ng genomic DNA. The primers TYLCV-F (5′-TTCGTCTAGATATTCCCTATATGAGGAGGTA-3′) and TYLCV-R (5′-GGGAAGCCCATTCAAATTAAAGG-3′) were used to amplify TYLCV, and actin-F (5′-GAAATAGCATAAGATGGCAGACG-3′) and actin-R (5′-ATACCCACCATCACACCAGTAT-3′) to amplify actin (Powell et al., 2012; Verlaan et al., 2013). The qPCR was performed in an Applied Biosystems AB7500, using the following cycling conditions: 2 min 95°C, 40 cycles of 15 s 95°C and 1 min 60°C, followed by a melting curve with 0.5°C steps from 60 to 95°C to determine PCR specificity. Relative viral titers were calculated using the DeltaDeltaCt method (Livak and Schmittgen, 2001). Values were normalized relative to the internal control and calibrated to levels of the TYLCV infected MM plants, which were set as 1.

Tomato lines containing Ty-1 were previously shown to contain increased amounts of vsiRNAs during TYLCV infection compared to susceptible MM plants. However, nothing was known yet on the exact role of Ty-1 in the biogenesis of vsiRNAs, i.e., on the size distribution of vsiRNAs and their genomic distribution profile. To determine the effect of the resistance gene Ty-1 on vsiRNA profiles in more detail, sRNAs were purified from infected and uninfected susceptible and Ty-1 encoding tomato plants and then analyzed via deep sequencing. The used Ty-1 bearing tomato plants are an introgression line, and it was previously shown that silencing of Ty-1 breaks the resistance, ruling out the involvement of another resistance gene (Verlaan et al., 2013). First, five 3-week old Ty-1 bearing tomato seedlings and susceptible tomato MM seedlings were either infected with TYLCV by agroinoculation or, as the negative control, infiltrated with untransformed agrobacteria. As expected, during the progression of the disease MM plants infected with TYLCV showed clear yellowing and curling of the top leaves (Figure 1A), while infected Ty-1 plants did not show any symptoms, just like uninfected control plants. At 29 days postinfection (dpi) systemic leaves were harvested. From each group, three plants were used to collect three independent biological replicates. To confirm a successful infection, DNA was extracted and virus titers determined via qPCR. The fold difference in titer was calculated using the DeltaDeltaCt method (Livak and Schmittgen, 2001) (Figure 1B). While viral titers in infected MM plants were high, those in the Ty-1 introgression line were reduced ~20-fold. No virus was detected in the mock control plants, which all were in accordance with previous results (Voorburg et al., 2020).

After having verified viral titers in the mock treated and virus-infected MM and Ty-1 tomato introgression line, systemic leaves collected from each treatment were used for the purification of sRNAs and subsequent library preparation. Deep sequencing of the sRNAs from three biological replicates of each treatment resulted in libraries containing in total between 9.5 and 11.8 million clean reads, varying in length from 15–70 nts (Table 1). The libraries of mock infected samples consisted almost fully of reads mapping to the tomato genome (average 97.7%) and with 0% of the reads mapping to the viral genome (Figure 2A, Table 1). siRNAs mapping to the viral genome (vsiRNAs) in samples of TYLCV infected plants represented on average 4.0% and 3.1% of total reads in MM and Ty-1 plants, respectively (Figure 2A). The amount of viral reads was comparable between the MM and Ty-1 plants, while previously an increased amount of vsiRNAs in Ty-1 plants compared to MM was reported (Butterbach et al., 2014). However, it is important to note the lower viral titer determined in Ty-1 plants. When the ratio of viral reads relative to the viral titer was determined, a strong increase in vsiRNAs per viral genome was observed in Ty-1 plants (Table 1, Figure 2B).

To further characterize the sRNAs, the size distribution was determined within the total pool of sRNAs. The most abundant sRNAs, as observed in samples from all treatments, were from the 24-nt size class, covering about 42–54% of all sRNA reads (Figure 3A). The sRNA reads were next split into those sRNAs mapping to the tomato genome and those mapping to the TYLCV genome (vsiRNAs). Also, among the tomato specific sRNAs, the 24-nt size class was by far the most prevalent, representing 54% of all sRNAs in uninfected MM, 53% in uninfected Ty-1 plants, 48% in TYLCV-infected MM plants, and 46% in TYLCV infected Ty-1 plants (Figure 3B). Besides the 24-nt sRNAs, the other most frequently found sRNAs were from 23-, 22-, and 21-nt size classes, which were all present in relative comparable amounts (8–16%) (Figure 3B). vsiRNAs were only collected from the TYLCV infected leaf samples and not found, as expected, in the mock infected plants. The most predominant class of vsiRNAs collected from susceptible MM, were from 21-nt (37%), followed by 22-nt (30%) and the less abundant classes of 24-nt (8%) and 23-nt (5%) (Figure 3B). In Ty-1 plants, this pattern was changed, with the largest size class presented by 22-nt (38%), followed by 21-nt (18%), 24-nt (14%), and 23-nt (9%) (Figure 3B). So, the class of the 21-nt vsiRNAs was reduced by ~50% in Ty-1 tomato compared with the susceptible MM, whereas the 22-, 23-, and 24-nt size classes were clearly increased (~25% increase for 22-nt and ~75% increase for 23- and 24-nt vsiRNAs).

Having determined the size distribution of the vsiRNAs, the polarity of the reads was analyzed. This was done for the total viral reads, and the 21-, 22-, and 24-nt size classes. In MM plants, the majority of all viral reads (53%) and those of the size class 21-nt (54%) and 22-nt (54%) were in sense polarity (clockwise orientation of the viral genome), whereas of the 24-nt size class 46% was in sense polarity (Figure 4A). In Ty-1 tomato, the amount of sense reads showed a small reduction of 3–6% in all the studied size classes when compared to MM plants, with 47% of all viral reads, 48% of 21-nt reads 51% of 22-nt reads, and 40% of 24-nt reads being in sense polarity (Figure 4A).

In light of the importance of the 5′ terminal nucleotide of siRNAs toward loading into different AGO proteins (Kim, 2008; Mi et al., 2008), the 5′ terminal nucleotide of viral reads was analyzed. In both MM and Ty-1 tomato reads from the 21-nt and 22-nt classes were enriched in a 5′ terminal uridine (U), with a prevalence of around 40% in both groups. The second most dominant 5′ terminal nucleotide was adenosine (A), with an average prevalence of around 27% in both groups (Figure 4B). Within the 24-nt reads, A was observed most prevalent (~40%), followed by U (~30%) (Figure 4B). The same patterns were observed when reads of viral sense and antisense polarity were separately analyzed (Supplementary Figure S1), with the exception of two cases. Within the 21-nt and 22-nt reads of antisense polarity U and A were more equally presented at the 5′ end, whereas U was dominant within the sense reads. The 24-nt reads of sense polarity showed a comparable proportion of U and A as 5′ terminal nucleotide, while A was more dominant in the reads of antisense polarity (Supplementary Figure S1). No clear differences were observed in 5′ terminal nucleotides between vsiRNAs isolated from MM and Ty-1 plants.

To examine the genomic distribution of vsiRNAs produced and analyze for possible changes in the presence of Ty-1, the reads derived from infected susceptible MM and resistant Ty-1 tomato were plotted against the TYLCV genome. The results showed that in both MM and Ty-1 plants the vsiRNAs cover the entire genome in viral sense and antisense orientation (Figure 5). For both MM and Ty-1 plants, the genomic distribution was somewhat heterogenous, with several regions being covered by more and others by less vsiRNAs. This was reflected in all vsiRNA size classes (Figure 5) and observed highly consistent in all three biological replicates (Supplementary Figure S2). Hotspots, regions with a higher density of vsiRNAs per nucleotide, were visible in the viral genome distribution of both MM and Ty-1 tomato, but these hotspots appeared more distinct in MM. Hotspots coincided with regions of overlapping ORFs, i.e., V1/V2, C1/C2/C3, and C1/C4 respectively (Figure 5, Table 2A), and were present in both the viral sense and antisense orientation and observed for all vsiRNA size classes. Only in the profile of 24-nt vsiRNAs, the hotspot at the C1/C2/C3 overlapping region was less pronounced.

In both MM and Ty-1 tomato, the intergenic region (IR) had the lowest abundance of vsiRNAs (Table 2A). Other regions with a low abundance in MM plants included the V2 start (5′ end of the V2 ORF until the start of V1 ORF), the V1–C3 junction (where the 3′ ends of the V1 and C3 transcripts are assumed to overlap), and the middle part of C1 (between the end of C4 and the start of C2) (Figure 5, Table 2A).

Not only the distribution of vsiRNA in all size classes was changed in the presence of Ty-1, the amount of reads per size class in MM and Ty-1 plants was also changed (Figure 3B). This complicated the comparison of the genome distribution between the two groups. Therefore, per size class, the percentage of mapped vsiRNAs per nucleotide was plotted (Supplementary Figure S3). The three hotspots earlier observed at the regions with overlapping ORFs were again visible. However, when the proportion of vsiRNAs within the size class 21, 22, and 24-nts mapping to these hotspots was calculated, a higher amount was obtained for MM (55%) compared with Ty-1 tomato (41%), indicating that in Ty-1 tomato the hotspots were less targeted. To further quantify the observed differences in vsiRNA genome distribution between MM and Ty-1 tomato, the number of reads/nucleotide was calculated for several genomic regions and compared with the reads/nucleotide in the situation of a completely homogenous distribution (Table 2A). The results strengthened the observation of less pronounced hotspots at the ORF overlapping regions of TYLCV in Ty-1 tomato compared with the situation in MM. Interestingly, the region at the 3′-end of the V1 and C3 ORF (V1–C3 junction), which has a low abundance of vsiRNA in MM, showed to have a higher abundance than average in Ty-1 plants (Table 2A, region “V1–C3 junction”). The increased density at V1 was in agreement with an earlier study, where the ratio of vsiRNAs mapping to V1/C3 changed in the presence of Ty-1 (Butterbach et al., 2014). When the ratio was calculated for the same regions as used in the previous study, a comparable shift was observed (Table 2B).

To better visualize the change in genome distribution and obtain a broader overview image, the amount of reads per 100 nucleotides were plotted (Figure 6). Also, the percentage of mapped vsiRNAs per 100 nucleotides was plotted to compare the two groups without a possible bias by a difference in absolute reads (Supplementary Figure S4). Those plots confirmed the previous finding of the high density of vsiRNAs at the regions with the overlapping ORFs V1/V2 and C1/C4 in both MM and Ty-1 plants. It also revealed that in Ty-1 plants less vsiRNAs accumulated at the region with the overlapping ORFs C1/C2/C3, while this region is a clear hotspot in MM (Supplementary Figure S4; Table 2A). In addition, the amount of vsiRNAs that were mapped to the IR and V1 was increased in Ty-1 plants, especially in an antisense orientation (Supplementary Figure S4; Table 2A).

Comparison of the genomic distribution profiles of the vsiRNAs collected from susceptible MM and resistant Ty-1 tomato indicated that the IR and the region “V1–C3 junction” were relatively more targeted in Ty-1 tomato. The percentage of vsiRNAs to those regions was increased 3 and 1.5 fold, respectively (Table 2A). A previous study reported that in MM plants the IR was targeted relatively more by vsiRNAs of 24-nt than by those of 21-nt and 22-nt (Piedra-Aguilera et al., 2019). A closer look at the data from this study confirmed these observations (Table 2C). This study in addition showed that in Ty-1 tomato, the amount of vsiRNAs mapping to the IR increased for all three size classes (21-, 22-, and 24-nt) in both polarities, but with the highest amounts of antisense vsiRNAs (Table 2C). In both MM and Ty-1 tomato, the antisense reads of 24-nts vsiRNAs were most prevalent (3.6 and 6.1%, respectively, Table 2C).

To analyze whether the increased amount of vsiRNAs from the IR and “V1–C3 junction” region mapped to specific sequences (e.g., secondary structures) or were distributed comparable to the situation in MM, genomic distribution graphs were made zoomed in on these regions (Figure 7). In the region “V1–C3 junction” the increased amounts of vsiRNAs from Ty-1 tomato distributed all over the region (Figure 7B). In contrast, vsiRNAs from the IR appeared to concentrate in two “sub-regions,” namely around nucleotides 2624–2665 and nucleotides 41–84 of the TYLCV genome (Figure 7A). Earlier, reports of resistance breaking strains in Ty-1 tomato described recombination events between TYLCV and Tomato yellow leaf curl Sardinia virus (TYLCSV) taking place in the IR region (Belabess et al., 2015, 2016, 2018; Torre et al., 2018; Urbino et al., 2020), and with the smallest recombination event reported for the resistance breaking TYLCV-IS76 (Urbino et al., 2020). To analyze whether this region coincided with the hotspot “sub-regions” of the IR found in this study, the region of recombination in TYLCV-IS76 was compared with the sub-regions of increased vsiRNAs in Ty-1 plants. The results showed that one region (nucleotide 41–84) indeed overlapped with the recombination event (Figure 7A).

Recently, an alternative Dicer-independent pathway resulting in genome methylation was discovered (Ye et al., 2016). In this process, longer RNA molecules bind to AGO4 and are subsequently trimmed, independent of Dicer, by 3′- to 5′-exonucleases into siRNAs. The population of siRNAs resulting from this shows up as a “ladder” pattern (~20–60 nt in size) and shares the same 5′ end but differ at their 3′ end with one nucleotide. Since not much is known yet on the role of members from the RDRγ in (viral) siRNA biogenesis, and whether these would assist in the production of dicer independent siRNAs (sidRNA), the reads collected from MM and Ty-1 tomato were examined for molecules larger than 24-nt and matching features of sidRNAs. In both MM and Ty-1 tomato TYLCV-specific vsiRNAs >24 nt were found, but the amounts of these were not substantially increased in the presence of Ty-1 (Supplementary Figure S5A). Spots with potential sidRNAs were found by screening for vsiRNAs 25-33 nt in size and sharing the same 5′ end (for an example of a potential sidRNA refer to Supplementary Figure S5C). The potential sidRNA spots were mostly in antisense orientation and some were even consistently found in biological replicates, although the positions identified in Ty-1 plants were different from those in MM plants (Supplementary Table S1, Supplementary Figures S5D–F). However, in general, their number was relatively low, i.e., 11 and 3 from MM and Ty-1 plants, respectively (Supplementary Figure S5B).