Analysis of Virus-Derived siRNAs in Strawberry Plants Co-Infected with Multiple Viruses and Their Genotypes

Plants can be infected with multiple viruses. High-throughput sequencing tools have enabled numerous discoveries of multi-strain infections, when more than one viral strain or divergent genomic variant infects a single plant. Here, we investigated small interfering RNAs (siRNAs) in a single strawberry plant co-infected with several strains of strawberry mottle virus (SMoV), strawberry crinkle virus (SCV) and strawberry virus 1 (StrV-1). A range of plants infected with subsets of the initial viral species and strains that were obtained by aphid-mediated transmission were also evaluated. Using high-throughput sequencing, we characterized the small RNA fractions associated with different genotypes of these three viruses and determined small RNA hotspot regions in viral genomes. A comparison of virus-specific siRNA (vsiRNA) abundance with relative viral concentrations did not reveal any consistent agreement. Strawberry mottle virus strains exhibiting considerable variations in concentrations were found to be associated with comparable quantities of vsiRNAs. Additionally, by estimating the specificity of siRNAs to different viral strains, we observed that a substantial pool of vsiRNAs could target all SMoV strains, while strain-specific vsiRNAs predominantly targeted rhabdoviruses, SCV and StrV-1. This highlights the intricate nature and potential interference of the antiviral response within a single infected plant when multiple viruses are present.


Introduction
RNA silencing is a conserved mechanism invoked not only to counteract invading viral and viroid species, but also to regulate gene expression [1]. Double-stranded RNA (dsRNA) synthesized during viral replication, as well as highly structured single-stranded RNA of the viral genome can trigger RNA silencing pathways [1][2][3]. Host DICER-like enzymes cleave triggering viral RNA templates into 21-24 nucleotide-long virus-derived small interfering RNA (vsiRNA) products [4]. As a component of RNA-induced silencing complexes, after vsiRNAs are loaded, they are used as blueprints for sequence recognition, ensuring precise degradation and/or translational inhibition of the targeted RNAs [1,4]. Furthermore, vsiRNAs bind to target RNAs during viral infection and direct RNA-dependent RNA polymerases to synthesize new dsRNA products; this in turn produces transitive or secondary vsiRNAs and amplifies silencing potency in plants [5,6].
Analyzing the vsiRNA pattern of virus-infected plants has shown that specific genomic regions give rise to greater numbers of vsiRNAs, termed RNA silencing 'hotspots' [7]. Stolons from the original plant were used for vegetative propagation. Three randomly selected daughter plants were used for small RNA sequencing library preparation. From each plant, three to four libraries were prepared from different leaves, resulting in eleven datasets (1B-four datasets (1B_1-1B_4), 2A-three datasets (2A_2-2A_4), and 3K-four datasets (3K_1-3K_4)). The samples were verified in parallel with RT-qPCR (Table S1, detailed view). Additionally, aphid-mediated transmission was used to produce strawberry plants infected with only a subset of viruses. From these plants, eleven additional small RNA sequencing libraries were prepared to evaluate plant vsiRNA specificity against individual viral strains (Table S1, datasets). The obtained data were deposited at the NCBI repository under accession number SAMN34995356-77.
Among 10 million randomly selected trimmed reads, 29,832-69,742 could be mapped to the genome of any of the infecting viral strains, with virus-specific reads comprising less than 1% of the total small RNA population (Table S1). Of these viral reads, 82-89% were found to be strain-specific and mapped to only one of the presenting viral genomes.
Certain reads were not unique and mapped to more than one strain. The exclusion of such reads resulted in a smaller set of strain-specific reads (Table 1). This difference was the largest in the case of SMoV, which has long 3′ untranslated regions (3′ UTRs) that are moderately conserved not only between its two genomic segments, but also between the different strains included in the current study ( Figure 1A,B, lower triangular parts). These regions showed a notably higher vsiRNA coverage than the remaining genomic parts. To some extent, this might be caused by the presence of secondary structures, which are known to play an important role in siRNA genesis [49] and have an indispensable role in the viral replication of RNA viruses by serving as 3′ cap-independent translation enhancers [50,51]. Stolons from the original plant were used for vegetative propagation. Three randomly selected daughter plants were used for small RNA sequencing library preparation. From each plant, three to four libraries were prepared from different leaves, resulting in eleven datasets (1B-four datasets (1B_1-1B_4), 2A-three datasets (2A_2-2A_4), and 3K-four datasets (3K_1-3K_4)). The samples were verified in parallel with RT-qPCR (Table S1, detailed view). Additionally, aphid-mediated transmission was used to produce strawberry plants infected with only a subset of viruses. From these plants, eleven additional small RNA sequencing libraries were prepared to evaluate plant vsiRNA specificity against individual viral strains (Table S1, datasets). The obtained data were deposited at the NCBI repository under accession number SAMN34995356-77.
Among 10 million randomly selected trimmed reads, 29,832-69,742 could be mapped to the genome of any of the infecting viral strains, with virus-specific reads comprising less than 1% of the total small RNA population (Table S1). Of these viral reads, 82-89% were found to be strain-specific and mapped to only one of the presenting viral genomes.
Certain reads were not unique and mapped to more than one strain. The exclusion of such reads resulted in a smaller set of strain-specific reads (Table 1). This difference was the largest in the case of SMoV, which has long 3 untranslated regions (3 UTRs) that are moderately conserved not only between its two genomic segments, but also between the different strains included in the current study ( Figure 1A,B, lower triangular parts). These regions showed a notably higher vsiRNA coverage than the remaining genomic parts. To some extent, this might be caused by the presence of secondary structures, which are known to play an important role in siRNA genesis [49] and have an indispensable role in the viral replication of RNA viruses by serving as 3 cap-independent translation enhancers [50,51].

Comparison of Viral Concentrations and vsiRNA Abundances
Using biological replicates of CRM plants, significant differences were established between the relative RNA concentrations of some, but not all of the strains of the studied viruses (Figure 2A). At the same time, the differences in the numbers of vsiRNA ( Figure 2B) did not fully correspond to the differences in viral concentrations.  Table S2.
The numbers of mapped reads for SCV and StrV-1 did not differ much when comparing total and strain-specific mapped reads. However, due to the long conserved 3′ UTRs, the difference for SMoV was twofold on average (Tables 1 and S1). This was clear based on the comparison of plots for all mapped ( Figure 3A), as well as only strain-specific ( Figure 3C) reads. Only statistically significant differences between strains of the same species are denoted by square brackets with asterisks (* p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). Relative viral concentration data are listed in Table S2.
For SMoV RNAs, significant differences in both the concentrations and plant vsiRNA quantities were observed only between viral RNA1 A and RNA1 B . The relative concentrations of RNA2 segments were higher than those of RNA1s segments, but the corresponding small RNAs showed the opposite trend: the number of vsiRNA reads generated from RNA1 variants was higher than that generated from RNA2 variants. This opposite correlation can be explained by the fact that vsiRNAs are not only produced from, but also target viral RNAs. Hence, their increased level might result in a greater abundance of SMoV-specific loaded Argonaute proteins, which target and cleave the corresponding viral RNAs, resulting in a decrease in their concentrations.
Both strains of SCV had comparable viral concentrations and numbers of strainspecific vsiRNAs. In contrast, for the other rhabdovirus, StrV-1, significant differences in the concentrations and vsiRNAs between all three strains were observed, with the only exception being nonsignificant changes in the levels of vsiRNAs between strains A and C. The numbers of mapped reads for SCV and StrV-1 did not differ much when comparing total and strain-specific mapped reads. However, due to the long conserved 3 UTRs, the difference for SMoV was twofold on average (Tables 1 and S1). This was clear based on the comparison of plots for all mapped ( Figure 3A), as well as only strain-specific ( Figure 3C) reads.

Figure 2.
Boxplots of relative viral concentrations ((A), log 10 scale) and mapped reads (B) calculated for eleven viral references (axis of abscissas) from three CRM daughter plants. Only statistically significant differences between strains of the same species are denoted by square brackets with asterisks (* p ≤ 0.05, ** p ≤ 0.01, *** p ≤ 0.001). Relative viral concentration data are listed in Table S2.
The numbers of mapped reads for SCV and StrV-1 did not differ much when comparing total and strain-specific mapped reads. However, due to the long conserved 3′ UTRs, the difference for SMoV was twofold on average (Tables 1 and S1). This was clear based on the comparison of plots for all mapped ( Figure 3A), as well as only strain-specific ( Figure 3C) reads.  Another irregularity was observed in the StrV-1 and SCV mappings. The StrV-1B strain shows a hotspot region at the 5 end of the genomic RNA ( Figure 4) that is absent in the A and C strains, but invariably present in all samples ( Figures S1 and S3).
A closer inspection of the peak showed that it predominantly contains 27 nt long reads mapped to the terminus. In an RNA silencing study of tomato yellow leaf curl Sardinia virus, a begomovirus with a DNA genome, Miozzi et al. identified host-derived siRNAs that cross-react with the major virus hotspot [52], and these siRNAs have been suggested to originate from integrated ancient remains of geminiviral-related DNA that is part of the RNA silencing machinery [52,53]. To determine whether this might explain our results, a BLASTn search against a Fragaria GenBank database (taxid:3746) was performed but did not return a 100% match, indicating a non-Fragaria origin of this sequence (searched on 12 November 2022). For StrV-1, further analyses showed that in addition to the plants with multiple infections, some StrV-1-negative samples (see Section 2) also contained these 27 nt long reads, yet in much lower quantities ( Figures S1 and S3). It may be speculated that such reads were assigned to the sample due to index hopping or read misidentification [54], as library construction was performed with a single index [54,55]. Nevertheless, to obtain a clear explanation, further investigation is needed.
(mapping) coverage. The genomic organization of the reference is shown at the bottom side of the figure. The boxed rectangle denotes the open reading frame. Distribution plots for all samples are shown in Figure S1.
Another irregularity was observed in the StrV-1 and SCV mappings. The StrV-1B strain shows a hotspot region at the 5′ end of the genomic RNA ( Figure 4) that is absent in the A and C strains, but invariably present in all samples ( Figures S1 and S3). A closer inspection of the peak showed that it predominantly contains 27 nt long reads mapped to the terminus. In an RNA silencing study of tomato yellow leaf curl Sardinia virus, a begomovirus with a DNA genome, Miozzi et al. identified host-derived siR-NAs that cross-react with the major virus hotspot [52], and these siRNAs have been suggested to originate from integrated ancient remains of geminiviral-related DNA that is part of the RNA silencing machinery [52,53]. To determine whether this might explain our results, a BLASTn search against a Fragaria GenBank database (taxid:3746) was performed but did not return a 100% match, indicating a non-Fragaria origin of this sequence (searched on 12 November 2022). For StrV-1, further analyses showed that in addition to the plants with multiple infections, some StrV-1-negative samples (see Section 2) also contained these 27 nt long reads, yet in much lower quantities (Figures S1 and S3). It may be speculated that such reads were assigned to the sample due to index hopping or read misidentification [54], as library construction was performed with a single index [54,55]. Nevertheless, to obtain a clear explanation, further investigation is needed. Unlike StrV-1, SCV exhibits hotspot vsiRNA regions near the 3 end of the genomic RNA, as detected for both strains ( Figure 5) across all samples ( Figure S1). This region is not conserved between A and B strains of SCV, meaning that the mapped reads are unique to both strains. Furthermore, there are no mapped reads in the SCV-negative samples.

The Abundance of vsiRNAs against Individual Strains of SMoV Does Not Correspond to Their Relative Concentration
During our previous study on aphid-and petiole-wedge grafting-mediated transmission of SMoV, SCV and StrV-1, we obtained a set of Fragaria vesca plants infected with either individual viruses and strains or their various combinations. During identification of viral strain composition after aphid-mediated transmission, we identified two plants infected with a combination of viral SMoV strains: RNA1 BC and RNA2 AC , and found that RNA1 B showed a significantly lower concentration than RNA1 C ( Figure 6A). Although we observed differences between different SMoV RNAs during analyses of other samples, in this case, the relative concentration of RNA1 B was on average 500-fold lower (Table S2, Figure 6A). Unlike StrV-1, SCV exhibits hotspot vsiRNA regions near the 3′ end of the genomic RNA, as detected for both strains ( Figure 5) across all samples ( Figure S1).   Figure S1. either individual viruses and strains or their various combinations. During identification of viral strain composition after aphid-mediated transmission, we identified two plants infected with a combination of viral SMoV strains: RNA1 BC and RNA2 AC , and found that RNA1 B showed a significantly lower concentration than RNA1 C ( Figure 6A). Although we observed differences between different SMoV RNAs during analyses of other samples, in this case, the relative concentration of RNA1 B was on average 500-fold lower (Table S2, Figure 6A). Therefore, we hypothesized that the RNA1 B -and RNA1 C -specific vsiRNAs should also quantitatively vary. However, there was no significant difference in the numbers of mapped reads ( Figure 6B Taking into account the considerable difference in the titers of RNA1 B and RNA1 C , it would be problematic to detect RNA1 B without prior knowledge of its existence. Nevertheless, as the actual variability of SMoV strains was uncovered only a few years ago, reliable tools for the identification of the full spectrum of SMoV strains may still be lacking. The existing SMoV detection system relies on the highly conserved stretch within the Therefore, we hypothesized that the RNA1 B -and RNA1 C -specific vsiRNAs should also quantitatively vary. However, there was no significant difference in the numbers of mapped reads ( Figure 6B Thus, future studies of SMoV strain variability should further examine these findings. Taking into account the considerable difference in the titers of RNA1 B and RNA1 C , it would be problematic to detect RNA1 B without prior knowledge of its existence. Nevertheless, as the actual variability of SMoV strains was uncovered only a few years ago, reliable tools for the identification of the full spectrum of SMoV strains may still be lacking. The existing SMoV detection system relies on the highly conserved stretch within the 3 UTR of both RNA1 and RNA2, which is a superb region from a high conservation perspective, but a poor choice for strain discrimination.

Analysis of vsiRNA against StrV-1 Strains in StrV-1-Infected Strawberry
Considering the observed 27 nt hotspot in the StrV-1_B strain, we analyzed the small RNA population in plants infected with either StrV-1_A or StrV-1_B, or a combination thereof ( Figure S3). On average, more than 90% of the mapped reads in this study were recognized as either A-or B-strain-specific. Interestingly, although 21 nt and 22 nt vsiRNAs were dominant, all three datasets from single StrV-1_B infections showed a considerable 27 nt fraction among the mapped reads. After a close examination, it was revealed to consist of a single sequence 5 -ATTGATCGTATAGATGTTATCATCCGT-3 aligned against the StrV-1 genomic strand at its 5 end (positions 2-28). An analysis of potential secondary structures of the sequence showed only a weak structure with a four-base-long loop (underlined). The genesis of vsiRNA involves the production of mainly 21 nt and 22 nt long products, and the presence of 24 nt vsiRNAs has been reported only for some DNA and RNA viruses [56,57]. During the peer-reviewing of the current manuscript, one of the reviewers pointed out that this RNA might be a defective viral RNA. Defective interfering RNAs (DI RNAs) are produced during several passage of viral infection [58]; hence, it would be interesting to further investigate the viral derived RNA population for the presence of such DI elements. Moreover, siRNAs produced from DIs can saturate the siRNA binding capacity of viral VSRs, leading to the attenuation of viral symptoms [59], which would further alter the complex picture under co-infection. However, HTS in combination with validation via Sanger sequencing would be an appropriate tool for DI identification, as the building of chimeric contigs could generate false DI annotations. To determine if the detected 27 nt long small RNAs are or could act as DI RNAs is an interesting question that could be addressed in the future. Therefore, assuming that the 27 nt small RNAs are not sequencing artifacts, they might be produced in other pathways or could be an intermediate product of vsiRNA processing.
Furthermore, comparing the number of vsiRNAs per virus in the different samples showed that there were 14,000 mapped reads on average in a single StrV-1 infection, while the corresponding average when other viruses were present was much lower, at 3500 (Table S1). SMoV encodes two weak viral suppressors of RNA silencing, which were able to alter the potato virus X concentration in a co-expression study [60]. Thus, it is possible that SMoV coinfection with rhabdoviruses influences the RNA silencing reaction against them through this mechanism.

Summary and Conclusions
Our data show that plants infected with more than one viral strain produces strainspecific vsiRNAs against each viral RNA. For SCV, the vsiRNA hotspot was found at the 3 terminus of the genomic RNA, whereas no universal hotspot was found for StrV-1. Furthermore, although both the SCV and StrV-1 vsiRNAs were quite specific, a large fraction of SMoV vsiRNAs matched more than one strain due to its long conserved 3 UTR regions. It seems that in the development of a plant protection method based on exogenous siRNAs against SMoV, 3 UTRs might be a good target. However, there was no consistent agreement between the relative levels of viral RNA and the corresponding levels of vsiRNAs. Despite the very low abundance of SMoV RNA1B in Fv-SMoV_BC_AC-1 and -2 plants, the levels of RNA1B-specific vsiRNAs were not significantly lower than those of vsiRNAs targeting the other strain or other viruses. When loaded into the RISC, vsiRNAs in multiple infections can target related viral strains, which is why it is very difficult to find a direct correlation between the concentration of different viral strains and the number of vsiRNAs. The presence of viral silencing suppressors capable of interfering with antiviral silencing in different steps of RNAi can further alter this complex pattern, which should be investigated in future studies.

Plant Materials
The CRM3 isolate of F. ananassa, cultivar Cacanska rana, served as a plant source for stolon-mediated propagation ( Table 2,  The CRM daughter plants served as virus sources during the aphid-mediated transmission of SMoV and StrV-1 to Alpine strawberry, F. vesca semperflorens; for details, see Koloniuk et al. [44]. Briefly, F. vesca plants infected with some of the abovementioned strains of SMoV and StrV-1 were obtained ( Table 2). The mixed infection of StrV-1-A and StrV-1-B strains was obtained using the consequent grafting-mediated infection [44].
The plants were maintained in an insect-proof greenhouse with a controlled temperature, under a 16 h light/8 h dark photoperiod.

Total and Small RNA Isolation
Combined total plant RNA and small RNA extractions were performed using a single leaf that was snap-frozen in liquid nitrogen, homogenized, divided into two approximately 50 mg aliquots and processed with a Thermo Scientific GeneJET Plant RNA Purification Mini Kit (Thermo Scientific, Waltham, MA, USA) and mirPremier microRNA Isolation Kit (Sigma-Aldrich, St. Louis, MI, USA), respectively, following manufacturers' recommendations. The quantification and quality control of the RNA extracts were performed using a Nanodrop 1000 UV-Vis spectrophotometer and Qubit HS RNA and IQ assays (Invitrogen, Carlsbad, CA, USA).

cDNA Synthesis and Two-Step Reverse-Transcription qPCR (RT-qPCR)
Total RNA (200 ng-400 ng of the total plant RNA) was reverse-transcribed to cDNA using the Maxima H Minus First Strand cDNA Synthesis Kit with dsDNase in 20 µL reactions following the manufacturers' recommendations. The cDNA was then diluted to a ratio of 1:10 with Milli-Q-grade water and subjected to qPCR assays.
The RT-qPCR assays were conducted using a CFX96 real-time PCR detection system (Bio-Rad, Hercules, CA, USA). The 10 µL reaction was prepared from 5 µL of tenfolddiluted cDNA, 0.25 µL forward and reverse primers (10 mM, final concentration 250 nM), 2.75 µL of nuclease-free water and 2 µL of 5x HOT FIREPol EvaGreen qPCR Mix Plus (Solis BioDyne, Taru, Estonia). The reaction conditions were 95 • C for 12 min, followed by 40 cycles of 95 • C for 10 s, 60 • C for 20 s and 72 • C for 20 s. The dissociation curve analysis was performed by ramping from 65 • C to 95 • C (with increments of 0.5 • C for 5 s) to verify the specificity of primer amplification and the presence of potential primer dimers based on a single peak. No template, positive, and if necessary, no reverse transcriptase controls were included to check for potential cross-contamination and the presence of genomic DNA. The amplification efficiency (E , Table S3) was assessed using a standard curve based on serial dilutions of the cDNA template. Each reaction was carried out in triplicate.
Relative viral concentrations were calculated using the formula efficiencyˆ(Cq ref − Cq virus ), where efficiency is the experimentally calculated efficiency of the corresponding primers (Table S3), Cq ref is the geometric mean of the Cq values of three endogenous reference mRNAs (encoding DNA-binding protein, histone H4 and pyruvate decarboxylase) and Cq virus is the Cq value of the viral strain.

High-Throughput Sequencing and Data Analysis
The plants and respective sequencing libraries are listed in Table 2. The sequencing libraries were prepared from small RNA extracts following the manufacturer's recommendations. The NEBNext Multiplex Small RNA Library Prep Kit for Illumina was used. The ready-to-load libraries were processed using either the NovaSeq6000 or HiSeq 2500 system. Raw reads were quality-, library adapter-(5 -AGACGTGTGCTCTTCCGATCT-3 ), and length-trimmed (18-30 nt selection) and a random selection of either 7 or 10 million reads was performed (Table S1, detailed view). The obtained reads were then mapped to eleven viral sequences using the 'Map to reference' tool of Geneious (Biomatters, Inc., Auckland, New Zealand).
Specifically, the minimum identity was set to 100%, and the reads matching more than one reference sequence were mapped either to all or none of them, depending on the required result. Assuming the presence of more than one strain in the data, we set mapping conditions that were rather stringent, without allowed mismatches. The mapping settings allowed us to establish rules for a read matching more than one position (reference): (a) discard the read, (b) map it to all matching positions or (c) map it to only one position through a random algorithm. In the current study, we used the first two options. When such reads were mapped to all references, they were referred to as 'mapped reads'; when such reads were discarded, thus retaining only reference-unique reads, they were referred to as 'strain-specific' ( Table 1, Figures S1 and S3).

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants12132564/s1, Table S1: Overview of the mapped reads to viral references: strawberry crinkle virus, strawberry mottle virus and strawberry virus 1; Table S2: Overview of Cq, ∆Cq and fold difference values between the viral targets and endogenous references; Table S3: List of primers used in the study; Figure S1: Mapping and distribution plots of vsiRNAs (reads) from the 1B, 2A and 3K samples; Figure S2: Mapping and distribution plots of vsiRNAs (reads) from the Fv-SMoV_BC_AC samples; Figure S3: Mapping and distribution plots of vsiRNAs (reads) from the Fv_StrV-1_A, Fv_StrV-1_B and Fv_StrV-1_AB samples. Reference [62] is cited in the supplementary materials.