Khushwant SlNGH,Jana JARO?OVá,Jan FOUSEK,CHEN Huan,Jiban Kumar KUNDU

1 Plant Virus and Vector Interactions Group,Division of Crop Protection and Plant Health,Crop Research Institute,Prague 16106,Czech Republic

2 Zhejiang Tianke High-Technology Development Co.,Ltd.,Hangzhou 310012,P.R.China

Abstract High-throughput deep-sequencing technology and bioinformatics analysis of the small RNA (sRNA) population isolated from plants allows universal virus detection and complete virome reconstruction for a given sample.In the present sRNA deep-sequencing analysis of virus-infected wheat samples in the Czech Republic,samples were firstly tested for barley yellow dwarf viruses (BYDVs),wheat streak mosaic virus (WSMV) and wheat dwarf virus (WDV) using ELISA,RT-PCR and PCR.Subsequent sRNA sequencing of these samples yielded more than～60 million single-end 50-bp reads with high confidence for nine field samples of wheat.Overall,16.5% of reads were virus-specific and 83.5% were mapped to the host.More 21-nt reads (～7.7E+06 reads) were found than 24-nt (～6.20E+06 reads) or 22-nt (～4.30E+06 reads) reads.De novo assembly of the high-quality contigs revealed the presence of three earlier reported viruses in the Czech Republic:BYDVs (31.48%),WSMV (24.23%) and WDV (26.66%).We also showed the presence of cereal yellow dwarf virus (14.33%；two species CYDV-RPS and CYDV-RPV (family Luteoviridae/Polerovirus) and wheat yellow dwarf virus (WYDV,3.30%；Luteoviridae).Phylogenetic analysis showed CYDV and WYDV grouped separately from BYDVs.Furthermore,several recombination breakpoints were found among the groups of yellow dwarf viruses (BYDVs,CYDV,and WYDV).Using RNA deep sequencing,we confirmed the presence of the three known viruses (BYDVs,WSMV,and WDV) and the first record of two species of CYDV and WYDV in wheat in the Czech Republic.

Keywords:NGS,sRNA,BYDVs,CYDV,WSMV,WYDV,WDV

1.lntroduction

Barley yellow dwarf viruses (BYDVs),wheat dwarf virus(WDV) and wheat streak mosaic virus (WSMV) are among the major viral pathogens causing serious diseases in cereal crops (Kunduet al.2009b；Jaro?ováet al.2013；Singhet al.2018).At least six viral pathogens,including BYDVs(species PAS,PAV,and MAV),WDV (wheat strain and barley strain),WSMV,ryegrass mosaic virus (RgMV),brome mosaic virus (BrMV) and lolium latent virus (LoLV),have been detected in cereal crops in the Czech Republic (Drabet al.2014).The genomes of these viruses can be highly diverse within a single species (Robinson and Murray 2013),and recombination among the isolates is frequent (Singhet al.2018).Thus,a high-throughput multiple diagnostic tool has been needed to discriminate and quantify viruses and enhance virus detection in cereals and their wild relatives,which are often natural reservoirs for viruses and play a significant role in virus dispersal and epidemic.

Next-generation sequencing (NGS),using deep RNA sequencing (Adamset al.2009；Kreuzeet al.2009) has allowed great advances not only for universal virus detection and genome reconstruction but also for complete virome reconstruction (Pooggin 2018；Turcoet al.2018).RNA sequencing have been successfully applied for analyzing the virome (Massartet al.2019),detecting viruses in cereal crops fromLuteoviridaefamily such as maize yellow mosaic virus (MaYMV；Chenet al.2016) and wheat leaf yellowing-associated virus (WLYaV；Zhanget al.2017),and characterizing a novel virus species or genus (Pecmanet al.2017；Yanget al.2017；Massartet al.2019).In particular,RNA-sequencing (RNA-Seq) of small RNAs(sRNAs),the products of RNA interference (RNAi),an antiviral defense system that targets dsRNA in plants (Ding and Voinnet 2007),allows all potential viral populations and species in a given sample to be detected in a single assay without any previous knowledge of the viral genomes(Roossincket al.2015) and is very effective for virome identification (Kreuzeet al.2009；Massartet al.2017) in crops (Kashifet al.2012；Adamset al.2013；Wylieet al.2014；Ibabaet al.2017) and wild plants (Biet al.2012).Most plant viruses have RNA genomes,and DNA genome viruses transcribe the genome to yield mRNA during virus replication；hence,incorporatingde novoassembly makes RNA-Seq a high-throughput technique to identify both known and novel RNA and DNA viruses (Barbaet al.2014).It is also 10 times more sensitive than RT-PCR (Santala and Valkonen 2018).In addition,sRNA sequencing has been used to detect many cereal-infecting viruses from different virus genera such asMachlomovirus,Potyvirus(Adamset al.2013),Cytorhabdovirus(Yanget al.2017),Mastrevirus(Chenet al.2015) andPolerovirus(Chenet al.2016；Zhanget al.2017).

In this report,we thus used sRNA-Seq to analyze the virome in wheat samples from various fields in the Czech Republic.Besides the commonly distributed viruses such as BYDVs (PAS,PAV),WDV and WSMV,we detected cereal yellow dwarf virus (CYDV and CYDV-RPS and CYDV-RPV)and wheat yellow dwarf virus (WYDV),the first record of these viruses in the Czech Republic.We also analyzed the evolutionary trajectory of yellow dwarf viruses (YDVs)isolates including BYDVs (PAV,PAS,MAV,and GAV),CYDV and WYDV.Apart from revealing their genetic diversity,the analysis of genetic recombination showed several recombination events in CYDV and WYDV isolates from the Czech Republic.

2.Materials and methods

2.1.Sample collection

In 2016,symptoms of infection by wheat streak mosaic virus (WSMV),barley yellow dwarf viruses (BYDVs) and wheat dwarf virus (WDV) were observed on plants ofTriticum aestivumin fields in the Czech Republic.Leave samples with visible symptoms of aforementioned viruses were collected from various fields in the Czech Republic.Collected samples were stored at -80°C until further use.

2.2.Detection of viruses using ELlSA and PCR

For serological detection of viruses including BYDVs,WDV and WSMV,polyclonal antibodies (SEDIAG；http://www.sediag.com/) were used in DAS-ELISA to detect the viruses.Samples for ELISA were prepared by grinding 1 g leaf tissue at a 1:20 ratio in phosphate-buffered saline,pH 7.4,with 2% polyvinylpyrrolidone and 0.2% egg albumin.Optical density at 405 nm of samples in microplates was read using a MR 5000 Dynatech Reader (Dynex Technologies,Chantilly,VA,USA).

For PCR amplification using primers of the respective coat protein of BYDVs (Kunduet al.2009b),WDV (Kunduet al.2009a) and WSMV (Chalupnikováet al.2017).For CYDV and WYDV primers described by Deb and Anderson(2008) and Zhanget al.(2017) were used,respectively.Total RNA was extracted from wheat leave samples using a Spectrum Plant Total RNA Kit (Sigma Aldrich,St.Louis,MO,USA) as previously described (Singh and Kundu 2017).Complementary DNA was synthesized using a Reverse Transcription System (Promega,Madison,WI,USA).A reaction mixture composed of 1 μg total RNA and 0.5 μg random primers was incubated at 70°C for 5 min and chilled on ice for 2 min.Then,the following components were added in the given order:5× reaction buffer,20 U mL-1of RNasin RNase Inhibitor (Promega),10 mmol L-1dNTP mix and 200 U mL-1of M-MLV reverse transcriptase enzyme.The mixture was incubated at 37°C for 1 h,then the reaction was stopped by heating at 70°C for 10 min,then chilling on ice.

2.3.RNA extraction,library preparation and sRNA sequencing

Total RNA was extracted from leaves using a Spectrum Plant Total RNA Kit (Sigma Aldrich,St.Louis,MO,USA)as described above.sRNA-Seq libraries were prepared using NEBNext Small RNA Library Prep Set for Illumina and NEBNext Multiplex Oligos for Illumina (New England Biolabs,Ipswich,MA,USA).The library was quantified using a KAPA Library Quantification Kit (Roche,Basel,Switzerland),and an Agilent Bioanalyzer 2100 (Agilent Technologies,Santa Clara,CA,USA) was used for quality control.To identify potential viruses causing the symptoms,we used a NovaSeq6000 System (Illumina,San Diego,CA,USA) for deep sequencing of the sRNAs.High-quality raw reads were received in a FASTQ format.Sequencing was performed commercially by SEQme s.r.o.,(Dobris,Czech Republic).

2.4.Data analysis using bioinformatics

Raw reads were cleaned of potential contamination.Adapter sequences were removed using Cutadapt v1.9.1(https://cutadapt.readthedocs.io/en/stable/),and reads less than 17 nucleotides (nt) or longer than 30 nt were filtered out.Quality control of raw sequence data was checked using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/).The remaining high quality reads were aligned with sequences in the Plant rRNA Database(http://www.plantrdnadatabase.com/) using Bowtie version 1.2.2,allowing up to three mismatches (http://bowtie-bio.sourceforge.net/index.shtml).Unaligned reads were used for virus discovery using VirusDetect version 1.7(Zhenget al.2017).The non-redundant (95%) plant virus database was downloaded from GenBank and used as the reference.The Burrows-Wheeler Aligner (BWA)version 0.7.12-r1039 was used to align sRNA reads to the reference virus sequences,host sequences (T.aestivum),or assembled virus contigs (Li and Durbin 2009).For reference-guided assembly,SAMtools version 1.9 (http://samtools.sourceforge.net/) was used to process BWA alignments and generate per-position alignment information,which was used to guide the construction of virus contigs.De novoassembly of viral sRNAs was performed using Velvet (v1.2.09) (https://www.ebi.ac.uk/).Resulting contigs were mapped to the host genome,and unmapped contigs were mapped using standalone blast version 2.8.1 (https://blast.ncbi.nlm.nih.gov/) to the virus reference nucleotide(BLASTN) and protein sequence databases (BLASTX).The sequences were edited using Sequencer 4.8 (Singh and Kundu 2017) and submitted to NCBI,accession numbers:CYDV (MN179491) and WYDY (MN179492).

2.5.Multiple sequence alignment and phylogenetic analysis

Multiple sequence alignment of nucleotide sequences of YDVs was performed as described by Singh and Kundu(2017).Phylogenetic trees were constructed using Clustal X(version 2.1) (www.clustal.org/clustal2/) with default parameters and the neighbor-joining method in MEGA 7.0(www.megasoftware.net/).The reliability of branches was inferred from bootstrap analysis of 1 000 replicates.The final phylogenetic tree was then edited using ITOL (http://itol.embl.de/).

2.6.Recombination event analysis

Recombinant sequence analysis and localization of recombination breakpoints in the genome of WYDV and CYDV nucleotide sequences from wheat samples from the Czech Republic was done using RDP4 Software as described earlier (Singh and Kundu 2017) and various algorithms (RDP,GENECONV,Chimaera,MaxChi,Bootscan,SiScan,3seq,LARD).For reference,WYDVGPV-UK (WYDV reference) and CYDV-RPV-MI-USA (CYDV reference) were used.Recombination events withP＜0.01 that were detected by at least three of the eight algorithms in the RDP Software are shown.Numbers of events withP＞0.01 are given in parentheses.Nucleotide and amino acid divergence was estimated using MEGA 7.

3.Results

3.1.Viruses detected

The viruses detected in symptomatic wheat samples from various locations in the Czech Republic using ELISA and RT-PCR for BYDVs,WDV and WSMV are listed in Table 1.

3.2.Composition of sRNA populations in wheat samples

In the total sRNA profile from nine wheat samples determined using high-throughput sequencing,we obtained～63 million single-end 50-bp reads with high confidence.The number of reads varied from～5 to～10.4 million reads(average～7.3 million reads per sample).Approximately 41 million (～61%) clean reads were processed for the analysis (Table 2).When clean reads were mapped to the genome ofT.aestivum(bread wheat),83.5% of the reads(～52.5 million reads) were found to be host-specific,and 16.5% were virus-specific (～10.6 million reads) (Fig.1-A).The length of the clean reads ranged from 18 to 30 nt(Fig.1-B)；21-nt reads were the most abundant (～7.7E+06 reads),followed by 24-nt (～6.20E+06 reads) and 22-nt(～4.30E+06 reads).

3.3.Virome in the samples

De novoassembly of the reads that were not mapped to the wheat genome using the Velvet program generated4 298 contigs.In the comparison and annotation of the assembled contigs with the virus reference database using BLAST (Appendices A and B),more than 80% contigs had high identity with cereal viruses from (i) familyLuteoviridaeincluding barley yellow dwarf viruses(Luteovirus；BYDVPAS,BYDV-PAV,BYDV-MAV,and BYDV-GAV),cereal yellow dwarf virus-RPS (Polerovirus；CYDV-RPS) and wheat yellow dwarf virus-GPV (WDV)；(ii) familyGeminiviridae:WDV (Mastrevirus)；(iii) familyPotyviridae:wheat streak mosaic virus (WSMV；Tritimovirus) (Table 3).

Table 1 Wheat (Triticum aestivum) samples collected from the Czech Republic for small RNA (sRNA) deep sequencing and ELISA and PCR results for BYDV,WDV and WSMV1)

Table 2 Small RNA (sRNA) profile of wheat (Triticum aestivum) samples from the Czech Republic

Fig.1 Small RNA (sRNA) profile of wheat samples from the Czech Republic.A,percentage of host-specific and virus-specific reads.B,length distribution of total sRNAs.The 21-nt sRNAs were more abundant than the 24-and 22-nt sRNAs.

3.4.Phylogenetic analysis of yellow dwarf viruses

In the phylogenetic tree using the contigs for CYDV and WYDV and known yellow dwarf viruses sequences,three clusters formed:(i) cluster 1:BYDV-PAV,BYDV-MAV,BYDV-PAS and BYDV-GAV from the United States,Europe(Sweden),and Asia (Pakistan,China,Japan,and Iran)；(ii)cluster 2:WYDV-Czech and WYDV-GPV from China and the United Kingdom；(iii) cluster 3:CYDV-RPS and CYDVRPV from the United States and CYDV-Czech (Fig.2).Thedetection of CYDV and MYDV were further confirmed by RT-PCR as described earlier by Deb and Anderson (2008)and Zhanget al.(2017),respectively.

Table 3 Viruses identified from the viromes in wheat (Triticum aestivum) samples from the Czech Republic using small RNA(sRNA) deep sequencing1)

3.5.Recombinant analysis revealed possible variations among WYDV and CYDV

The predicted recombination events from the recombinant event analysis using different yellow dwarf virus isolates including WYDV and CYDV are shown in Fig.3.Recombination breakpoints were distributed along the entire genome.Although we found many recombination breakpoints in the boundaries of ORFs,we also found many putative points within coding sequences.An isolate of WYDV-GPV from the United Kingdom underwent recombination with a WYDV-Czech isolate at nucleotide positions 450-1 110 in ORF1 (encoding the RdRp),while it recombined at the 3′ end with isolate CYDV-RPV from the United States at nt 4 200-4 750 in ORF3 and ORF5 (encoding RTP).Isolate WYDV-Czech showed a recombination breakpoint at nt 1 100-1 700 (ORF2) and nt 3 350-4 320 (ORF3 and ORF4) with the WYDV-GPV-UK isolate (Fig.3).On the other hand,isolate CYDV-RPV from the United States recombined with isolate CYDV-Czech at the 5′ end at nt 680-2 100 (ORF1 and ORF2,encoding RdRp) and at the 3′ end with WYDV-GPV-UK at nt 4 600-5 100 (ORF5).CYDV-Czech showed a recombination event at a central position with CYDV-RPV from the UnitedStates at nt 2 320-3 710 (ORF2,ORF3a and ORF4),and CYDV-RPS isolate RPV-Mex-1 from the United States at nt 2 200-2 800 (ORF2).CYDV-RPS isolate RPV-Mex-1 from the United States underwent recombination events with CYDV-RPV isolate from the United States at nt 330-4 700(ORF3a,ORF3,and ORF4) (Fig.3).

Fig.2 Phylogenetic tree reconstructed using the neighbor-joining method based on alignment of whole genome sequences from yellow dwarf viruses (YDVs).Sequences used:BYDV-PAV-064-USA (EF521850.1)；BYDV-PAS-KSPAS-1-USA (KY593456.1)；BYDV-PAV-0109-USA (EF521828.1)；BYDV-PAS-KSPAS-2-USA (KY593457.1)；BYDV-PAV-KS-SHKR-USA (KU170668.1)；BYDV-PAV-KS-PAV-USA (KY593458.1)；BYDV-PAV-M14-Pakistan (HE985229.1)；BYDV-PAV-Kerman-Iran (KP771878.1)；BYDVPAV-052-Sweden (EF521841.1)；BYDV-PAV-0100-USA (EF521843.1)；BYDV-PAV-Japan (D85783.1)；BYDV-PAV-PAV-CN-China(AY855920.1)；BYDV-GAV-YL4-China (KF523380.1)；BYDV-MAV-USA (NC_003680.1)；WYDV-GPV-China (NC_012931.1)；WYDV-GPV-UK (FM865413.1)；CYDV-RPV-005RPV-USA (EF521827.1)；CYDV-RPV-MI-USA (KY553235.1)；CYDV-RPS-RPVMex-1 (AF235168.2)；CYDV-RPV-046-USA (EF521839.1)；CYDV-RPV-USA (L25299.1)；CYDV-RPV-010-USA (EF521830.1)；CYDV-RPV-092-USA (EF521848.1)；and CYDV-RPS-MI-USA (KY623680.1).The tree was constructed using 1 000 bootstrap replicates.Nodes with bootstrap values ＜60% were collapsed.An isolate of oat necrotic mottle virus (ONMV；accession number AF454461) was used as an outgroup sequence to root the tree.

4.Discussion

ELISA is routinely used to monitor cereal crops for specific viruses,usually the most common ones；less severe,less frequent or unknown viruses will be missed since more than 100 viruses from various genera are reported to infect cereal crops (Lapierre and Signoret 2004).Several PCR-and Sanger-sequencing-based methods have also been used to detect and quantify (Kunduet al.2009a,b；Jaro?ová and Kundu 2010；Boonhamet al.2014；Drábet al.2014；Jaro?ováet al.2018) viruses in cereal crops and volunteer grasses,but PCR methods require primers/probes designed for a known sequence of a particular viral genome (Jaro?ováet al.2018).

In present study of the wheat virome using deep sequencing of sRNAs and targeting 17-30 nt,the length distribution of the sRNAs showed that 21-nt sRNAs were the most abundant,followed by 24-and 22-nt reads,in accordance with reports elsewhere (Donaireet al.2009；Zhanget al.2017) (Fig.1).The assembled contigs associated with identified viruses in our study were sufficient to cover most regions of the target viruses including BYDVs,WDV and WSMV.Previous monitoring of various cereal hosts in the Czech Republic has identified BYDV species PAS,PAV and MAV (Kunduet al.2009b；Jaro?ováet al.2013；Beoniet al.2016) belonging to the genusLuteovirus,but no viruses belonging to the genusPolerovirus,although the high prevalence of aphid vectors,which efficiently transmit many virus species from the familyLuteoviridae,were found in a wheat field (Jaro?ováet al.2013).Ourde novoassembly and analysis of the viromeidentified not only BYDVs,WSMV,and WDV,but also two viruses new to the Czech Republic,CYDV and WYDV (former BYDV-GPV；Zhanget al.2009),both belonging toPolerovirus(Kruegeret al.2013).

Fig.3 Recombination analysis of whole genome sequences of yellow dwarf viruses (YDVs) using various algorithms in the RDP Software package.Wheat yellow dwarf virus (WYDV)-GPV-UK (WYDV reference) and cereal yellow dwarf virus (CYDV)-RPV-MIUSA (CYDV reference) were used as references.The order of the designation for the events is as follows:isolate name,accession number,algorithm used (R,RDP；G,GENECONV；C,Chimaera；M,MaxChi；B,Bootscan；SS,SiScan；3s,3seq；L,LARD),P-value and nucleotide positions.Only recombination events with P＜0.01 detected by at least three algorithms are shown.Numbers of events with P＞0.01 are given in parentheses.Recombination events were found using RDP,GENECONV,Bootscan,Chimaera,SiScan,and MaxChi.

These two poleroviruses,reported from many wheatgrowing regions worldwide (Hadiet al.2012；Milgateet al.2019；Wallset al.2019),are transmitted byRhopalosiphumpadiandSchizaphis graminum,very common aphids in the Czech Republic.Our sequence comparison and phylogenetic analysis of CYDV and WYDV indicated that they are distinct from the commonly detected BYDVs (Fig.2)and genetically closely related to each other,as reported by Boulila (2011).Several recombination events among CYDV and WYDV isolates were also found (Fig.3).In yellow dwarf viruses,most recombination breakpoints occur at the boundaries of ORFs and within coding sequences (Pagán and Holmes 2010).In WSMV genome sequences,great diversity has been found within a single species from various hosts,indicating recombination among the isolates (Singhet al.2018).Therefore,it would be intriguing to analyze the mutations in the sequences of CYDV and WYDV from various hosts.

Luteoviruses and poleroviruses are among the most economically important causal agents of disease in cereal crops,causing significant yield losses and reductions in grain quality (McKirdy and Jones 2002).The viruses share common aphid vectors with distinct transmission efficiency(Gray and Gildow 2003).As a result of climate change,the abundance and distribution of aphid species (Tr?bickiet al.2016) are also changing,thus influencing the incidence and distribution of viruses in cereal crops (Trebickiet al.2015).Therefore,further monitoring studies using ELISA and following RNA sequencing will help to determine the impact of climate change to the dynamics and diversity of viruses in cereal crops.

5.Conclusion

The sRNA deep sequencing used here significantly improved virus diagnostics for wheat.Thede novovirome reconstruction approach allowed us to identify cereal yellow dwarf virus and wheat yellow dwarf virus (bothPolerovirus)(for the first time in the Czech wheat crops.The phylogenetic analysis showed that the newly detected CYDV and WYDV clustered separately from BYDVs.Further field monitoring of CYDV and WYDV will be aimed at elucidating the dynamics,diversity,incidence and spread of the viruses in cereal crops and identifying and confirming the virus vectors.We anticipate that findings from this study will lead to a better understanding of the incidence and distribution of different wheat viruses in the Czech Republic and thus the design of suitable control measures.

Acknowledgements

We thank Ms.Michaela Bro?enská (Plant Virus and Vector Interactions Group,Division of Crop Protection and Plant Health,Crop Research Institute,Czech Republic) for her excellent technical assistance.We also thank Dr.Beth E.Hazen (Willows End Scientific Editing and Writing Cortland,NY,USA) for editing the English of the manuscript.This work was supported by a grant from the Technology Agency of the Czech Republic (TF02000056).

Appendicesassociated with this paper can be available on http://www.ChinaAgriSci.com/V2/En/appendix.htm