The plant materials used in this study were sampled at an experimental field in the HOKUREN Agriculture Research Institute, Naganuma, Hokkaido, Japan (43.3°N) in May and June of 2018 and 2019 (a total of 14 wheat samples or pools; Table 1). The leaves of wheat plants showing yellow mosaic and asymptomatic leaves, including from the cultivars “Kitahonami” (KTH) in 2018−2019, and “Kitanokaori” (KTN) and “Yumechikara” (YM, a WYMV resistant cultivar; Kojima et al., 2015) in 2019, were collected (Supplementary Figure S1 and data not shown). The wheat plants collected in 2019 were grown on plots with or without the fungicide fluazinam (Flu; Ishihara Sangyo Kaisha, LTD.) soil treatment. Wheat (KTH) and rye (Secale cereale cv. Fuyumidori) samples (leaves and roots) were also collected from the same Naganuma field in 2020 and 2021 (Table 1). The collected plant materials were stored at −80°C until their analysis.

Meta-RNA-seq analysis was basically conducted as has been described previously (Kondo et al., 2016; Lin et al., 2019). The total RNA from each wheat-leaf sample was extracted using NucleoSpin RNA Plant and Fungi Kit (Macherey and Nagel, Düren, Germany), following the manufacturer’s instructions. The obtained RNA fractions from the 2018 (three samples: KTH-18-6, −7, and −8) and 2019 (three samples: KTH-19-2, −4, and −6) samples were pooled into two groups based on the collection year – pool-18L (total 15.8 μg; RNA integrity number, RIN = 7.9) and pool-19L (total 7.8 μg; RIN = 7.7), respectively (Table 1). The two sample pools were subjected to ribosomal RNA (rRNA) depletion using a Ribo-Zero kit (Illumina, San Diego, CA, United States), and were subsequently used for the synthesis of a cDNA library using TruSeq Stranded Total RNA LT Sample Prep Kit (Plant; Illumina). The two cDNA libraries were then subjected to RNA-seq using the Illumina HiSeq 2000 platform (Illumina, 100-bp pair-end reads) performed by Macrogen Corp. Japan (Tokyo, Japan). After RNA-seq, the sequence reads were trimmed and de novo assembled using the CLC Genomics Workbench version 11 (CLC Bio-Qiagen, Aarhus, Denmark) using default parameters (minimum contig length = 200; mismatch cost = 2; length fraction = 0.5; similarity fraction = 0.8). Subsequently, assembled contigs were used as queries for basic local alignment search tool (BLAST) analyses (all contigs or contigs larger than 1.0 kb were subjected to BLAST-N or BLAST-X search, respectively) against the viral genome reference sequence (RefSeq) dataset of the National Center for Biotechnology Information (NCBI¹; E-value cut-off set >0.05). Quinvirus-like sequences from public datasets were also used as queries for BLAST searches against assembled contigs generated in this study to detect further virus-related sequences. The sequence reads (or assembled contigs) were mapped to the assembled virus or virus-like contigs using the Read Mapping algorithm using default or more stringent mapping parameters (match score = 1; mismatch cost = 2; length fraction = 0.5 or 0.9; similarity fraction = 0.8 or 0.95) in the CLC Genomics Workbench.

The total RNA was extracted from the plant materials using RNAiso Plus Reagent (TaKaRa Biotech. Co., Shiga, Japan) according to the manufacturer’s instructions. The cDNA was synthesized using Moloney-murine leukemia virus (MMLV) reverse transcriptase (Thermo Fisher Scientific, Waltham, MA, United States) with random primers [nonadeoxyribonucleotide mixture; pd. (N)9; TaKaRa Bio], and was then used as a template for PCR amplification with QuickTaq (Toyobo, Osaka, Japan). The conditions of PCR were as follows: 94°C for 2 min; then 30 or 35 cycles of 94°C for 10 s, 53°C or 57°C for 30 s, and 72°C for 1 or 2 min; and 72°C for 10 min. Alternatively, the RNA samples were subjected to one-step reverse transcriptase PCR (RT-PCR) analysis by using PrimeScript One-Step RT-PCR Kit Ver. 2 (TaKaRa), following the manufacturer’s instructions. The 5' and 3' rapid amplification of cDNA ends (5' and 3' RACE) analyses were conducted using the 5'-Full RACE Core Set and 3'-Full RACE Core Set (TaKaRa Bio), respectively, with virus-specific primers (Kondo et al., 2006; Supplementary Table S1). The wheat 18S rRNA gene was used as the reference target of an endogenous gene for the RT-PCR (Jarošová and Kundu, 2010). The primers used in the RT-PCR analyses are listed in Supplementary Table S1 and are available upon request. The selected PCR products were purified and subjected to Sanger sequencing to confirm their nucleotide sequences. For 5' and 3' RACE analyses, the obtained DNA products were ligated into the pGEM-T easy PCR cloning vector (Promega) and transformed into Escherichia coli strain DH5 alpha (TaKaRa Bio), and then plasmids were subjected to DNA sequencing.

Sequence data were analyzed using EnzymeX ver. 3.3.3² and 4peaks v1.8.³ BLAST or reverse BLAST searches were conducted using the GenBank database through the NCBI web site running on the non-redundant (nr) DNA and protein databases from the NCBI (nucleotide collection – nr/nt; transcriptome shotgun assembly – TSA; and expression sequence tag – EST; see Footnote 1). Sequence identities were also calculated using BLAST program available from NCBI (BLAST-N suite-2 sequence program; see Footnote 1). The conserved domains were searched using the NCBI conserved domain database.⁴ Pairwise sequence comparisons (PASCs) were made using the Sequence Demarcation Tool (SDT) ver. 1.2 (Muhire et al., 2014). A genome-based web tool for virus classification (pairwise sequence comparison – PASC) was used for the novel virus sequence analysis⁵ (Bao et al., 2014).

For the phylogenetic analyses, maximum-likelihood (ML) tree was constructed, as described in Kondo et al. (2019, 2020), with minor modifications. Multiple amino-acid sequence alignments were generated by using MAFFT online ver. 7⁶ (Katoh and Standley, 2013), with poorly aligned sites being removed using Gblocks online version 0.91b⁷ (Talavera and Castresana, 2007). Phylogenetic trees were then constructed using the PhyML 3.0 online program⁸ (Guindon et al., 2010; Lefort et al., 2017). Neighbor-joining (NJ) trees (Saitou and Nei, 1987) were constructed based on the MAFFT alignments masked with Gblocks. The reliability of the branches was obtained from 100 bootstrap replicates. The trees were refined using FigTree ver. 1.3.1 software.⁹

Small RNA deep sequencing was performed as described by Shahi et al. (2019). Two total RNA samples from 2018 (KTH-18-1 and −2) were mixed together (total 14.6 μg; RIN = 7.8) and subjected to small RNA seq (Table 1). After cDNA library preparation using a TruSeq Small RNA Library Prep Kit (Illumina), deep RNA-seq was conducted using an Illumina HiSeq 2500 (Illumina, 150-bp single-end reads) by Macrogen Corp., Japan. The raw sequence reads (total read number, 46,026,536; total read bases, 6,950-Mb) were trimmed of adapters and filtered for low-quality sequences and size range (15 to 32 nt, in length) using the CLC Genomic Workbench. The clean reads were subsequently mapped into the virus genomes. The virus-derived small RNA reads were used for further analysis using the program MISIS-2 (Seguin et al., 2016).

For soil transmission, the soils that were collected from the field with WYMV-infected wheat plants (KTH-20-1) were used as a virus inoculation source. Commercial culture soil (Nihon Hiryo Co., Ltd., Tokyo, Japan) was also used for non-infested healthy soils. Wheat seeds (cv. “Kitahonami”) were sown in plastic pots (7.5 cm in diameter) containing a mixture of the infested soil and the culture soil, and were grown in a growth cabinet (around 16°C, 12 h light/12 h dark) or in a greenhouse (non-controlled temperature conditions, around 15–25°C or occasionally slightly higher). At different periods after sowing, total RNA extracted from the roots of the plants was subjected to RT-PCR.

RNA sequencing analysis of the two pooled, rRNA-depleted RNA preparations of wheat-leaf samples (Table 1; Supplementary Figure S1) resulted in totals of 55,585,610 (5,614 Mb; pool-18L) and 64,002,754 (6,464 Mb; pool-19L) raw reads, respectively. De novo assembly using the CLC Genomics Workbench yielded 124,777 (pool-18L) and 131,004 (pool-19L) sequence contigs. These assembled contigs were then subjected to local BLAST-N (all contigs as queries) and BLAST-X (contigs larger than 1.0 kb as queries) analysis against the viral RefSeq collection. At least, 86 contigs (44 from pool-18L and 42 from pool-19L, ranging from 230 to 8,590 nt) were indicated as candidates for virus-derived sequences. The 41 contigs had a sequence similarity to WYMV (a bymovirus; 11 contigs), northern cereal mosaic virus (NCMV, a plant rhabdovirus; 2 contigs), and betaflexiviruses (at least 28 contigs; Supplementary Table S2). The remaining 43 contigs were related to fungal viruses of the families Narnaviridae, Botyourmiaviridae, and Mymonaviridae and some others (data not shown). Approximately, 1.1% of the reads (636,069 reads) in pool-18L and 0.9% of the reads (550,442 reads) in pool-19 were assigned to virus-related reads in each library. Almost all the virus-associated reads were derived from plant viruses – 43.5% (515,962 reads) for WYMV, 0.4% (4,216 reads) for NCMV, and 53.3% (632,472 reads; several unassembled reads related to a betaflexivirus strain might also be presented in the two pools, see below) for betaflexiviruses, with 2.9% (33,861 reads) representing others, likely derived from fungal-associated viruses (Supplementary Table S2 and data not shown).

We obtained at least 11 contigs of WYMV RNA segments from the pool-18L and pool-19L libraries. Among these, two contigs (Wh18L_c88 and Wh19L_c362) were nearly complete sequences of the WYMV RNA1 segment, while the RNA2 coding-complete sequence was generated using partially overlapping contigs in each dataset (Wh18L_c38/1970/1462 and Wh19L_c740/2865/2864, respectively; Supplementary Figure S2A and Supplementary Table S2). The WYMV sequences obtained from both libraries were almost identical, and the sequences were confirmed by mapping the reads (the representative contig sequences of RNA1 and RNA2 have been deposited in DDBJ under Accession Nos., LC63269 and LC63270, respectively; Table 2). The WYMV segments showed their highest sequence identities with a Japanese isolate Morioka (Accession No., AB627810) for RNA1 (99.6% nucleotide sequence identity) and a Chineses isolate (Shandong, KY354407) for RNA2 (99.0% nucleotide sequence identity).

To verify the presence of WYMV in the RNA samples used for the meta-RNA-seq analyses (Table 1), we performed RT-PCR, using the specific primer sets for RNA1 sequences (Supplementary Table S1). The typical yellow mosaic symptoms on the leaves of the cultivar “Kitahonami” samples are shown in Supplementary Figure S1. WYMV was detected in most of the 2018 “Kitahonami” samples, while one of two 2019 “Kitahonami” samples was positive (Figure 1). A phylogenetic analysis of the WYMV isolates, based on their RNA1 sequences, revealed that the obtained sequences belonged to genotype B, whose members are widely distributed in wheat-growing regions in northern Japan, including Hokkaido island (Ohki et al., 2014; Supplementary Figure S3A). In Japan, the WYMV isolates have been divided into three distinct pathotypes (I–III) based on their pathogenicity to wheat varieties, with pathotype II (genotype B) presenting in northern Japan, including Hokkaido island (Ohto et al., 2006). The WYMV isolate from the Naganuma field likely belongs to pathotype II, and symptoms could not be confirmed in the resistant wheat variety “Yumechikara,” in which a major quantitative trait loci, designated as Q. Ymym, against WYMV on chromosome 2D is known (Kobayashi et al., 2020; Figure 1; Supplementary Figure S3A and data not shown).

Two negative-sense RNA virus-like contigs (Wh18L_c253 and Wh19L_c20097) were identified from the pool-18L library (Table 2 and Supplementary Table S2). These two contig sequences are related to the NCMV (a planthopper-vectored rhabdovirus, genus Cytorhabdovirus, family Rhabdoviridae) and cover nearly the entire region of its genomic RNA (Supplementary Figure S2B). The concatenated contig sequence (Wh18L_c253/20097; deposited in DDBJ under Accession No., LC632071; Table 2) showed a high level of sequence identities with those of the Japanese (Accession No. AB030277) and Chinese (Hebei, GU985153) isolates of NCMV (98.5 and 93.8% nucleotide sequence identities, respectively). The distribution pattern of the mapped reads on the NCMV reference genome (4,186 mapped reads; Table 2) likely reflects the transcription gradient (a 3'–5' polar gradient of mRNA production), in which the transcripts of the N protein mRNA are the most abundant, with those of the downstream genes being at gradually lower levels (Dietzgen et al., 2017; Horie et al., 2021; Supplementary Figure S2B). The NJ tree, based on the L protein (RdRP), showed that the NCMV isolates are placed within a large sub-clade of cytorhabdoviruses (tentatively named subclade I), and formed a tight cluster together with other cereal rhabdoviruses, including barley yellow striate mosaic virus and maize yellow striate mosaic virus (Supplementary Figure S3B). Among the tested RNA samples in the pool-18L, NCMV was detected in one “Kitanokaori” sample (18_4-KTH; Figure 1). NCMV is known to be a viral agent of the stunted rosette symptom in wheat in East Asian countries (Tanno et al., 2000); however, the “Kitanokaori” plant (18_4-KTH) that was infected with the virus showed no stunting symptoms typical of the NCMV infection (data not shown). This may be due to the lower titer of NCMV in the host plants or other unknown reasons.

We identified 28 contigs in the two datasets (14 contigs each from the two sample pools, ranging from 203 to 3,828 nt) related to quinviruses (genera Foveavirus and Carlavirus, subfamily Quinvirinae) using the local BLAST analysis on the two datasets (Supplementary Table S2). Two concatenated sequences related to quinviruses were obtained by using the overlapped contigs, which share high nucleotide sequence identities at their overlapping regions with neighboring ones (Figure 2 and Supplementary Table S2). Subsequently, the entire genomic RNA region was identified by RT-PCR and Sanger sequencing using KTH-19-1 samples, following the 5' and 3' RACE analyses (Figure 2; Supplementary Table S1 and data not shown). The complete genomic sequences of the putative quinvirus were 8,412 and 8,411 nt in length, excluding the poly(A) tail; these were designated “strain a” and “strain b,” respectively (Table 2, Accession No., LC632066 and LC632067, respectively). In the 5' RACE clones, the GG (strain a, 9/16 clones) and G (strain b, 6/7 clones) nucleotides were mainly identified at the 5'-terminal ends of two genomic RNAs (data not shown). Both 5'-GG and 5'-G nucleotides at the 5'-terminal ends were most likely to be derived from 5'-end cap structure, as previously reported (Hirzmann et al., 1993; Kondo et al., 2014).

The genome organization of strain a was similar to that of foveaviruses and some others, with five ORFs, encoding replicase, TGB proteins, and CP, while strain b has an additional small ORF (ORF6) overlapping the CP gene (Figure 2). The presence of strain b with sequence variations within the 3'-terminal region were identified by the sequencing of the cloned 3' RACE products. The additional “A” residues presented at the potion 8,351 in some clones (UUAAAAA₈₃₅₁AnCCC…, additional A residues different from deposited one was underlined; of 12, 2, 1, and 1 3’-RACE clones had A₁, A₄, and A₅; data not shown). These sequence variants lacked the extra ORF (ORF6) and it seems to be a minor population in the tested wheat sample (KTH-19-1), as the majority of the RACE clones (8/12), which would retain ORF6, had no additional A. The replicase of WVQ does not contain an alkylation B (AlkB)-like domain, which is commonly present in betaflexiviruses that infect perennial plants, such as fruit crops (Kondo et al., 2013; Moore and Meng, 2019; Figure 2). They showed moderate levels of nucleotide (79.5% for the entire genome) and amino acid sequence identities (68.6∼91.7%) to each other (Table 3). Since different species of quinviruses should have less than about 72% nt identity or 80% aa identity for quinviruses (King et al., 2011), these two sequences likely represent the strains of a single quinvirus species. It should be noted that there were other overlapping fragments (contigs 133 and 431 in pool-18L, marked with a white circle in Figure 2) with high nucleotide sequence identities (>97.0%) to the strain a reference sequence (Supplementary Table S2 and data not shown), suggesting that this strain consists of multiple variants.

Along with these quinvirus sequences, some related contigs showed moderate levels of sequence similarity with the two strains (share 72.1–85.0% nucleotide sequence identity with the reference sequences), indicating the presence of a third strain (designated “strain c”) in the sample pools (Table 2; Supplementary Figure S2C; Supplementary Table S2 and data not shown). To confirm the presence of these three strains in the RNA samples, RT-PCR and subsequent direct Sanger sequencing were performed using strain-specific primer sets (Supplementary Table S1). A partial nucleotide sequence has been deposited as a representative of strain c sequence in DDBJ under Accession No., LC632068 (Table 2). All three strains were amplified in most of the tested samples in both pool-18L and pool-19L, although the amount of amplicons derived from strain c were relatively lower compared to those of the other two strains in some of the tested wheat plants (Figure 1). Unfortunately, the entire genomic sequence of strain c could not be determined via combinational RT-PCR and direct sequencing. Some overlapping contigs (contigs 936 and 1832 or 935 in pool-19L), likely derived from strain c, show high nucleotide sequence identities >95.9% to each other (Supplementary Table S2 and data not shown). Therefore, the mixed variant population of strain c likely existed in these pooled samples. This mixed population may have affected the direct sequencing analyses for the strain c. Further analyses are required to determine its entire genome.

BLAST-N analyses using the novel quinvirus sequences (strains a and b) as queries revealed that no significant genome-wide hits were detected (only for less than 12.0% for query coverage using other quinviruses; data not shown). PASC analysis of the genomic sequence of strain a showed 41.1% identity with that of a foveavirus – grapevine rupestris stem pitting-associated virus – which was the top hit sequence (Supplementary Figure S4A). Viruses among different genera in the Quinvirinae usually have less than about 45% nt identity in their genes (King et al., 2011), suggesting that WVQ represents a novel genus in the family. In the BLAST-P analyses, each predicted protein of the novel quinvirus (except for ORF6 in strain b) showed moderate levels of deduced amino-acid sequence identities (∼38.1–47.5%) with the corresponding proteins encoded by other quinviruses (mainly foveaviruses infecting fruit trees; Table 4). PASCs using STD also showed similar results (Figure 3 and Supplementary Figures S4B,C). Taken together, we propose that this newly discovered virus represents a virus species belonging to a novel genus in the family Quinvirinae.

To understand the phylogenetic relationships between WVQ and the known quinviruses, we constructed an ML tree using the amino-acid alignment of their entire replicase sequences. The constructed tree showed that the WVQ strains form a separate clade within the subfamily Quinvirinae. WVQ is distantly related to members of the three genera Foveavirus, Carlavirus, and Robigovirus, and other floating or unassigned quinviruses, such as SCSMaV (proposed genus “Sustrivirus”), banana mild mosaic virus (proposed genus “Banmivirus”), garlic yellow mosaic-associated virus, and yam virus Y (Gambley and Thomas, 2001; Thompson and Randles, 2001; Da Silva et al., 2019; Silva et al., 2019; Figure 4). Similar trends were also observed in the NJ trees based on CP and TGB, in which WVQ was placed separately from the three established genera (see also Ryu and Song, 2021; Supplementary Figures S5A,B).

Small RNA fractions of the wheat leaves (a pooled sample, KTH-18-1 and −2; Table 1) were deep sequenced to investigate the antiviral RNA silencing response. The small RNA reads (ranging 15–32 nt in length) mapped to the WYMV and WVQ were accounted for 0.5% (207,007 reads) and 0.2% (110,492 reads) of the total small RNAs, respectively. The small RNAs derived from both viruses were nearly equal positive and negative strands (45.9–59.6% for positive small RNA strands; Supplementary Figure S6A), predominantly 21 and 22 nt lengths for both strands (Figure 5A and Supplementary Figure S6B) and their 5'-terminal nucleotides were biased toward adenine (A) or uracil (U), which were particularly more prominent for 21 and 22 nt small RNAs in WVQ than those of WYMV (Figure 5B). Moreover, viral 21 and 22 nt small RNAs were distributed throughout the WYMV and WVQ genomes, with several distinct hotspots (Supplementary Figure S6C). Taken together, these analyses showed that WYMV- and WVQ-derived small RNAs possess the typical characteristics of viral small interfering RNAs (siRNAs; Li et al., 2017).

Our preliminary finding that WVQ was detected in WYMV-infected roots in the early spring (data not shown) when the temperature was relatively cool and thus the insects usually do not spread widely in the fields, suggests the possibility that WVQ may be transmitted to the roots through soil. Therefore, we aimed to examine the transmission of WVQ via soil under different conditions. Wheat seeds (cv. “Kitahonami”) were sown in pots containing soils from WYMV-infested fields or commercial soils as a control. One group of pots (three pots for each treatment) being kept under growth-cabinet control at approximately 16°C, which is generally allowed for WYMV infection, while the other group of pots was kept in the greenhouse (non-controlled temperature), where WYMV infection is usually more difficult (Zhang et al., 2021). Under the growth cabinet condition, both WVQ and WYMV were detected by RT-PCR in the roots of wheat after 3 months of growth, but not after 1 month (Figure 6A). After 3 months of growth in the greenhouse, WVQ, but not WYMV, was detected in the roots of the plants (Figure 6A). No systemic infection or virus symptoms were observed in the aerial parts of the plants in either trial (data not shown). Thus, it was revealed that WVQ is potentially soil-transmissible, and this virus may require lower temperature and/or other unknown conditions for systemic infection but not for root infection unlike WYMV.

WVQ strains were constantly detected by RT-PCR in the roots and leaves of the wheat plants (cv. “Kitahonami”) that were grown in WYMV-infested fields at Naganuma over 4 consecutive years (2018–2021; Table 1; Figures 1, 6B and data not shown). To investigate the occurrence of WVQ in other crops, the rye samples were also collected from the same field in 2020 and 2021 and were subjected to RT-PCR. The results showed that WVQ strains (except for the strain c) were also detected in both leaves and roots of rye plants along with co-infecting WYMV (Figure 6B), indicating that WVQ is able to infect rye plants.