In July 2015, leaf samples with mosaic and curling symptoms were collected from watermelon plants of cv. Hēixiùlóngjǔanfēng grown in the field in Hénán Province, China. Subsequently, in June 2016, watermelon samples with virus-like symptoms were randomly collected from Kāifēng County, Hénán Province. Samples were stored at −70°C.

Symptomatic leaves from cv. Hēixiùlóngjǔanfēng (sample KF-1) and healthy leaves (sample KF-0) were selected for small RNA-sequencing (sRNA-seq) in 2015; leaves coded KF-15 were used for RNA-sequencing (RNA-seq) in 2016. Total RNAs were extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol but modified to repeat the chloroform extraction and subsequent centrifugation step (Liu et al., 2015). A BioAnalyzer 2100 (Agilent Technology, Santa Clara, CA, USA) was used to assess the quality of the purified RNA, then the purified RNA was quantified using a Nanodrop ND-2000 spectrophotometer (Infinigen Biotechnology, City of Industry, CA, USA). The Ribo-Zero Magnetic Kit (Plant Leaf) (Epicentre, Madison, WI, USA) was used to deplete the ribosomal RNA from transcriptomes for RNA-seq, then the library was constructed with a TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA), and a TruSeq Small RNA Sample Prep Kit (Illumina) was used to construct the sRNA libraries. Subsequently, the Illumina Hiseq2500 platform (1 × 50 bp read lengths) was used for sRNA-seq (SinoGenoMax, Beijing, China), and the Hiseq 4000 platform was used for RNA-seq with PE150 bp (Biomarker Biology Technology, Beijing, China). Raw reads from the Illumina platform were trimmed of adaptor sequences and low-quality reads by the CLC Genomics Workbench 9.5 (Qiagen, Valencia, CA, USA). The Velvet program was used for de novo assembly with a k-mer of 17 for small RNA data (Zerbino and Birney, 2008; Wu et al., 2010). The CLC Genomics Workbench 9.5 was used for de novo assembly for RNA-seq data. The assembled contigs were subsequently screened against the NCBI databases using a BLASTn and BLASTx search with standard parameters.

To verify the presence of viruses in field samples, we used RT-PCR and specific primers designed on the sequences of assembled contigs (Table S4). The sequences of the extreme ends of the genomic RNAs were determined employing the 5′- and 3′-RACE System for Rapid Amplification of cDNA Ends kits (Thermo Fisher Scientific, Waltham, MA, USA) (Than et al., 2016). The primers and other oligonucleotides used for RT-PCR and RACE analyses are listed in Table S3. The PCR products were separated by electrophoresis in 1.0% agarose gels and visualized by ethidium bromide staining, then purified with the Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI, USA). The final products were cloned into the PEASY-T5 cloning vector (TransGen Biotech, Beijing, China), and then inserted into Trans-T1 Chemically Competent Cells (TransGen). At least eight clones from each subclone were sequenced (Beijing Genomics Institute, Beijing, China).

Sequences were analyzed and assembled using Vector NTI 11.5, then submitted to the GenBank database in NCBI using the software Sequin. The reconstructed genome was subjected to standard sequence analyses: (I) prediction of ORFs using the ORF Finder program in NCBI (http://www.ncbi.nlm.nih.gov/projects/gorf/); (II) identification of conserved and functional domains of the predicted proteins in WCLaV-1 and WCLaV-2 using the Conserved Domain Database (CDD) in NCBI (https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml) (Marchler-Bauer et al., 2015) and the SMART tool (http://smart.embl-heidelberg.de/) (Letunic et al., 2015); analyses of the core RdRp and NP domains in other (-ss) RNA viruses using the CDD; (III) multiple sequence alignment of the core motifs in RdRp and the conserved nucleotide sequence at the 5′- and 3′-end by CLC Genomics Workbench 9.5; (IV) identity analyses using the needle program in WebLab (http://weblab.cbi.pku.edu.cn/); (V) multiple sequence alignments first using the Clustal W method and then MEGA 6.0 (Tamura et al., 2013) for phylogenetic tree construction using the neighbor-joining method with 1,000 bootstrap replicates, the Poisson model and pairwise gap deletion options. Abbreviations and accession numbers of viruses used in this study are listed in Table S5.

The crude sap, extracted from leaves (about 1 g) infected by WCLaV-1 and homogenized in 10 mL 0.01 M phosphate buffer (PB) (pH = 7.2), was used for mechanical transmission tests on four-leaf-stage seedlings of watermelon cv. Hēixiùlóngjǔanfēng. Plants non-inoculated and inoculated only by PB were used as negative controls. After inoculation, the seedlings were maintained in a clean growth chamber at 25°C and 80% relative humidity with 16-h light/8-h dark. The mechanically inoculated plants were carefully observed for symptoms and subsequently tested for virus infection by RT-PCR assays with the virus-specific primer pairs listed in Table S4 and sequenced at 17, 42, and 67 days post inoculation (dpi). The mechanical transmission experiments were repeated three times with 20, 23, and 46 healthy watermelon plants inoculated, respectively. The cDNA was synthesized using Moloney murine leukemia virus (MoMLV) reverse transcriptase (Promega), and the PCR protocol employed 32 cycles after an initial denaturation at 95°C for 5 min. Each cycle consisted of denaturation at 95°C for 45 s, primer annealing for 45 s at 55°C and an extension for 1 min at 72°C, followed by an elongation step of 10 min at 72°C. The products were separated using 1.0% agarose gel electrophoresis.

Our field investigation revealed diseased plants of watermelon with various virus-like symptoms, including leaf crinkling, mosaic, and bunchy top (Figures 1A,B); no symptoms were observed on the seemingly healthy plants (Figure 1C). To identify the viruses present in symptomatic samples, we subjected samples KF-1 and KF-15 to sRNA-seq and RNA-seq, respectively. Healthy plant KF-0 was also small RNA-sequenced. For sRNA-seq, a total of 22,235,704, and 20,667,689 reads with length between 18 and 27 nucleotides (nt) were obtained in sample KF-0 and KF-1, respectively. After de novo assembly and BlastX and BlastN searches against nr and nt databases, ZYMV (genus Potyvirus), WMV (genus Potyvirus), melon aphid-borne yellows virus (MABYV, genus Polerovirus), CiLCV (genus Deltapartitivirus), and four contigs homologous to viruses in the genera Phlebovirus and Ophiovirus were identified in sample KF-1 (Table S1). In contrast, contigs from KF-0 did not show any significant similarity with previously reported viruses. For RNA-seq, a total of 16,547,389 reads with a length of 150 nt were obtained after filtering the reads mapping to the watermelon genome (Guo et al., 2013). Interestingly, we identified CiLCV, CGMMV (genus Tobamovirus), watermelon virus A (WVA, family Betaflexiviridae), and three contigs homologous (high similarity with that identified in KF-1) to phleboviruses and ophioviruses in the sample KF-15 collected in 2016 (Table S1). Reverse transcription-polymerase chain reaction (RT-PCR) using specific primers (Table S4) and sequencing of the amplified products confirmed the presence of these contigs.

The genomes were obtained by RT-PCR and RACE-PCR (Figure S1) using virus-specific primers designed on the sequences of these contigs (Tables S3, S4). The assembled genomes are composed of three RNA segments, with complementary-sense RNAs (vcRNAs) of these two viruses predicted to contain a single ORF (Figures 2A,B). The nucleotide sequences reported here have been deposited in GenBank (respectively, accession KY781184, KY781185, KY781186 for RNA1, RNA2, and RNA3 of WCLaV-1; KY781187, KY781188, KY781189 for RNA1, RNA2, and RNA3 of WCLaV-2).

RNA1 for both viruses is the largest segment and has only one ORF (ORF1): 6,636 nt and 6,680 nt for WCLaV-1 and -2, which start, respectively, at position 6,584–6,586 and 6,634–6,636 and end at position 38–40 and 40–42. ORF1 encodes a putative RNA-dependent RNA polymerase (RdRp, p1) of 2,182 amino acids (aa) and 2,198 aa with a predicted molecular weight of 252.76 and 253.63 kDa, respectively, for WCLaV-1 and -2 (Figures 2A,B). A Blastp search with the algorithm PSI-BLAST using the deduced aa sequences of p1 show that they have about 25% identity to the L proteins of viruses in the genus Phlebovirus (Table S2). Additionally, SMART analyses recognized the conserved domain Bunya_RdRp (Pfame-value: 5.7e-51, 1.1e-47) in the p1 protein of these two viruses (Table S2). Interestingly, the conserved domain RRM_SF super family (RNA recognition motif superfamily, also known as RNA binding domain) was predicted by the Conserved Domain Database for p1 in WCLaV-2 (interval: 1,386–1,443; E-value: 4.22e-04), but was not found in WCLaV-1 (Marchler-Bauer et al., 2015). The alignment of RdRp of WCLaV-1, WCLaV-2, and other (-ss) RNA viruses indicates that the amino acid sequences of their RdRps contain the six conserved motifs (premotif A and motifs A–E), which represent highly conserved regions of the RdRp of the members in the order Bunyavirales (Figure 2C). Both WCLaV-1 and WCLaV-2 have the motif A (DATKWC), motif B (QGILHYTSS), and motif D (KS). Motif C includes the SDD sequence. The motif E tetrapeptide E(F/Y)xS, which is specifically conserved in the polymerases of segmented negative-sense RNA viruses (van Poelwijk et al., 1997; Kormelink et al., 2011), was also found in the proteins encoded by RNA1 of these two viruses (Figure 2C). In addition, we identified three basic residues (K, R, and R/K) in premotif A and a glutamic acid (E) downstream of premotif A (Figure 2C), which are also conserved in bunyavirus RdRps (Bruenn, 2003; Elbeaino et al., 2009). Moreover, the endonuclease domains H₆₆-D₇₈-PD₉₅₋₉₆-E₁₀₇-K₁₂₆, similar to that in phlebo- and tenuiviruses and differing from the H-D-PD-DxK-T in orthotospoviruses, were also observed in the N-terminus of WCLaV-1 and WCLaV-2. These domains are probably involved in cap-snatching, a genome expression strategy to cap viral mRNAs proposed for members of the previous family Bunyaviridae and genus Tenuivirus (Biswas and Nayak, 1994; Reguera et al., 2010). In addition, the core RdRp domain of WCLaV-1 and WCLaV-2 is closely related to that of the tenuiviruses and phleboviruses, respectively; WCLaV-1 shares 27.6% aa sequence identity with Salehabad virus (SALV, genus Phlebovirus) and WCLaV-2 shares 27.2% aa sequence identity with rice stripe virus (RSV, genus Tenuivirus), both remarkably higher than with viruses in other genera, but far lower than the 75% identity between these two viruses.

RNA2 of WCLaV-1 (1,409 nt) and WCLaV-2 (1,423 nt) may, respectively, encode a protein (p2) of 44.29 and 46.11 kDa, spanning position 1,356–1,358 and 1,366–1,368 (AUG start codon) to position 196–198 and 155–157 (UGA stop codon) of vcRNA (Figures 2A,B). Interestingly, stretches highly rich in U were found in the 5′-terminus of RNA2 of both viruses. Blastp searches showed that the deduced aa sequences of RNA2 have 24% identity (28/116; E = 1.4) with the 54 kDa protein of citrus psorosis virus (CPsV, genus Ophiovirus), demonstrating a putative movement protein (MP) function, for WCLaV-1 and 19% identity (25/133; E = 1.4) to the 58 kDa MP of blueberry mosaic associated virus (BIMaV) in the genus Ophiovirus (Robles Luna et al., 2013; Borniego et al., 2016) for WCLaV-2 (Table S1). Therefore, the deduced protein p2 of WCLaV-1 and WCLaV-2 may be the putative MP. Nevertheless, no significant conserved domains in p2 of these two viruses have been identified using CDD and SMART analyses (Table S2). Pairwise comparisons showed that the aa sequence identities between the p2 of WCLaV-1, WCLaV-2 and other ssRNA viruses ranged, respectively, from 2.2 to 18.3% and 9.7 to 18.5%; the highest aa identities, 18.3 and 18.5%, respectively, were found with proteins encoded by RNA2 of groundnut bud necrosis virus (GBNV, the genus Orthotospovirus) (Table 1). Compared with the aa identity of 46.1% between these two viruses, the identities with known viruses are very low. Altogether, these data revealed that p2 encoded by RNA2 of these two viruses has unique features worth further study.

RNA3 of WCLaV-1 is 1,310 nt long and 1,260 nt for WCLaV-2. The ORFs start at position 1,251 and 1,198 nt, respectively, with an AUG start codon and terminate with an UAA stop codon at position 208 and 158 nt, respectively, putatively encoding a nucleoprotein (NP, p3) of 347 a and 346 aa with a molecular mass of 38.75 and 38.68 kDa, respectively (Figures 2A,B). It is worth noting that U-rich regions are also observed at the 5′-terminus in RNA3 of these two viruses, similar to the case in RNA2. According to the CDD and SMART analyses, RNA3 encodes proteins that contain the conserved domains Tenui_N super family (accession cl05345), the typical domains in nucleocapsid proteins of viruses in the genera Tenuivirus and Phlebovirus, with a Pfame-value of 5.4D-16, 3.2e-17, respectively, distinct from those in other genera. Moreover, the highest amino acid identities for the nucleocapsid protein, 17.2 and 17.4%, respectively, were found with SALV and sandfly fever Naples virus (SFNV, genus Phlebovirus) for WCLaV-1 and WCLaV-2, far lower than the 47.9% aa identity between these two viruses for nucleocapsid protein (Table 1). They share the highest aa identities, 21.2 and 23.2%, respectively, with the core NP domains of rice grassy stunt virus (RGSV, genus Tenuivirus) for WCLaV-1 and Punta Toro virus (PTV, genus Phlebovirus) for WCLaV-2 (Table 1).

In addition, the highly conserved nucleotides 5′-ACACAAAG-3′ at the 5′-terminus were found in each RNA1 segment, and the conserved nucleotides 3′-UGUGUUUC-5′ at the 3′-terminus were detected in each RNA2 and RNA3 segment (Figures 2D,E). The 5′- and 3′-terminal regions exhibited obvious complementarity, which is a common characteristic of genomic segments in tenuiviruses and phleboviruses (Elliott and Blakqori, 2011).

To reveal the relationships of WCLaV-1 and WCLaV-2 with known (-ss) RNA viruses, we constructed phylogenetic trees based on the amino acid sequences of the putative RdRp (Figure 3), MP (Figure S2A) and NP (Figure S2B). According to the generated trees based on amino acid sequences of RdRp, WCLaV-1, and WCLaV-2 are most closely related to viruses in the genera Tenuivirus, Phlebovirus, Goukovirus, Phasivirus, the family Phenuiviridae (Figure 3). However, phylogenetic analyses based on the amino acid sequences of RNA2-encoded proteins indicate these two viruses consistently cluster with emaraviruses and ophioviruses (Figure S2A). WCLaV-1 and WCLaV-2 form a distant and separate clade with the selected (-ss) viruses in all the phylograms, revealing that they may belong to a novel genus.

The mechanical transmission experiment with WCLaV-1 was done three times with 20, 23, and 46 healthy watermelon plants inoculated, respectively. Each time, only one plant, respectively, of the inoculated samples was positive for WCLaV-1 infection by RT-PCR, thus conclusively showing infection rates of 5.0, 4.3, and 2.2%, respectively. For example, at 17 dpi, WCLaV-1 was detected only in watermelon plant No. 4 using RT-PCR with primer pair NP-1F/R, and the positive band was further identified by sequencing (Figure 4A, Table S4). Moreover, four uninoculated leaves of watermelon plant No. 4 were confirmed as infected by WCLaV-1 by RT-PCR detection with primer pairs NP-1F/R and MP-1F/R (Figure 4B), thus demonstrating that the watermelon plant was systemically infected by WCLaV-1. WCLaV-1 was not detected in the control plants non-inoculated and inoculated with phosphate buffer (PB). By 17 dpi, new leaves of watermelon seedlings had developed obvious mosaic and crinkling similar to the symptoms observed on plants in the field (Figure 4C). Subsequently, at 42 dpi, watermelon plant No. 4 displayed obvious crinkling and dwarf symptoms, which were lacking in the non-inoculated plant (Figure 4D). At 67 dpi, leaf crinkling was also observed (Figure 4E).

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.