Sugar beet plants at the 5–6 leaf-stage inside cages were exposed to beet leafhoppers (BLH) by releasing viruliferous BLHs (approximately 6–8 BLH/plant) carrying predominantly the BCTV Severe (BCTV-Svr) and BCTV CA/Logan strains. The uninfected control plants were similarly transferred into cages but without any BLHs. Growth chamber conditions were 28°C (day)/21°C (night), 16 h (day)/8 h (night) photoperiod, and 20% relative humidity. Newly emerged leaves from both control (uninfected) and BLH infected plants were collected at 1, 2, and 6 day post inoculation (dpi), flash frozen in liquid N, and stored at –80°C until further use. After sample collections at each time point, the plants were sprayed with Admire^® Pro (Bayer CropScience LLC, NC, United States) insecticide to eliminate any BLHs and future pest infections. This was followed by spraying with water to remove any residual insecticide, and then moved to the green house for symptom development in the infected plants in the subsequent weeks. The plants were evaluated for curly top symptoms at 3 week post inoculation.

Total RNA was extracted using “Plant/Fungi Total RNA Purification Kit” (Norgen Biotek Corp, ON, Canada) according to the manufacturer’s protocol. The quality and quantity of total RNA were analyzed using Bioanalyzer 2100 (Agilent Technologies, CA, United States) and acceptable RIN number for samples >7.0. High quality RNA was used for sRNAseq and mRNAseq library preparations. For sRNA libraries, approximately 1 μg of total RNA was used according to the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, United States). This was followed by 50 bp single-end sequencing on an Illumina Hiseq 2500 sequencing platform at the LC Sciences (Houston, TX, United States) following the vendor’s protocol.

Poly(A) RNA sequencing libraries were prepared according to the Illumina’s TruSeq-stranded-mRNA sample preparation protocol. Approximately 1 μg of total RNA was used for this purpose. Depletion of ribosomal RNA were performed following the protocol described in the Ribo-Zero™ rRNA Removal Kit (Illumina, San Diego, CA, United States). Poly(A) mRNAs were purified using oligo-(dT) magnetic beads following two rounds of purification. Fragmentation of poly(A) RNA was performed using divalent cation buffer under elevated temperature. Cleaved RNA fragments were then reverse-transcribed to produce cDNA that was used to produce U-labeled second strand DNA. Following end repair, 3′ adenylation, adapter ligation, and PCR, final libraries were produced. Quality control and quantification of the sequencing libraries were done using Bioanalyzer 2100 (Agilent Technologies, CA, United States). Paired-end (150 bp) sequencing was performed on Illumina’s NovaSeq 6000 sequencing platform.

Raw reads were processed using an in-house program, ACGT101-miR (LC Sciences), that removed sequences associated with adapter dimers, common RNA families [ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), small nucleolar RNA (snoRNA)], low quality reads, and repeats. Unique sequences of 18--25 bp in length were mapped to specific species precursors in the miRBase 22.0 database¹ using BLAST to identify known miRNAs and novel miRNAs. The alignment took into consideration any length variation at both 3′ and 5′ ends and one mismatch inside the sequence. Known miRNAs were identified as sequences mapping to specific species mature miRNAs in the hairpin arms. The sequences mapping to the other arm of known specific species precursor hairpin opposite to the annotated mature miRNA-containing arm, were considered as novel 5p- or 3p derived miRNA candidates. The remaining sequences were mapped to other selected species precursors (with the exclusion of specific species) in the miRBase 22.0 database using BLAST search, and the mapped pre-miRNAs were further BLAST searched against the sugar beet genome to determine genomic locations. The unmapped sequences were BLAST searched against the sugar beet genome, and the hairpin RNA structures containing sequences were predicted from the flanking 120 nt sequences using the software RNAfold.² A global normalization method was used to rectify copy numbers among different samples (Li et al., 2016). This helps to overcome any effect resulting from sequencing discrepancy on miRNA expression among samples. miRNAs that showed p-value <0.05 and fold change ≥1.5, were designated as Differentially Expressed (DE) miRNAs. Genes targeted by most abundant miRNAs was performed using computational target prediction algorithm PsRobot 1.2 (Wu et al., 2012). Predicted miRNA target genes were further validated using mRNAseq data by comparing expression of specific miRNAs to the expression of corresponding target genes.

Raw reads originating from mRNA-Seq data were processed using in house (LC Sciences) Perl scripts and Cutadapt (Martin, 2011) to remove any low-quality reads and adapter sequences. FastQC³ was performed to analyze sequence quality. Reads were mapped to the sugar beet (Beta vulgaris subsp. vulgaris) reference genome EL10⁴ using Hisat2 (Kim et al., 2015) and mapped reads of each sample were assembled using StringTie (Pertea et al., 2015).

All transcriptome data originating from different samples were merged to reconstruct a comprehensive transcriptome data using Perl scripts. Following merging, StringTie (Pertea et al., 2015) and edgeR (Robinson et al., 2010) software were used to estimate the expression levels of transcripts. Expression level of mRNAs was performed by calculating FPKM (fragments per kilobase of transcript per million mapped reads) using the StringTie software. EdgeR-R packages were used to estimate differential expression. A cutoff p-value of <0.05 and | log2 (fold-change)| > 1 was set to identify differentially expressed genes (DEGs).

Annotation of transcripts were performed using Blastx against the NCBI database. For Gene Ontology (GO) analysis, transcripts were blasted to the GO database to calculate the gene numbers for each term. Pathway enrichment was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG; Kanehisa et al., 2008).

RNAseq reads, detailed previously, were aligned to a single consensus sequence from four BCTV genomes (Strausbaugh et al., 2017) which represent the important strains of BCTV in sugar beet: Severe strain CTS06-021 (GenBank accession KX867019.1), Colorado strain CTS15-113 (KX867056.1), CA/Logan strain CTS06-104 (KX867032.1) and the Worland strain CTS15-095 (KX867055.1). Alignment of reads was carried out using BWA (Li, 2013). The resulting files (.sam) were converted to (.bam) and sorted using Samtools v1.9 (Li et al., 2009). Only reads with high quality alignments (Q > 30) to the virus genome were retained for downstream analysis. These reads were converted to (.fasta) format. Each sequence name was derived from the sequence ID generated from the sequencer. The reads in (.fasta) format were reverse complemented using a short python script. The program blastn (Altschul et al., 1990) was carried out against the EL10.1 CDNA sequences. Blast results were parsed for information about core sequence homology of putative sncRNA species to the host sugar beet genome EL10.1. Specific sncRNA species were investigated further based on sequence length and percent identity between the core sncRNA species and their sugar beet genome targets. mRNA expression was evaluated for targets of sncRNA species by looking at FPKM values of the target genes. Using simple bash scripts to extract this information, statistical tests between treatments, and correlations between sncRNA abundance and gene expression (FPKM) of host target genes were carried out.

The virus genome was visualized using the “R” software package. Virus gene positions were extracted from the Severe strain CTS06-021 (GenBank accession KX867019.1). Samtools depth was used to determine the coverage of reads across the genome for all time points. Read depth was plotted across the viral genome with respect to functional elements within the genome (e.g., viral genes and regions which produced specific sncRNA). Data were transformed into circular plots using functions found in the “Rcircos” source code (Zhang et al., 2013).

Both uninfected control and BCTV inoculated plants at 1, 2, and 6 dpi were evaluated for any BCTV symptom development in the leaves and representative plants from each timepoint were photographed. Qualitative/quantitative analysis of viral load in the plant samples were determined through RT-PCR. cDNA was prepared using the iScript™ gDNA Clear cDNA Synthesis Kit (Bio-Rad, CA, United States) according to the manufacturer’s protocol. The thermocycler conditions included a pre-incubation step at 98°C for 30 s followed by 33 cycles of 98°C for 5 s (enzyme activation), 60°C for 10 s (primer annealing), 72°C for 3 s (elongation), and followed by a final elongation step at 72°C for 5 min. Equal amount of cDNA template was added to each sample for PCR amplification followed by gel electrophoresis of equal volume of the PCR products in a 2% agarose gel. The primers used for the amplification of BCTV capsid protein gene (V1), qCP_F-5′-TCCCAGAAAAGAAAGGTGAATCCT-3′ and qCP_R-5′-CCATTGGTATTTCCTCGAAGTCGT-3′, were universal for all 4 BCTV strains. The sugar beet house-keeping gene [glutamine synthetase (GNL2; EU558132.1)] primers were qBv_F-5′-CCTTCAGGGTGTCGCCAAT-3′ and qBv_F-5′-CCTTCAGGGTGTCGCCAAT-3′.

Soluble sugars were extracted from approximately 100 mg of finely ground sugar beet leaf tissues (stored in –80°C) and analyzed by the method described here. Briefly, soluble sugars were extracted in 1 ml of 80% absolute ethanol at 65°C for 30 min. Samples were removed from heat, cooled for 5 min., then vortexed at medium speed for 2 min. and centrifuged at 13,000 rpm for 8 min to pellet tissue. Extracts were filtered through a 0.45 μm nylon syringe filter (Pall Corp., Port Washington, NY, United States) fitted onto a 3 ml syringe (Becton, Dickinson and Company, Franklin Lakes, NJ, United States). The profiles for 11 soluble sugars were determined using a PerkinElmer series 200 HPLC pump and autosampler (PerkinElmer Inc., Waltham, MA, United States) coupled with a Shimadzu RID-10A refractive index detector set at 30°C (Shimadzu Scientific Instruments Inc., Columbia, MD, United States). Sugars were separated on a Luna NH2 analytical column (heated to 25°C, 250 mm × 4.6 mm, 5 μm, 100 Å, Phenomenex Inc., Torrance, CA, United States) outfitted with an NH2 Guard column (4 mm × 3 mm, Phenomenex Inc.) using an isocratic mobile phase of 80% acetonitrile at a flow rate of 2 ml min^–1. Besides Sucrose, other identified sugars in sugar beet leaves were quantified using a 5-point external standard curve (0.125–2 mg ml^–1). For quantifying Sucrose, the standard curve was expanded from 2 to 18 mg ml^–1. To quantify glucose + galactose (two peaks that did not separate*), the areas and concentrations of each were added together to create a combined standard curve. The chromatographs were analyzed, and the data were processed using PerkinElmer TotalChrom software (version 6.2.1).

The raw data resulting from sRNAseq (BioProject ID: PRJNA764694) and mRNAseq (BioProject ID: PRJNA764690) were submitted to the NCBI SRA database. All codes used in these steps are documented and available at ‘‘GitHub’’.⁵

On an average we obtained >20 million raw reads/sample and >12 million valid reads/sample (Supplementary Table 1). The size distribution of sRNAs varied between 18 and 26 bp with the highest number of sRNAs were observed around 24 bp (Figure 1A). The nucleotide biasedness of sugar beet miRNAs revealed highest abundance of A, C, G, U at the nucleotide positions 1, 24, 23, and 2, respectively (Figure 1B).

The total number of DE miRNAs that were unique to a specific line decreased from 1 to 6 dpi in the susceptible line 19 and in the resistant line 4 (Figure 2). Whereas, in the resistant line 13, the number of DE miRNAs increased over time. The changes observed in the control plants were relatively smaller compared to the infected plants (Supplementary Figure 1). Heatmaps of DE miRNAs (p < 0.01) at 1, 2, 6 dpi are shown in Figure 3 and the control samples in Supplementary Figure 2. All DE miRNAs (p < 0.05) in the infected and uninfected plants from the S and R lines are presented in Supplementary Tables 2–4. Examples of moderately expressed miRNAs at 1 dpi only in the R lines and absent in the S line are ath-miR169b-5p_R + 3, PC-3p-24046_2377, and PC-3p-85965_890 (Figure 3A and Supplementary Table 2). Some of the high to moderately expressed miRNAs that were down-regulated (>2-fold) only in the R lines (vs. S) are csi-miR164a-3p_2ss9CT10TG, ath-miR172c_R + 1, and PC-5p-93017_834. At 2 dpi, examples of miRNAs whose expression were high to moderate and were upregulated (∼1.5 to 50-fold) in the R lines (vs. S) are ath-miR166a-3p_1ss20CT, stu-MIR8005a-p5_1ss5GA, PC-5p-3041_10997, PC-5p-14805_3385, and PC-3p-62242_1151 (Figure 3B and Supplementary Table 3). Micro RNAs whose expression were high to moderate and were downregulated (∼2 to 30-fold) in the R lines (vs. S) are csi-miR167a-5p_R + 2, csi-miR156b-3p_R-1, and PC-5p-132263_616 to name a few. At 6 dpi, examples of miRNAs that were moderately expressed only in the R lines include PC-3p-93312_831 and PC-3p-85965_890 (Figure 3C and Supplementary Table 4). miRNAs whose expression were high to moderate and were upregulated (∼5 to 7-fold) in the R lines (vs. S) include PC-3p-777_32390, PC-5p-14805_3385, and PC-5p-3041_10997 to name a few. Examples of miRNAs that were down-regulated (∼2.5 to 27-fold) in the R lines (vs. S) are csi-miR164a-3p_2ss9CT10TG, ath-miR160a-5p_R-1, and PC-5p-205203_404.

Gene ontology (GO) analysis of genes targeted by DE miRNAs (considering all time points) due to BCTV infection represented genes that were mainly associated with cellular processes including transcription, protein phosphorylation, oxidation-reduction to name a few (Figure 4A). Pathway enrichment analysis of genes targeted by DE miRNAs (considering all time points) belonged to the pathways predominantly associated with glycolysis/gluconeogenesis, galactose metabolism, starch and sucrose metabolism, ribosome biogenesis, pyruvate metabolism etc. (Figure 4B and Supplementary Figure 3).

We also analyzed potential sugar beet gene targets of DE sugar beet miRNAs. An extensive list of the target genes is presented in Supplementary Table 5. This list was further narrowed down to those target genes whose expression were altered by the DE miRNAs in the R vs. S lines. For mRNAseq, we obtained >32 million raw reads/sample and >31 million valid reads/sample (Supplementary Table 6). In several cases a single miRNA targeted multiple genes and reduced their expression. At 2 dpi, examples of miRNAs that were upregulated in the R lines (Table 1) and reduced target gene expression are lja-miR1511-3p_2ss5AT24CG targeting EL10Ac7g16579 (SPX domain-containing membrane protein; 11 to 20-fold downregulation), PC-5p-48966_1387 targeting EL10Ac6g14610 (cyclin-dependent kinase; 4 to 5-fold downregulation), stu-MIR8005a-p3_2ss15GT21AG targeting EL10Ac9g22557 (ubiquitin-like domain-containing CTD phosphatase; 2 to 400-fold downregulation), stu-MIR8005a-p5_1ss5GA targeting EL10Ac1g00203 (transcription factor bHLH14; 2 to 5-fold downregulation), and stu-MIR8005c-p5_2ss14AG23TA targeting EL10Ac2g03749 (transcription factor; >200-fold downregulation). At 2 dpi, some of the downregulated miRNAs in the R lines (Table 2) include ahy-miR167-3p_L-1_2_ss12CT15TC targeting EL10Ac4g09000 (hypothetical_protein; >20,000-fold upregulation), ath-miR396a-5p targeting EL10Ac5g11116 (chaperone protein; 2 to 3-fold upregulation), cas-miR172a-5p_3ss2_TC3GA9AT targeting EL10Ac3g05363 (ABC transporter; 4 to 4.5-fold upregulation), gma-MIR159e-p3_2ss16AG18TA targeting EL10Ac1g02046 (pyruvate dehydrogenase; >80,000-fold upregulation), and gma-MIR159e-p3_2ss16AG18TA targeting EL10Ac9g21768 (probable LRR receptor-like serine/threonine-protein kinase; ∼2-fold upregulation).

Some of the upregulated miRNAs in the R lines at 6 dpi (Table 3) that targeted sugar beet genes include PC-3p-48704_1393 targeting EL10Ac5g12315 (ethylene-responsive transcription factor; 1.5 to 4-fold downregulation), PC-3p-70441_1045 targeting EL10Ac5g12641 (long chain acyl-CoA synthetase; >70-fold downregulation), PC-3p-777_32390 targeting EL10Ac9g21076 (putative disease resistance protein RGA4; >350-fold downregulation), and PC-5p-178120_466 targeting EL10Ac2g02731 (probable carboxylesterase; ∼2.5-fold downregulation). At 6 dpi, miRNAs that were downregulated in the R lines (Table 4) and targeted sugar beet genes include ath-miR160a-5p_R-1 targeting EL10Ac1g01085 (15-cis-phytoene desaturase; ∼370 to 450-fold upregulation), ath-miR160a-5p_R-1 targeting EL10Ac2g02632 (putative disease resistance RPP13-like protein; ∼15 to 25-fold upregulation), ath-miR160a-5p_R-1 targeting EL10Ac7g17778 (probable LRR receptor-like serine/threonine-protein kinase; 15 to 17-fold upregulation), bra-miR164e-3p_1ss14CT targeting EL10Ac2g02793 (putative disease resistance RPP13-like protein; ∼30-fold upregulation), PC-3p-50231_1360 targeting EL10Ac1g01374 (serine/threonine-protein phosphatase; >20,000-fold upregulation), PC-3p-539842_121 targeting EL10Ac6g14637 (heat shock protein; >35,000-fold upregulation), PC-5p-205203_404 targeting EL10Ac1g01085 (15-cis-phytoene desaturase; ∼370 to 450-fold upregulation), and PC-5p-725678_75 targeting EL10Ac8g20096 (LRR-like serine/threonine-protein kinase; ∼6 to 9-fold upregulation).

As miRNAs targeted genes that were highly associated with carbohydrate metabolism, we therefore investigated leaf carbohydrate contents in control and BCTV infected plants of the S and R lines. Cellular content of glucose + galactose was lower (37–54%) in the control plants of the R lines (vs. S) without any change among the lines following infection (Figure 5). Fructose on the other hand increased by 30–35% only in the infected plants of the R lines. No specific trend in leaf sucrose content was observed within R lines post BCTV infection. Sucrose content decreased by 40% in the line 13 (R) whereas, line 4 (R) showed a 44% increase in sucrose content following infection.

We also investigated if any sugar beet miRNAs targeted BCTV genes and if higher expression of specific miRNAs have any effect on target viral gene expression. A complete list of miRNAs and their potential BCTV target genes are listed in Table 5. At 2 dpi miRNAs which were moderately expressed, detected only in the R lines and reduced BCTV target gene expression by >2 to 4-fold include PC-3p-28452_2099 and PC-3p-63465_1134. The miRNAs that showed high expression (FPKM 5,000–14,000) at 2 dpi and reduced target gene expression by 3 to 4-fold were ath-miR172a_1ss1AC and ath-miR162a-3p. At 6 dpi only two miRNAs were detected in sugar beet lines that showed moderate/low expression and potentially targeted BCTV genes. We found another miRNA, PC-3p-24046_2377 that showed 10 bp matches with the BCTV origin of replication and were present at both 3 dpi and 6 dpi only in the R lines. In general, we found a 10–12 bp matches between sugar beet miRNAs and their corresponding BCTV target genes and most of these sequence identities were at the seed region of the miRNAs. A majority of the sugar beet miRNAs targeted BCTV genes were associated with capsid protein (V1) and rolling circle replication initiator protein (C1).

Virus derived sncRNAs were investigated by aligning sequencing reads generated from this experiment to the BCTV genome sequences. The mock, or non-infected samples did not contain reads which aligned to the virus genome, indicating no presence of the virus. The infected samples were variable in regard to the abundance of reads and representation of putative virus sncRNA species. RNA was recovered for 1, 2, and 6 dpi with BCTV. Time and inoculation were large factors in the abundance of virus transcripts. No reads mapped to the virus genome from the non-inoculated samples (data not presented) were detected. Due to low viral RNA reads at 1 and 2 dpi, only data from 6 dpi are presented here.

Reads that mapped to any of the four virus strain genomes produced a total of 30 putative sncRNA sequences. These unique sncRNA sequences contained a core motif which showed significant homology to the sugar beet cDNA sequences. The core sequences were 18–22 nt long apart from BCTV_sncRNA 4, which contained a core sequence length of 26 nt (Table 6). Reads that generated 18–22 nt core sequence were longer (31–41nt). This exceeded the size expected from a traditional siRNA or virus-derived siRNA (vsiRNA) but included a core sequence whose reverse complement showed 100% sequence identity to the host genome and defined as a putative sncRNA (Table 6). For visualization purposes, only sncRNA alignments to the BCTV Severe strain CTS06-021 genome were used (Figure 6). This reduced the complexity of the multi fasta alignments and demonstrated relative distribution of the putative sncRNAs. Correlation analysis between sncRNA abundance and corresponding sugar beet target gene expression (FPKM) helped to validate their potential role in BCTV resistance in sugar beet. The BCTV sncRNAs that showed highest to moderate abundance were BCTV_sncRNA 4, 18, and 19 and had the most down-regulation (1.5 to 2-fold) of target sugar beet gene expression EL10Ac1g01206 (LRR), EL10Ac5g12605 (7-deoxyloganetic acid glucosyltransferase), and EL10Ac6g14074 (transmembrane emp24 domain containing) respectively in the S line (Table 7).

Both control and BCTV infected plants were evaluated for disease symptoms after three weeks post inoculation. The S line from all BCTV exposure times (1, 2, and 6 days) showed strong BCTV symptoms (3 weeks post inoculation) such as swelling of veins, leaf thickening, and leaf curling predominantly in the apical leaves whereas, the R lines had minimal to no disease symptoms at this time (Figure 7). The control plants that were kept in separate cages without exposure to any BLH did not show any disease symptoms (Supplementary Figure 4). Viral load was determined through RT-PCR estimation of the BCTV CP gene expression with respect to the expression of sugar beet house-keeping gene, glutamine synthetase, in BCTV inoculated plants. In general, the infected plants showed BCTV CP gene expression which had relatively higher signal strength in the S line (vs. R lines) at early time points and the trend was similar at 6 dpi though the overall signal strength of the CP gene expression were relatively higher compared to the early time points (Supplementary Figures 5–7). The sugar beet housekeeping gene showed strong expression across all samples.

The DE miRNAs belonged to different categories including miRNAs that are conserved in plants and miRNAs unique to sugar beets. Upregulation of specific genes (>70 to 200-fold) such as EL10Ac2g03749 (transcription factor) and EL10Ac5g12641 (long chain acyl-coA synthetase) in the S line (vs. R) at 2 dpi/6 dpi (Tables 1, 3) due to downregulation of corresponding miRNAs raises the question if such higher expression of target genes in the S line contribute to susceptibility. On the other hand, downregulation of specific miRNAs such as ath-miR396a-5p, cas-miR172a-5p_3ss2_TC3GA9AT, and bra-miR164e-3p_1ss14CT (Tables 2, 4) in the R lines at 2 dpi/6 dpi, upregulate (∼2 to 30-fold) the expression of key target genes such as EL10Ac5g11116 (chaperone protein), EL10Ac3g05363 (ABC transporter), and EL10Ac2g02793 (putative disease resistance RPP13-like protein) respectively in the R lines (vs. S line). These types of genes are well known in plants for their contribution in disease resistance (Cheng et al., 2018; Ul Haq et al., 2019; Devanna et al., 2021). The data presented here might indicate evolutionary advantage in the R lines where genomic/allelic variations lead to the regulation of key miRNAs that contribute to host plant resistance during plant–virus interaction. Presence of specific miRNAs only in the R lines (Tables 1, 3) also suggest that these resistance traits could be purely due to the presence or absence of regulatory elements in the genome. Future comparative genome analysis of the S and R lines (Galewski and Eujayl, 2021) in context to the presence/absence of specific miRNAs will provide better insights on BCTV resistance at the genome level and could possibly contribute to marker development for resistance traits. Besides specific sugar beet miRNA targeting a specific sugar beet gene, it is also evident that in some cases a single miRNA targeting multiple genes and modulating their expression (Tables 1–4). In plants, optimization of cleavage efficiency has been prioritized over maximization of sequence complementarity between miRNAs and their targets during miRNA evolution (Voinnet, 2009; Jones-Rhoades, 2012). This has resulted in the flexibility of miRNAs to be able to target multiple genes. In our study we see both upregulation and/or downregulation of specific ‘R’ genes in the resistant lines following BCTV infection. In plants, LRR and nucleotide binding (NB)-LRR immune receptor proteins are encoded by R genes that are involved in pathogen recognition and initiating resistance responses. Earlier studies have shown that plant miRNAs negatively regulate R genes through the production of trans-acting siRNAs (tasiRNAs) and phased siRNAs (phasiRNAs) which is suppressed during viral/bacterial infection (Zhai et al., 2011; Li et al., 2012). Up-regulation of R genes such as EL10Ac3g05809 (RGA2), EL10Ac2g02633 (RPP13) in the resistant lines might indicate efficient regulation of miRNAs toward resistance against BCTV.

Differentially expressed miRNAs predominantly targeted genes associated mainly with transcription, protein phosphorylation, oxidation-reduction related to gene ontology and carbohydrates, amino acid (tryptophan) metabolism in relation to pathway-enrichment (Figure 4). Many of these cellular processes, pathways, and metabolites are well known in plants for their role in stress tolerance (reviewed in Rojas et al., 2014; Castro-Moretti et al., 2020). Among different metabolites, carbohydrates have been implicated in plant immune response leading to host resistance against pathogens (reviewed in Bolouri Moghaddam and Van den Ende, 2012). Carbohydrate mediated host resistance is achieved mainly through effector-triggered immunity (ETI) and/or pathogen-associated molecular patterns (PAMP)-triggered immunity against diverse type of pathogens including viruses (reviewed in Gouveia et al., 2017). As an example, cellular contents of sucrose and fructose increased in the leaves of Arabidopsis following Tobacco rattle virus (TRV) inoculation at early infection stages (Fernández-Calvino et al., 2014). In our study, we saw an increase in fructose content in the leaves of R lines (vs. S) at 6 dpi and when compared between Line 13 (R) infected versus control plants. Whereas leaf sucrose content increased only in one of the R lines, Line 4, at 6 dpi but decreased in the other R line, Line 13. This may indicate genotypic difference in modulation of specific carbohydrates and no specific trend could be obtained between BCTV infection and leaf sucrose content in the sugar beet R lines (vs. S). The data presented here indicate differential regulation of specific carbohydrates in sugar beets where higher abundance might contribute to resistance, whereas others such as glucose and galactose might result in viral susceptibility. The involvement of other metabolic pathways (Figure 4B) such as tryptophan and sphingolipid metabolism observed here are also in line with host plant resistance against pathogens including viruses (Castro-Moretti et al., 2020).

The data presented here show several fold upregulations of specific miRNAs in the sugar beet R lines (vs. S) including presence of unique miRNAs such as PC-3p-28452_2099 and PC-3p-63465 observed only in the R lines following BCTV infection (Table 5). Further bioinformatics analysis to identify if any of these sugar beet miRNAs could potentially target BCTV genes revealed that a majority of them targeted genes/genome elements associated with capsid protein gene and viral replication (Table 5). In our case we found a 10–12 bp sequence complementarity between sugar beet miRNAs and BCTV target genes/other genome elements. The majority of these sequence complementarities were either at the seed region and/or 5′/3′ regions of the miRNAs. Several attributes of miRNAs that determine the stability of miRNA-mRNA complex and the efficacy of target gene silencing include complementarity at the seed region, 5′ region, specific nucleotide positions, and more recently even beyond the seed region (Chipman and Pasquinelli, 2019; Satish et al., 2019). However, seed region complementarity has been shown to not be an absolute necessity in plants for a miRNA to target a specific gene (Chipman and Pasquinelli, 2019). Critical plant miRNA target sites in viruses vary depending upon the virus species. A >3-fold reduction of the BCTV target genes by sugar beet miRNAs in some cases suggests that they possibly play an important role in down-regulating key virus genes that contribute to pathogenesis during early infection stages. Overall viral load was lower in the R lines in comparison to the susceptible line (Supplementary Figures 5–7) along with the R lines showing significant resistance against the virus (Figure 7). The duration of exposure to the virus did not show any striking difference in disease symptoms in the S line when evaluated after 3 weeks suggest that a small viral load can rapidly replicate inside a susceptible host to produce significant disease symptoms. An interesting observation from our studies show the presence of a miRNA, PC-3p-24046_2377, only in the R lines that targets the virus origin of replication. Thus, we postulate that sequence complementarity between this miRNA and the BCTV origin of replication could inhibit viral replication by spatially competing with the viral DNA polymerase and may contribute to the reduction of viral load. A similar mode of action by miRNAs has been reported in RNA viruses where binding of miRNAs to the viral genome affected viral replication and pathogenesis (reviewed in Trobaugh and Klimstra, 2017). Our observations are in line with the earlier reports where plant derived miRNAs targeted key viral genes thereby contributing to host resistance (reviewed in Liu et al., 2017). In addition, we validated miRNA target gene expression through mRNAseq (Table 5). Future functional characterization through over-expression of some of these miRNAs described above will help in delineating their precise role in BCTV resistance.

Alignment of RNAseq reads to the virus genome showed next generation sequencing (NGS) of sRNAs is an efficient tool to detect virus derived transcripts during early stages of plant infection. Further characterization of these transcripts revealed putative virus derived sncRNAs with the potential to interact with targets in the sugar beet host genome (Tables 6, 7). The recovery of transcripts that were between 31and 41 nt in length that align to the virus genome were of interest. Other pipelines recommended size filters to identify virus-derived siRNAs (vsiRNAs; Leonetti et al., 2021) but would have recovered no vsiRNAs in our dataset. Alternative hypotheses explaining these observations may include degraded mRNA, novel mechanism of processing, less fidelity and regulation for viral transcription or size prior to processing. In any case, it is important to mention that there could be a potential interaction between sncRNAs and host genomes due to the high sequence homology between putative sncRNA core sequences and host target genes. These matches appeared to be what we expected for siRNAs regarding length. The 31-41 nt sequences generated may suggest that the production of these sequences from virus genomes is different from what has been described in other systems. Potential modifications of these transcripts could also be forthcoming. In Entamoeba parasites, specific oligo-adenylated modifications of sRNA populations were observed and represent novel mechanisms of sRNA regulation (Zhang et al., 2020). Reads did span the core sequence on both sides hinting at potential degradation of virus mRNA which was captured in the sRNA fractions. High potential for interaction due to sequence homology, known siRNA match core sequence length, and transcriptional spikes in the genome suggest that these sncRNAs are produced and may regulate the interaction between virus and host (Prasad et al., 2019). When sequences were aligned to a single virus genome, for example, the Severe strain CTS06-021, only 26 sncRNA sequences were found (data not presented here). This showed both the presence and absence of different putative sncRNAs as well as differences in the homology between the sncRNA core sequences and gene targets in the host genome. Further research is needed to investigate BCTV strain specificity and their interaction with the host genome resulting in distinct resistance phenotypes. This should shed more light on how variation within host genomes as well as variation in virus strain genomes modulate host virus interactions thereby contributing to host plant resistance or susceptibility.

A population genomic survey and epidemiological approaches in the field suggest potential for functional divergence among the BCTV strains (CA/Logan, Colorado, Severe and Worland) (Strausbaugh et al., 2017). Our observations highlight a potential mechanism regulating host virus interactions involving putative sncRNAs. The approach allows for complex variables such as virus strain specificity and host genome variability to be investigated. Relative abundance of sncRNAs at different time points and differential regulation of target gene expression in the S vs. R lines suggest importance of such interactions resulting in sugar beet host susceptibility and/or resistance to the virus (Table 7). The R lines were able to maintain higher expression of key genes such as EL10Ac1g01206 (LRR), EL10Ac5g12605 (7-deoxyloganetic acid glucosyltransferase), and EL10Ac6g14074 (transmembrane emp24 domain containing) against the predominant and most effective BCTV sncRNAs. The overall lower relative abundance of virus sncRNAs especially BCTV_sncRNA 4, 18, and 19 in the R lines suggest possible resistance mechanism/s by which these sncRNAs are pacified in the resistant plants (Table 7).