During field surveys in April 2016, 12 wheat samples exhibiting flag leaf yellowing and little dwarfing were collected from Jinan, Shandong Province of China. Whole plants were put into plastic bags and their roots kept moist during transport to the laboratory. One hundred milligrams of fresh leaves were cut from the plants to extract total RNA; other leaves were stored at −70°C. Healthy wheat plants growing in an insect-proofed greenhouse served as controls.

Total RNA of wheat leaves was extracted using TRIzol reagent according to the manufacturer's instruction (Invitrogen, USA). The quality of RNA was tested using a Nanodrop 2,000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and RNA integrity was verified by gel electrophoresis.

BYDV-GAV, -PAV and WYDV-GPV were detected by RT-PCR using specific primers detailed in Table S1, which also gives the sequences of all primers used in this study. B/W-YDV-infected plants that were used as positive controls were growing in the growth chamber (16 h light at 20°C/8 h dark at 18°C). Reverse transcription was performed in a total volume of 20 μL with 1 μg of RNA, 2 μL of reverse (R) primer (10 μM), 2 μL of dNTP Mix (each 2.5 mM; TaKaRa, China), 4 μL of 5× M-MLV buffer (Promega, USA), 1 μL of M-MLV reverse transcriptase (200 U/μL; Promega), 0.5 μL of Recombinant RNase Inhibitor (40 U/μL; TaKaRa) and DEPC water. The PCR contained 2 μL of cDNA, 2.5 μL of 10× buffer (Mg²⁺, 15 mM, TaKaRa), 2 μL of dNTP Mix (each 2.5 mM; TaKaRa), 0.5 μL of forward (F) and reverse (R) primer (10 μM), 0.2 μL of rTaq (5 U/μL; TaKaRa) and 17.3 μL of ddH₂O. The PCR reaction was performed using thermal cycler (Bio RAD, Hercules, CA, USA) as follows: denaturation at 94°C for 3 min; 35 cycles at 94°C for 30 s, 58°C for 45 s and 72°C for 80 s; final extension at 72°C for 10 min (Zhao et al., 2010). PCR products were electrophoresed in 1% agarose gel and stained with ethidium bromide (EB).

The quality of extracted RNA was tested using the Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific) and Bio Analyzer 2100 (Agilent Technology, Santa Clara, CA, USA). An RNA sample with an RIN (RNA integrity number) value ≥6 was considered as high enough quality and used in the next steps. The small RNAs (sRNAs) were separated using PEG8000 and ligated to a 3′ adaptor. The 36~44-nt (ligated sRNAs) bands were purified using PAGE in a 15% denatured gel and linked to 5′ adaptor. The first-strand cDNA was synthesized using the pair-linked sRNAs, and 16 cycles of PCR amplification were performed. The 140~160-bp products, namely the sRNAs library, were purified by 3.5% agarose gel electrophoresis. The library was quantified by ECO (Illumina, San Diego, CA, USA) and submitted for sequencing on the Illumina 2,500 platform (Illumina).

Raw reads from the Illumina platform were processed to trim adaptor sequences and low-quality reads. Unique sequences were then generated as clean reads by collapsing the identical sequences. The abundance of sRNAs was determined using the software Bowtie2 and default parameters (Langmead and Salzberg, 2012). Then the clean reads were assembled using the Velvet program with a k-mer of 17 as the minimal overlapping length required for joining sRNAs into larger contigs (Zerbino and Birney, 2008). Assembled contigs were screened against the GenBank nucleotide collection (nt) and non-redundant protein sequences (nr) databases using a BlastN and BlastX search using standard parameters respectively (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

Contigs and gaps were confirmed and filled by RT-PCR using specific primer pairs 1, 2, 3, 4, 1-2, 2-3, and 3-4 F/R (Table S1) designed according to the consensus sequences of the assembled viral contigs with more than 70-nt overlapping. The sugarcane yellow leaf virus (ScYLV, AF157029.1) genome sequence was used for the alignment and positioning of the contigs. The terminal sequences of the viral genome were obtained using 5′ and 3′ rapid amplification of cDNA ends (RACE) kits according to the manufacture (Invitrogen). The PCR products were collected using Wizard SV Gel and PCR Clean-Up System (Promega). The purified PCR products were then cloned into the pEASY-T5 vector (TransGen Biotech, China) and used to transformed Trans-T1 competent cells (TransGen Biotech) following the manufacturers' instructions. The clones harboring the transformed vector were identified by PCR and Sanger sequencing (Sangon Biotech (Shanghai) Co., China). The results of sequencing were assembled using DNAMAN (version 6) program (Lynnon Biosoft, San Ramon, CA, USA) with more than 70-nt overlapping regions to form the full-length viral genome.

After obtaining the whole viral genome, open reading frames (ORFs) were predicted using SnapGene software (Ian, 2004). The distribution and coverage of vsiRNA was determined using the software Bowtie2 under default parameters, and the results were exported to Excel (Microsoft, Redmond, WA, USA) for further analysis. Contigs were mapped to the WLYaV genome using the CLC Genomics Workbench (Qiagen, Valencia, CA, USA).

The identities of nucleotide sequences of the whole genome, untranslated regions and deduced amino acid sequences of seven ORFs of WLYaV with other viruses belonging to the family Luteoviridae (Table S2) that were retrieved from NCBI (http://www.ncbi.nlm.nih.gov/) were calculated by the Needle program (Liu et al., 2009). Sequences were aligned with the ClustalW method, and phylogenetic trees were constructed by the neighbor joining method using MEGA 6 software (Tamura et al., 2013). The reliability of each branch was evaluated with bootstrap (1000 repeats).

The whole nucleotide sequence of WLYaV was divided into two overlapping parts, A (nt 1–3 059) and B (nt 3 035–5 772), a NcoI restriction enzyme site was within the overlapping region (nt 3 045–3 050). A and B were amplified by RT-PCR with specific primer pairs (A-Stu-F/A-Nco-R, B-Nco-F/B-Sal-R, respectively), adding restriction enzyme sites StuI and SalI, respectively. The amplification products of fragments A and B were purified and cloned into the pEASY-T5 vector and sequenced as above. pCB301 vector (kindly supplied by Prof. Xiaorong Tao, Nanjing Agricultural University), a binary vector with 2 × 35S promoter, ribozyme and NOS terminator (Shen et al., 2014; Wang et al., 2015), and plasmids containing fragment A and B without mutations were digested with restriction enzymes and ligated using T4 ligase (Promega, USA) to produce pCB301-WLYaV. Plasmid pCB301 (as a negative control) and pCB301-WLYaV were transferred to separate suspensions of Agrobacterium tumefaciens strain GV3101 (Wang et al., 2015). Then the abaxial surface of 2-week-old Nicotiana benthamiana leaves (3~4 leaves) was infiltrated with the bacterial cultures (Chen et al., 2016). At 14 days post inoculation (dpi), systemic virus dissemination was evaluated in non-inoculated leaves by RT-PCR using specific primer pair CP-F/R (Table S1).

In 12 wheat samples from Jinan, China that had yellow dwarf symptoms, three known B/W-YDVs (-GAV, -PAV and -GPV) in China were detected. Eight samples were infected with one of the BYDVs, and none had a mixed infection (Table S3). None of the above three viruses were detected in sample JN-U3, which exhibited dull yellowing on the leaf margins and little dwarfing (Figures 1B,C) comparing with the healthy sample (Figure 1A), while bright yellowing was found on samples infected by B/W-YDVs. Therefore, sample JN-U3 was used for the next experiment.

The sRNA library from sample JN-U3 was constructed and sequenced using Illumina Hiseq 2,500 platform. In total, 25,605,387 clean reads that had passed quality control were produced after removing adaptor sequences and low-quality reads. The length of the clean reads ranged from 18 to 25 nt (Figure 2); 21-nt reads were the most abundant, followed by 24 and 22-nt reads. After assembly using the Velvet program, 8,798 contigs were generated.

After blast analysis of the assembled contigs, four contigs of 1,151, 701, 718, and 1,714 nucleotides, respectively, had high identities with several viruses in the family Luteoviridae, and may be associated with the yellowing symptom. The higher identities were observed with ScYLV, with 83% identity on 99% coverage of the sequence for contig 3, for example. Four contigs were aligned to different positions of the ScYLV genome, and their relative positions and gaps were predicted. RT-PCR and RACE-PCR were performed to confirm the contig sequences, fill gaps and obtain the terminal sequences (Figure S1). The RACE-PCR was repeated three times with the similar results. Finally, the assembled genome of this virus was constituted by 5,772 nt, and deposited in GenBank as accession KY605226. The other 11 samples from Jinan were negative for WLYaV (Table S3) using RT-PCR and primer CP-F/R (Table S1).

Because the newly discovered virus was found in a wheat sample that had dull yellowing and little dwarfing symptoms, the virus has been temporarily named wheat leaf yellowing-associated virus (WLYaV). Six ORFs (ORF0, ORF1, ORF2, ORF3, ORF4, ORF5) were predicted by SnapGene software, and 5′ and 3′ UTRs, one intergenic-UTR (I-UTR) and ORF3a were also found in the viral genome (Figure 3A), suggesting the genome organization is similar to that of viruses in the genus Polerovirus (Miller et al., 1995). The 5′ UTR consists of 58 nt starting with ACUAAA, while the conserved sequence that start the 5′ end of poleroviruses is ACAAAA (Guilley et al., 1994; Moonan et al., 2000; Mo et al., 2010). ORF0 (nt 59~877) is the first ORF encoding a 31.4 kD protein (P0), which is a putative RNA-silencing suppressor (RSS) (Pfeffer et al., 2002, 2017; Kozlowska-Makulska et al., 2010). A putative F-box like motif (IPIIL) conserved among poleroviruses (LPxxL/I) has been found (Figure S2), which is required for RSS function (Pazhouhandeh et al., 2006). ORF1 (nt 228~2,099) encodes a 69.1 kD putative polyprotein (P1) with a proteinase motif that can generate three proteins including the viral genome-linked protein (VPg) (Figure S2) (Miller et al., 2002; Nickel et al., 2008). ORF2 (nt 1,658~3,364) and ORF1 compose a putative fusion protein of RdRp (P1-P2 fusion) by a -1 ribosomal frameshift at nt 1,658, which is regulated by a GGGAAAC sequence at nt 1,658~1,664 within the overlapping region of ORF1 and ORF2 and is involved in the conserved GDD amino acid core (Figure S2) (Prüfer et al., 1992; Nixon et al., 2002). ORF2 is followed by a 188 nt I-UTR (nt 3,365~3,552). This sequence contains the ORF3a starting with a putative CUG initiation codon (nt 3,435~3,569; 4.91 kD), which encodes a protein that involved in the viral long-distance movement (Smirnova et al., 2015). ORF3a, ORF3~ORF5 might be expressed through the 3′ co-terminal sgRNAs, yielding the putative capsid, CP-RTD (RTP) and movement proteins (Hwang et al., 2013). ORF3 (nt 3,553~4,143) encodes a 21.7 kD protein of putative coat protein (P3). ORF4 (nt 3,584~4,036) present within ORF3 and encodes a 17.0 kD putative movement protein (P4). ORF5 (nt 4,144~5,577) and ORF3 encode a fusion protein (P3–P5 fusion), a putative readthrough protein that may be involved in virus transmission (Brault et al., 1995). The 3′ UTR consists of 195 nt (nt 5,578~5,772) without a polyA tail.

After the complete viral genome was obtained, total contigs were remapped to the viral genome (Figure 3B). This time, we obtained 11 mapped contigs, which nearly covered the entire genome, while initially only 4 of the 11 contigs had homologies with Luteoviridae family. This difference may be due to the low identities of the seven unrecognized contigs with other known viruses.

To identify the hotspots of viral siRNAs (vsiRNAs), we analyzed the genome region targeted by RNAi by mapping the sRNAs on the viral genome (Figure 3C). A total of 1 301,797 reads (5.25% of total sRNAs reads) were mapped. They were well-scattered on the viral genome and anti-genome strands except for several peaks at ORF0 (RSS), ORF3 (CP) and ORF5 (RTD). The higher number of siRNAs corresponding to CP and RTD would correspond to a higher representativity of the RNA targets that correspond to both genomic and subgenomic RNA species.

To comprehensively understand vsiRNAs, we analyzed the size, polar distribution and base bias of vsiRNAs (Figure 4). The vsiRNAs were mostly 19~24-nt long; the 22- and 21-nt sRNA species had the highest percentages, together making up more than 90% of the total (Figure 4A). We also found that the number of vsiRNAs reads from the anti-genome strand was slightly higher than from the genome strand. It was also observed that the vsiRNAs were evenly distributed in the positive and negative strands of the genome, indicating that the vsiRNAs were formed from the replication intermediate. In addition, a bias toward base A and U especially in the anti-genome was observed (52.1% in the genome, 67.3% in the anti-genome) (Figure 4B).

The full-length genome sequence and the predicted UTR and ORFs of WLYaV were compared with those of 26 other viruses in the family Luteoviridae (Table 1). The nucleotide identities ranged from 42.3 to 64.9% for the full-length genome, 17.8 to 50.0% for the 5′UTR, 33.0–60.3% for the I-UTR, and 11.1–48.4% for the 3′UTR, respectively. Moreover, the deduced amino acid sequence identities were 2.1–25.4% for ORF0, 4.5–48.9% for ORF1, 10.8–62.3% for ORF1-2 fusion, 31.1–72.7% for ORF3a, 29.7–86.2% for ORF3, 20.2–74.0% for ORF4, 29.9–66.0% for ORF3-5 fusion, respectively. It is worth emphasizing that the CP produced by ORF3 shared the highest identities among the viruses, and ORF0 had the lowest identities (Table 1). WLYaV had the closest relationship with ScYLV (64.9% for full length nucleotide and 86.2% for deduced amino acid sequences of ORF3). When the full-length genome and ORF3 deduced amino acid sequences of 10 representative ScYLV isolates (Elsayed et al., 2017) were compared with WLYaV (Table S4), the identities were 63.7~65.3% and 85.7~86.2%, respectively, below the species demarcation thresholds (90%) in the family Luteoviridae (Domier, 2012; Simmonds et al., 2017).

The neighbor-joining phylogenetic analyses based on the nucleotide sequences of the full-length genome and deduced amino acid sequences of RdRp were performed with a bootstrap of 1,000 replications (Figure 5), which are widely used as important criteria for assigning virus species to a genus in the family Luteoviridae. Besides, the phylogenetic relationship of CP sequences was analyzed due to a possible recombinant site had been identified between RdRp and CP of ScYLV (Moonan et al., 2000; Smith et al., 2000). The three phylogenetic trees had similar topologies and showed that WLYaV was most closely related to ScYLV, genus Polerovirus. It's worth noted that in the CP tree, WLYaV, and ScYLV clustered with the members of genus Luteovirus, in contrast, two known luteoviruses, soybean dwarf virus (SbDV) and bean leafroll virus (BLRV) clustered with poleroviruses (Figure 5C). Therefore, a similar recombination event has probably occurred in the ancestral WLYaV genome. It could be speculated that this event may have occurred in the ancestral WLYaV and ScYLV, as similar to a drawing scenario by Domier et al. (2002) with the two luteoviruses, SbDV and BLRV. Furthermore, when the phylogenetic trees were constructed using the 10 representative isolates of ScYLV (Figure S3), WLYaV and ScYLVs were on different branches, giving the same result and indicating WLYaV and ScYLV were different. When the identities and phylogenetic analyses are considered together, this new virus could be defined as a new member of genus Polerovirus, family Luteoviridae according to criterion of the International Committee on Taxonomy of Viruses (ICTV).

The full-length cDNA clone of WLYaV was constructed using vector pCB301, namely pCB301-WLYaV (Figure 6A), and confirmed by colony PCR using primers A-Stu-F and B-Sal-R (Figure 6B) and Sanger sequencing. The positive clone was used to agroinoculate 2-week-old leaves of N. benthamiana plants. Typical systemic symptom of leaf yellowing was observed 2 weeks later and more obvious 4 weeks later (Figure 6C). Leaves were confirmed as systemically infected using RT-PCR and subsequent sequencing of the amplification product (Figure 6D). This test was repeated three times; each time about 60–70% plants were systemically infected.

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