Metagenomic Analysis of Flaviviridae in Mosquito Viromes Isolated From Yunnan Province in China Reveals Genes From Dengue and Zika Viruses

More than 6,000 mosquitoes of six species from six sites were collected and tested for their virome using metagenomics sequencing and bioinformatic analysis. The identified viral sequences belonged to more than 50 viral families. The results were verified by PCR of selected viruses in all mosquitoes, followed by phylogenetic analysis. In the present study, we identified the partial dengue virus (DENV), Zika virus (ZIKV), and Japanese encephalitis virus (JEV) sequences in mosquitoes. Metagenomic analysis and the PCR amplification revealed three DENV sequences, one of which encodes a partial envelope protein. Two ZIKV sequences both encoding partial nonstructural protein 3 and one JEV sequence encoding the complete envelope protein were identified. There was variability in the viral titers of the newly isolated virus JEV-China/YN2016-1 of different passage viruses. The newly identified Zika virus gene from ZIKV-China/YN2016-1 was an Asian genotype and shared the highest nucleotide sequence identity (97.1%) with a ZIKV sequence from Thailand isolated in 2004. Phylogenetic analysis of ZIKV-China/YN2016-1 and ZIKV-China/YN2016-2 with known Flavivirus genes indicated that ZIKV has propagated in Yunnan province, China.


INTRODUCTION
Mosquitoes are important hosts and insect vectors for a range of infectious viruses including Japanese encephalitis virus (JEV), Zika virus (ZIKV), dengue virus (DENV), and Getah virus (GETV), which pose a significant threat to human health and represent an economic burden worldwide (Pham et al., 2005;Klungthong et al., 2007). Mosquitoes acquire viruses when consuming blood from hosts undergoing viremia. The viruses then replicate and propagate in the insect host before being introduced into further victims during biting and blood feeding (Ritchie et al., 2007;Motooka et al., 2015). Mosquitoes feed on a wide range of hosts including humans and other vertebrates, invertebrates, and plants (Shi et al., 2014). Yunnan province in China harbors a diverse range of mosquito-borne viruses (Feng et al., 2015); therefore, regional surveillance is imperative. The detection of viruses in mosquitoes is usually performed by using reverse transcription polymerase chain reaction (RT-PCR) and nested PCR approaches (Almeida et al., 2005;Houghton, 2009;Li et al., 2015;Sim et al., 2015). However, compared with Illumina sequencing (Alquezar-Planas et al., 2013;Cholleti et al., 2016;Ergünay et al., 2017), these traditional methods are time-consuming and labor-intensive to detect low-level viromes in mosquito vectors (He et al., 2013;Miesen et al., 2016). Therefore, metagenomic analysis of mosquitoes is likely to be of great value to avoid missing the detection of viruses with high infectivity and pathogenicity as well as the detection of previously unknown viruses.
The present study aimed to build a valid surveillance method to monitor the distribution of Flaviviridae from mosquitoes in Yunnan province, China and to provide useful insights into viral isolation, prevention, and control. Diverse and abundant viromes from mosquitoes isolated in Yunnan province were investigated using metagenomic sequencing and nested PCR. The presence of DENV, JEV, and ZIKV viruses was confirmed, and JEV was isolated. This preliminary exploration of the metagenomes of mosquito-borne viruses lays the foundation for further research on the territorial distribution, diversity, and surveillance of mosquito-borne viruses in China and other countries.

Mosquito Collection
More than 6,000 living or freshly dead female mosquitoes were collected in Yunnan province, China, during August and September 2016 (Figure 1). Sample I comprised Culex tritaeniorhynchus (C. tritaeniorhynchus). Sample II was a mixture of Armigeres obturbans (Ar. obturbans) and Aedes albopictus (Ae. albopictus). Sample III was a collection of Anopheles sinensis (An. sinensis), Uranotaenia hebes (U. hebes), and Armigeres durhami (Ar. durhami) ( Table 1). The mosquito samples were grouped based on species. C. tritaeniorhynchus was the predominant transmitting vector of JEV, and Ar. obturbans and Ae. albopictus were the predominant transmitting vectors of DENV and ZIKV (Byrd et al., 2013). The identified species of mosquitoes were collected separately and stored at −80 • C. This research was approved by the Ethics Committee and the Research Ethics Committee of Yanbian University Medical College. All experiments involving active viruses were performed in a biosafety level 3 laboratory.

Preparation of Mosquito Samples
Briefly, each mosquito sample was mixed and ground using a tissue grinder with SM buffer (50 mM Tris, 10 mM MgSO4, 0.1 M NaCl, pH 7.5) (MeilunBio, Dalian, China). To remove mosquito debris and other materials, the mixed samples were centrifuged at 13,000× g for 20 min, and the supernatants were used to extract the viral nucleic acid.

Extraction of Viral Nucleic Acid and Reverse Transcription
To remove the free nucleic acid and eliminate the contaminating host genomic DNA, 14 U Turbo DNase (Ambion, Austin, TX, USA), 25 U Benzonase Nuclease (Novagen, San Diego, CA, USA), 20 U RNase I (Fermentas, Ontario, Canada), and 10× DNase buffer (Ambion) were added to 127 µL of the supernatants to a final volume of 150 µL, followed by digestion at 37 • C for 1 h. The viral nucleic acid in the obtained products was extracted using a virus nucleic acid isolation kit (Bioer Technology, Hangzhou, China) according to the manufacturer's instructions. The total viral nucleic acids were reverse-transcribed using anchored random primers and Superscript III reverse transcriptase (Invitrogen, Carlsbad, CA, USA). The anchored random primers ( Table 2) were added separately to the viral nucleic acid and incubated at 75 • C for 5 min and then placed on ice for 5 min for denaturation. To obtain the reverse-transcribed product, 40 U of RNase OUT (Invitrogen, Carlsbad, USA), 200 U of SuperScript III reverse transcriptase (Invitrogen), 1 µL of 0.1 M dithiothreitol (DTT) (Invitrogen), 1 µL of 10 mM dNTPs (TaKaRa, Dalian China), 4 µL of 5× first-strand buffer (Invitrogen), and RNasefree H 2 O (TaKaRa) were added to a final volume of 20 µL and incubated at 25 • C for 10 min, followed by 50 • C for 60 min and then 75 • C for 10 min.

Synthesis of Double-Strand cDNA (dscDNA)
Before the synthesis of dscDNA, 1 µL of RNase H (TaKaRa) was added to the obtained products to degrade the free RNA. To synthesize dscDNA, anchored random primers were added and incubated 75 • C for 5 min and then on ice for 5 min for denaturation. Then, 1 µL of Klenow fragment (TaKaRa), 1 µL of 10 mM dNTPs (TaKaRa), 2 µL of 10× Klenow buffer (TaKaRa), and 6 µL of dd H 2 O (TaKaRa) were added and the samples were incubated at 37 • C for 60 min, followed by 75 • C for 10 min. Then, 0.5 µL of exonuclease I (TaKaRa), 1 µL of shrimp alkaline phosphatase (SAP, TaKaRa, Dalian, China), 5 µL of 10× phosphatase buffer (TaKaRa), and 24 µL of DEPC H 2 O (TaKaRa) were added and incubated at 37 • C for 60 min followed by 75 • C for 10 min to eliminate the phosphates and the free single-strand nucleic acid in the dscDNA reaction.

Sequence-Independent Single-Primer Amplification (SISPA) and Purification of PCR Products
To increase the quantity of the viral nucleic acids, SISPA was applied to amplify the dscDNA. The 50 µL reaction comprised 10 µL of dscDNA, 2 µL of barcode primer (

Metaviral Sequencing
The purified PCR products from the three samples were sent to the Wuhan Genome Institute (BGI, Shenzhen, China) for Illumina sequencing. Briefly, to obtain ∼180 bp DNA fragments, the purified PCR products were ultrasonicated, and dATP and Klenow fragment were added to produce 3 ′ -dA overhangs. To establish genomic DNA libraries, DNA fragments were bound to Illumina adaptors and amplified using PCR with adaptor primers. Amplicons were ligated to flow cells to which fluorescently labeled dNTPs were added. DNA sequences were identified by using the sequencing-by-synthesis method (SBS, Illumina). Base calling was monitored by the program GAPipeline (BGI), with default settings. No-calling reads and adaptor sequences were eliminated. The remaining sequences were assembled into contigs using SOAPdenovo software (BGI, Shenzhen, China). Contigs and sequences longer than 100 bp were defined as significant data for further in silico analysis.

Computational Analysis
The contigs and sequences were aligned using the blastx and blastn with the nonredundant and viral reference sequences in the GenBank database. (https://www.ncbi.nlm.nih.gov/ genbank/) BLAST hits with an E-value ≤ 10e −5 were considered significant. After eliminating the bacterial and eukaryotic sequences, we analyzed the virus-like sequences.

Identification of Detected Viruses
Based on the alignment outcomes of the viral contigs and the match position of the viral contigs with the corresponding viruses in GenBank, specific primers were designed and synthesized to identify the detected viruses. Viral nucleic acids were extracted by Bioer Technology (Hangzhou, China) and amplified using the designed primers and a PCR Master Mix (Tiangen, Beijing, China). To avoid false positives in Illumina sequencing, viral identification was performed three times in an independent manner.

Phylogenetic Analysis
The obtained PCR products were sequenced, and the sequences were aligned against sequences of the representative viruses using CLUSTAL W version 2.0. (Accession numbers are shown in the phylogenetic trees). To achieve better accuracy in phylogenetic analysis, neighbor-joining phylogenetic trees were produced using MEGA 7 with 1,000 bootstrap replicates Li, 2015).

Isolation of Viruses
Viral isolation was conducted with the positive ground mosquito supernatants. Briefly, the ground mosquito supernatants were diluted seven-fold with DMEM containing 2% FBS and cultured with BHK-21 cells for 5-7 days. Cultures were examined daily for evidence of a virally induced cytopathic effect (CPE). Cultures without a CPE were blind-passaged three times.
Identification of the Antiviral Effect Using an Indirect Immunofluorescence Assay (IFA) BHK-21 cells were seeded in 12-well plates, and after 24 h of culture, they were rinsed with phosphate-buffered saline (PBS) (Solarbio, Beijing, China) and fixed using methanol (Solarbio) for 10 min at room temperature. The cells were washed three times with PBS and then incubated in blocking buffer [3% bovine serum albumin (BSA) (BOSTER, Beijing, China) in Trisbuffered saline (Solarbio) and Tween20 (Solarbio)] for 1 h at room temperature. The cells were then incubated with an anti-E monoclonal antibody (Abcam, Cambridge, UK) (1:20) in PBS with 3% BSA overnight at 4 • C. After the cells were washed three times in PBS, they were incubated with a fluorescently labeled secondary antibody-fluorescein isothiocyanate-conjugated goat antimouse antibody (ZSGB-Bio, Beijing, China)-in the dark at 37 • C for 1 h. The nuclei were washed with PBS three times and stained with 4,6-diamidino-2-phenylindole (DAPI) (Solarbio). Slides were imaged under a fluorescence microscope (Olympus, Tokyo, Japan).

Observation by Negative-Stain Electron Microscopy
BHK-21 cells with suspected JEV-induced CPE were collected by repeated freezing and thawing three times, followed by centrifugation at 10,000× g for 15 min. The obtained virus suspension was centrifuged at 60,000× g for 4 h and the supernatant was removed gently and discarded. Subsequently, the viral precipitate was resuspended with an isometric mixture of 6.1% (v/v) pH 7.2 glutaraldehyde (HEDEBIO, Beijing, China) fixative and DMEM, 25 µL of which were added to the copper grid. After desiccation, one drop of 3% phosphotungstic acid (JINDU, Shanghai, China) was added for negative staining.
Before observation under an electron microscope (FEI, Hillsboro, USA), the grid was placed in an incubator (SANYO, Osaka, Japan) at 37 • C for desiccation.

Amplification and Phylogenetic Analysis of JEV Complete Genome
The JEV strain was isolated, and then the complete sequence was amplified using RT-PCR. The designed primers are displayed in Table 3. The complete JEV genome was assembled by PCR products using SeqMan version 7.1 (DNAStar, Madison, USA) with the default value. Subsequently, phylogenetic analysis was performed based on the complete sequence of JEV with other representative JEV strains reported in China and neighboring countries.

Detection of the Viral Titer of Different Passages
The viral titers of different passages of JEV (passage 5, passage 10, passage 20, and passage 30) were evaluated in BHK-21 cells in triplicate, with DMEM as a control. Cells were incubated with serial ten-fold dilutions of JEV or DMEM for 120 h at 37 • C. The viral titer TCID 50 (50% tissue culture infective dose) was calculated by using the Reed-Muench method (Krah, 1991), according to the cytopathic effect (CPE) in cells caused by viruses. It was evaluated based on the changes in cell culture morphology under a light microscope (OLYMPUS, Tokyo, Japan) that were observed at 100× magnification. The CPE of JEV or DMEM at each dilution was evaluated in eight replicates.

Variability of Envelope (E) Genes of Different Passage Viruses
Viral RNA was extracted from the supernatant using an RNA viral kit according to the manufacturer's protocol (Bioer Technology) after a CPE was induced in BHK-21 cells by different passage viruses (P5, P10, P20, and P30). The purified RNA was used as the template for cDNA synthesis using the SuperScript TM III first-strand synthesis system (Invitrogen) with the reverse primer JWR1 5 ′ -AGATCCTGTGTTCTTCCTCACCACCA-3 ′ , according to the manufacturer's instructions. The envelope (E) genes of different passage viruses were amplified using the paired primers shown in Table 4 and sequenced. The alignment of the four E genes was performed using MEGA 7.0.

Statistical Analysis
The data were analyzed using SPSS 17.0 software. Group comparisons for viral titers were carried out by using the t-test for independent means. Experiments were performed in triplicate and a minimum of three independent experiments were evaluated. The value of P < 0.05 was considered statistically significant.

GenBank Accession Numbers
The

Metaviral Sequencing and the Virome of Mosquitoes
To obtain the clean data of the virome of the mosquitoes, contaminating host sequences and barcode DNA were eliminated. A total of 22,964,305 reads were acquired by Illumina sequencing, with averaging read lengths of 186 nt ( Table 1). The amount of viral sequences was 17.00, 0.24, and 13.85% in sample I, II, and III respectively and the viruses could be clearly classified (Figure 2). Among them, 87.4% were  vertebrate viruses, 5.4% were plant viruses, and 1% were fungal viruses; insect viruses and phages accounted for 0.8 and 0.6%, respectively (Figure 2A). Most viral sequences in Sample I were from C. tritaeniorhynchus, the predominant JEV-transmitting vector in China, followed by A. obturbans. A. albopictus is the predominant ZIKV-transmitting vector. Viral sequences were classified at the family level ( Figure 2B), and Flaviviridae was the predominant viral family, together with certain other important viral families such as Circoviridae, Anelloviridae, Parvoviridae, and Peribunyaviridea. Additionally, numerous unclassified viral sequences, presumably belonging to novel viruses not previously identified in mosquitoes, were identified and are worthy of further study.

PCR Verification of the Metavirome Results
Viral sequences were assembled into contigs using SOAPdenovo. From the 13,645 assembled viral contigs, 32 DENV-like contigs with a read coverage of 184× (236-1408 nt), 13 ZIKVlike contigs with a read coverage of 34× (194-652 nt), and 26 JEV-like contigs with a read coverage of 52× (306-1,564 nt) were obtained. The DENV-like contigs shared 92-98% nt identity with known DENV sequences. The ZIKVlike contigs shared 95-99% nt identity with known ZIKV sequences. The JEV-like contigs shared 96-99% nt identity with known JEV sequences. To confirm the outcome of metavirome sequencing, specific primers were designed and synthesized to amplify the identified viruses ( Table 4). The following were the results of PCR identification in mosquito samples.

PCR Amplification of Dengue Virus
Dengue virus belongs to Flaviviridae, and 3,044 cases of dengue were reported in China in 2005-2012, among which 134 cases were reported in Yunnan province. Currently, there is no DENV-5 case reported in China. There have been DENV 1-4 outbreaks in China repeatedly and alternately. However, DENV-1 is predominant. (Xiong and Chen, 2014). Based on the alignment results of viral contigs, viral PCR amplicons from samples II shared a high identity with the dengue virus (DENV) genes, indicating that they belong to insect viral species. The results were verified by nested RT-PCR, which identified a 1,231 nt sequence (DENV-China/YN2016-1) encoding a partial nonstructural protein 4A, 4B, and 5 gene; a 491 nt sequence (DENV-China/YN2016-2) encoding a partial E protein, and a 939 nt sequence (DENV-China/YN2016-3) encoding a partial nonstructural protein 5 gene from DENV. Phylogenetic analysis showed that the newly identified DENV sequences were genotype IV. In addition, DENV-China/YN2016-2 was highly homologous to a DENV sequence from Thailand that was isolated in 1991, sharing 95.3% nt identity. This suggested that the viruses represented a DENV variant ( Figure 3A). Interestingly, the partial E gene in DENV-China/YN2016-2 may represent a new genotype IV gene in the DENV variant strain.

PCR Amplification of Zika Virus
The first imported case of ZIKV, pertaining to Flaviviridae, was confirmed in China in February 2016 Zhang et al., 2016). A total of 28 cases were reported in Zhejiang province, China in September 2016 (Deng et al., 2016;Li et al., 2016). In addition, 19 cases were reported in Guangdong province, China in 2016 (Sun et al., 2017). In line with the alignment results of viral contigs, viral PCR amplicons from sample II shared a high identity with ZIKV, also belonging to insect viral species. These results were also confirmed by means of nested RT-PCR, which identified 694 and 470 nt sequences, both encoding partial nonstructural protein 3 from ZIKV. These sequences were named ZIKV-China/YN2016-1 and ZIKV-China/YN2016-2, respectively. Phylogenetic analysis showed that the newly identified Zika virus sequence ZIKV-China/YN2016-1 was from an Asian genotype. The nucleotide identity analysis revealed that it shared the highest nucleotide sequence identity (97.1%) with a ZIKV sequence from Thailand that was isolated in 2004 ( Figure 3B).

PCR Amplification of JEV
Since 1951, 33,900 cases of JEV, which belongs to Flaviviridae, were reported in China. The most obvious feature for the territory distribution of JE in China in the 1970s or 2010 is that a small number of cases were reported in west and north China, whereas a high number of cases occurred in east and southwest China. Yunnan province is still a highly epidemic area for JEV. First reported in 1949, JEV has existed in China for more than 60 years. The prevalence of JEV is genotype I and III in China, and genotype I is predominant in Yunnan province (Zheng et al., 2012). According to the alignment results of viral contigs, viral PCR amplicons from sample I shared a high identity with JEV genes, indicating that they belong to insect viral species. The results were verified by means of nested RT-PCR, which identified a 1,500 nt sequence (JEV-China/YN2016-1) encoding a complete E protein. Phylogenetic analysis showed that the newly identified JEV sequence, JEV-China/YN2016-1, was genotype I (Figure 4). Nucleotide identity analysis revealed that it shared the highest nucleotide sequence identity (98.8%) with a JEV sequence from Laos that was isolated in 2009.

Viral Identification of JEV and Phylogenetic Analysis of its Complete Genome
After inoculation into BHK-21 cells, the obtained JEV strain induced an obvious CPE. Monolayer cells appeared rounded, with bulging, aggregation, shedding, and incomplete destruction of a single loose layer, with gaps between cells. The severity of the CPE was time-dependent. JEV was propagated in BHK-21 cells for 24 h. The levels of intracellular virus were measured by IFA using an E-specific monoclonal antibody. As expected, mockinfected cells showed no E expression, while JEV-infected cells were positive for E expression. Viral particles were observed using negative-stain electron microscopy. The particles were rounded with a diameter of ∼30-40 nm. Moreover, there were tiny protrusions on the surface, similar to JEV particles ( Figure 5). The phylogenetic analysis of the complete JEV genome showed that JEV-China/YN2016-1 still shared the highest nucleotide sequence identity (97.6%) with a JEV sequence from Laos that was isolated in 2009 (Figure 6).

DISCUSSION
Mosquitoes are the intermediate hosts of numerous viruses that infect humans, animals, insects, plants, and other species, and play a significant role in the prevalence of many infectious diseases Fisher et al., 2015;Shi et al., 2015). Before 2007, virus identification by traditional virological methods required several decades and involved laboratories worldwide (Krah, 1991). Compared with the traditional method, metagenomic sequencing by Illumina sequencing combined with a high throughput analysis is more efficient, resulting in the identification of many viruses from different species in recent years by several laboratories (Wilson et al., 2015). Moreover, metagenomic sequencing has revealed the high abundance of viruses in mosquitoes (Bolling et al., 2015;Fauver et al., 2016;Shi et al., 2017). Metagenomic sequencing has also contributed to the discovery of novel viruses and their identification and characterization (Carissimo et al., 2016;Frey et al., 2016). In the present study, more than 50 viral families were detected by Illumina sequencing in mosquito samples. However, there was a considerable difference in the percentage of viral sequences obtained from the three mosquito groups. It could be related to the biology of the vectors (Atoni et al., 2018;Xia et al., 2018;FIGURE 6 | Phylogenetic trees of the complete JEV genome. Phylogenetic trees based on the complete genome of JEV. The trees were constructed using the p-distance-based neighbor-joining method in MEGA 7.0 software. Bootstrap values were calculated with 1,000 replicates. Black solid circles indicate the complete genome of JEV isolated in this study. Zakrzewski et al., 2018). Still, new findings of viral distribution and evolution in Yunnan province were obtained. RT-PCR and nested-PCR were used to confirm the results of Illumina sequencing. The detection of DENV in the mosquito samples showed that the virus still exists in Yunnan province, and many strains of DENV have been uncovered by previous studies Shi et al., 2014;Zhang et al., 2014;Wang et al., 2016). However, the detection of new DENV strains is important, as it suggests a new prevalence of DENV in this area. Moreover, the identification of the Zika virus genes is significant, suggesting that increased environmental surveillance is required because of its strong infectivity to humans, especially pregnant women and infants. Meanwhile, the presence of unidentified viruses in the Illumina sequencing results may ascribed to inadequate mosquito sampling and limited collection sites.
ZIKV and DENV were detected simultaneously, indicating that these viruses co-circulate in Yunnan province and suggesting that co-infection by the viruses is possible. Zika virus pathogenesis can be enhanced by preexisting dengue virus (Bardina et al., 2017). The simultaneous detection of the genes of two viruses may represent a new challenge to the control of ZIKV. We performed the isolation of DENV and ZIKV viruses in Vero and C6/36 cells besides BHK-21 cells as well as the suckling-mouse brain was also used in viral isolation. It was unfortunate that we could not isolate DENV and ZIKV.
The detection and isolation of JEV indicated that the virus still exists in Yunnan province. Analyzing the variant sites in the E gene and the consequences for the sequence of the encoded protein, we discovered that amino acid 138 in the E proteins of viruses at P5 and P30 is K (lysine). However, amino acid 138 in the E proteins of viruses at P10 and P20 is E (glutamic). Lysine at position 138 in the E protein is an attenuated viral site (Ni et al., 1994). The results showed that JEV-China/YN2016-1 has undergone mutations when passaged in BHK-21 cells, suggesting the possibility of increased threat to human health induced by JEV variation. The phylogenetic analysis based on the complete genome and the E gene of JEV-China/YN2016-1 showed that it shared the highest nucleotide sequence identity with a JEV sequence from Laos that was isolated in 2009, indicating the possible threat of JEV transboundary transmission.
This is the first report of the detection of Zika virus genes in Yunnan province, China. The newly identified DENV-China/YN2016-2 is a novel genotype IV gene in the DENV mutant strain. Metagenomic analysis of the viromes of mosquitoes in Yunnan province, combined with nested PCR, revealed that the distribution of viruses is dependent on both the mosquito species and the geographical location. The present work explored only a small portion of the virome of mosquitoes prevailing in this region of China. A much larger study on mosquitoes in other provinces in China and from other countries is needed to assess the diversity of mosquito-borne viruses. In conclusion, our findings provide a useful insight into the viral isolation and the characterization of Zika and dengue viruses and could hence be applied to larger studies in the future.

AUTHOR CONTRIBUTIONS
PX, HJL, and NJ conceived and designed the experiments. PX, JH, CL, SW, YL, MW, and NJ performed the experiments. PX, HL, JR, and NJ analyzed the data. YZ, XG, HZ, HJL, and MT contributed reagents, materials, and analysis tools. PX and NJ wrote the paper. All authors read and approved the final version of the manuscript.