With the development of urbanization and globalization, viruses previously confined to small, remote geographic areas is keeping spreading all over the world through global business^[1]. Classical viral diseases are invading novel populations. For example, the transmission of Middle East respiratory syndrome (MERS) coronavirus from camel to human in Arabic countries, and to remote areas, as far as Asia^[2-7], the transmission of yellow fever virus from Africa to Asia^[8-12], and the transmission of Zika virus from South America to the whole world^{[1, 13-17]}. Even in the area with inconvenient transportation, well-defined viruses can cause new headaches for public health, like the outbreak of Ebola in 2014-2016^[18-25].

The identification of epidemic pathogens is therefore challenging. Traditionally, serological tests and PCR-based methods are used as first-line identification methods in outbreak of epidemic diseases^[26]. However, these methods rely on pre-defined specific antibody, primer, probes, and do not work with new pathogens. The development of next generation sequencing (NGS), with unprecedented sequencing throughput, has promised to address these problems^[27-31]. With the use of random oligo nucleic acids instead of predefined primers, unbiased sequencing is able to identify any virus sequence, known or unknown, local or imported, or even the pathogens that have never been reported^[32].

However, unbiased sequencing or shotgun sequencing is susceptible to interference from environmental sequences. According to reported results as well as our own experience in processing traditional shotgun sequences, over 95% of success reads are commonly found to be host or environmental ribosomal RNAs^[33]. Consequently, background depletion techniques are particularly useful for samples where viral load is much lower than environment sequences. Some studies reported different methods for background depletion, but virome enrichment techniques are usually needed to access a broader set of viruses^[34].

In the present study, we aimed to develop a Viral Sequence Independent Targeted Amplification (VSITA) approach using a set of non-ribosomal and virus-enriched octamers (V8) and compare with traditionally used random hexamers (N6) in the reverse transcription step followed by library construction and NGS analysis, attempting to reduce the interference from unwanted RNA sequences and improve NGS-based virus identification in clinical samples.

First, we removed a total of 3, 985 hexamers that primed reverse-transcription of human or bacterial ribosomal RNA (Accession Number: NR_046235.3 of human 45S pre-ribosomal sequence, J01859.1 for Escherichia coli 16S ribosomal RNA, and NR_037007.2 for Staphylococcus aureus 16S ribosomal RNA) from a candidate dataset of 4, 096 (4⁶) random hexamers. Second, we constructed a single FASTA sequence file with all the virus-related sequences available from NCBI (8, 584 viral genome sequences) (release data 8 July 2014) and corresponding complementary sequences (17, 168 virus-related sequences in total). Third, we aligned the remaining 111 hexamers with both virus-related sequences, and scored each of them according to the numbers of perfect matching to either sequence. Forth, we then obtained the first 30 hexamers with the highest scores in matching to virus-related sequences (sense and antisense). Finally, we intended to add 2 random nucleotides to the 5' tail of the hexame to increase flexibility and diversity of targeted priming. Thus these 6+2N octamers (V8) (listed in Table 1) were used in reverse-transcription step for VSITA.

Forty-five archived clinical samples of different types (serum, feces, and throat swabs) were included in this study. These samples were previously tested by quantitative RT-PCR assay (qRT-PCR) (Supplementary Table 1 available in www.besjournal.com) and found to be positive for either of 14 different virus types/subtypes [Japanese encephalitis virus, Rotavirus, Norovirus, Human Adenovirus (hADV), Dengue virus, Influenzavirus A H1N1, Influenzavirus A H3N2, Influenzavirus A swirl H1N1 (swlH1N1), Enterovirus 71 (EV71), Coxsackievirus A16 (CA16), Cytomegalovirus (CMV), Human immunodeficiency virus (HIV), Hepatitis B virus (HBV), Hepatitis C virus (HCV)]. These 45 samples were used in parallel to compare the V8 and N6 enrichment performance of viral sequences and removal performance of ribosomal sequences in the step of reverse transcription followed by quantitative PCR (qPCR). Additional two separate groups, one group containing 10 serum samples from patients with fever of unknown origin (F1-F10) and one containing 10 feces samples from patients with diarrhea of unknown origin (D1-D10), were also included in this study. These 20 samples were used in comparison of V8 and N6 enrichment performance following NGS analysis.

For all the 65 samples, total RNA was extracted with RNease Mini Plus Kit (QIAGEN) following manufacturer's instruction. To evaluate the performance of VSITA, both 4, 096 random hexamers (N6) and VSITA primer set (V8) were used in parallel to synthesize first strand of cDNA with Superscript Ⅲ (Invitrogen) for each sample following manufacture's instruction. For second strand synthesis of each sample, Klenow fragment was added to first strand synthesis system for 2-hour incubation at 37 ℃.

For the 45 previously validated positive samples, both N6-primed and V8 primed cDNAs were evaluated by qPCR assays of ribosomal RNA and viral target sequence in parallel according to the protocols described in previous reports (Supplementary Table 1). Previously reported sequences of primers and probes were used in this study to ensure reliability of real-time PCR results (as listed in Supplementary Table 1). For comparison, relative quantity of target sequence using V8-primed cDNA vs. N6 counterpart was calculated as follows:

R: relative quantity of V8-primed target cDNA vs. N6-primed target cDNA; Ct_N6: Average Ct value of quantitative real-time PCR using N6-synthesized cDNA of four replicates; Ct_V8: Average Ct value of quantitative real-time PCR using V8-synthesized cDNA of four replicates.

For the 20 unidentified clinical samples, both V8 and N6 primed cDNAs were pre-enriched by multiple displacement amplification (MDA) using phi29 DNA polymerase (QIAGEN) according to manufacturer's instruction prior to library construction to guarantee enough amounts of templates for sequencing. Enriched cDNA products were evaluated with Agilent high sensitivity DNA assay on Agilent 2100 bioanalyzer. Fragmentation, end repair, adaptor ligation, size selection and purification were performed successively following instruction of Ion Torrent Hi-Q kit (Thermo Fisher Scientific). Adaptors with sequence barcodes were used to identify different libraries. After library evaluation with Agilent 2100 and real-time PCR, equal amount of each of the 40 libraries (20 samples for either V8 or N6 treatment) were grouped into four libraries and added to four 316v2 chips for sequencing. Sequencing was performed by 316 Chip Kit v2 on Ion Torrent PGM systems (Thermo Fisher Scientific) following instruction of Ion Torrent.

Raw data from sequencing was generated by software of Ion Torrent PGM system and downloaded in BAM format. The raw data was analyzed with our in-house software virus identification pipeline (VIP)^[35].

Relative quantity (the ratio) of V8-synthesized cDNA of target sequences vs. those N6-synthesized cDNA of the same sample was calculated to evaluate the ability of V8 primers in enriching target viral sequences and eliminating unwanted ribosomal sequences. Compared to N6, most V8-synthesized cDNA (39/45 for 18S, 31/45 for 28S) showed obvious decrease in the quantities of human 18S and 28S ribosomal RNA sequences (Figure 1A and 1B). The discrepancy occurred in partial samples infected with Rotavirus, hADV, swlH1N1, CMV, and HIV and all the samples infected with HBV and HCV. On the other hand, viral sequences were obviously enriched by V8 in most (38/45) samples (Figure 1C). The inconsistency occurred in partial samples infected with Norovirus, hADV, H3N2, coxsackievirus 16, HIV, HBV and HCV. In general, the relative quantities of 18S and 28S sequences were under 1 (Figure 1D), indicating the better performance of VSITA to inhibit the amplification of sequences from 18S and 28S. Meanwhile, the relative quantities of viral sequences were over 1 (Figure 1D), representing the better enrichment of viral genome sequences by VSITA.

Figure 1.  V8 versus N6 ratio of different RNA levels quantified by qPCR. (A) Human 18S ribosomal RNA level; (B) Human 28S ribosomal RNA level; (C) Viral RNA level; (D) Overview of relative quantity. The calculation of relative quantity of V8 versus N6 was showed in material and methods.

Fragment distribution of the enriched cDNA was evaluated by Agilent 2100 Bioanalyzer system. The obvious difference in the fragment distribution was observed between V8 and N6 priming. Taking the first-strand synthesis of 2 samples as examples for comparison (Figure 2A vs. 2B, 2C vs. 2D), the N6 priming showed higher diversity of the fragments, the fragments were more evenly distributed along the length axis (Figure 2A, 2C) while V8 priming tended to generate larger fragments more concentrated in some certain regions (Figure 2B, 2D). Results from other samples were showed in Supplementary Figure 1 (available in www.besjournal.com).

Figure 2.  Fragment length analysis with Agilent 2100 bioanalyzer with: (A) N6 library from sample F1; (B) V8 library from sample F1; (C) N6 library from sample D1; (D) V8 library from sample D1.

Figure Supplementary Figure 1.  Fragment length analysis with Agilent 2100 bioanalyzerwith samples: (A) N6 library from sample F1; (B) V8 library from sample F1; (C) N6 library from sample F2; (D) V8 library from sample F2; (E) N6 library from sample F3; (F) V8 library from sample F3; (G) N6 library from sample F7; (H) V8 library from sample F7; (I) N6 library from sample F10; (J) V8 library from sample F10; (K) N6 library from sample D1; (L) V8 library from sample D1; (M) N6 library from sample D5; (N) V8 library from sample D5; (O) N6 library from sample D6; (P) V8 library from sample D6; (Q) N6 library from sample D8; (R) V8 library from sample D8; (S) N6 library from sample D10; (T) V8 library from sample D10.

Of all the 10 samples from patients with fever of unknown origin, N6 primers approach succeeded to identify 185 reads of dengue virus (DENV Ⅲ) in sample F10 and covered 25.1% of dengue virus genome, while V8 approach identified 2, 257 reads and covered 94.88% of the same dengue virus genome (Table 2). The presence of dengue virus genome in sample F10 was later validated by Real-time PCR (Supplementary Table 1). Of 4 samples (F1, F2, F3, and F7, Table 2), a few reads were generated by V8 but not by N6 approach. However, very few reads from F1, F2, and F3 were found to match with Hantavirus and Phlebovirus and real-time PCR assay failed to amplify any sequence from these samples. For sample F7, V8 approach was only able to identify 46 reads (22.31% of genome) due to low quantity of dengue virus genome (real-time PCR Ct value of 35.14). Both methods failed in the detection of the remaining 5 samples.

Virus sequences were identified from 5 out of 10 diarrhea samples of unknown origin (D1, D5, D6, D8, and D10) by both V8 and N6 approaches. Both methods found co-infection of Echovirus 6 and Norovirus in sample D1 with high possibility (with genome coverage of 92.33% and 90.13% for Echovirus 6, 86.53% and 86.93% for Norovirus, respectively). Of the other four samples, infection of Coxsackievirus A6, Poliovirus, Coxsackievirus A4, Coxsackievirus B2 was found by both methods and validated by real-time PCR. For virus identification with NGS, V8 showed more viral reads hit and higher genome coverage in all the samples tested (Table 2).

Shotgun NGS sequencing with random primer- based library construction has become the preferred choice for investigation of unknown infections in many laboratories^[32]. However, clinical sample is always a combination of different genomes from the host, the pathogen, and environmental species. Under such circumstance, viral sequence is usually at disadvantage in random library construction, regardless of actual viral load^[36]. Ribosomal RNA, above all, is predominantly amplified in random sequencing, and therefore became major interference in shotgun NGS-based pathogen identification, especially in viral infection cases^[37]. Unlike prokaryotes and eukaryotes, there is little similarity between virus genomes. Accordingly, efforts have been made to enrich viral sequences in NGS investigation of complex clinical samples^{[34, 36]}. Endoh et al. has reported a set of 96 non-ribosomal hexanucleotides, which preferentially primes the virus genome RNA^[38], however, the viral enrichment abilities of both random hexamers and 96 non-ribosomal hexanucleotides primer sets were nearly comparative in model experiments as described in our recent report^[36]. We speculate that the reduced diversity in non-ribosomal hexamers (from 4⁶ = 4, 096 down to 96) might explain the limited coverage or sensitivity of viral pathogens found by 96 non-ribosomal hexamers.

To address this situation, we developed a novel VSITA approach by designing V8 primer set for virus-targeted selection in this study. After the removal of ribosomal RNA-targeted hexamers, we further narrowed virus-targeted hexamer sets down to an even smaller group containing only 30 virus-targeted hexamers. To expand the diversity of 30 primer set and improve detection sensitivity of viral pathogen, we intentionally added two random nucleotides (NN) to the 5' end of all the hexamers. As a result, the final V8 primer set was composed of 480 (30 × 4² = 480) non-ribosomal virus-targeted octamers. We infer that the virus-targeted specificity of V8 will not be obviously affected as the specificity of a primer is predominantly determined by its 3' end sequence.

Comparison of the virus-enrichment in cDNA synthesis between the V8 and random hexamers (N6) approaches was carried out by qPCR assays of 18S rRNA, 28S rRNA, and viral sequences. Our results showed V8 approach revealed improved virus enrichment and removal of human ribosomal RNA sequences in the majority of 45 archived clinical samples. Among a few samples (7/45), V8 approach showed inconsistent results (lower than N6) in virus (Norovirus, hADV, H3N2, coxsackievirus, HIV, HBV and HCV) enrichment. This phenomenon might be due to the difference in matrix constitution and targeted viral concentrations in extracted RNA of these clinical samples. Also, among some samples (6/45 for 18S, 14/45 for 28S), the performance of random hexamers was better in depletion of human ribosomal RNA than V8. Especially in the case of 28S sequences, N6 showed better performance than V8 in all of the 3 HBV-positive and 4 HCV-positive samples. The reason might be the nonspecific 28S amplification using cyber green-based qPCR analysis because subsequent melting curve analysis indicated they were nonspecific products. Though dsDNA viruses such as hADV and HBV were not primed in the step of reverse transcription by VSITA, they were involved in the virus-targeted enrichment in the PCR step due to the specificity of V8. More efforts will be needed to fine-tune the balance between virus-targeted sensitivity and specificity of V8 primer set in the future work. Nevertheless, given a large variety of virus-positive clinical samples of different types (serum, feces, and throat swabs) were tested, our data suggested VSITA is a worthwhile alternative to improve the enrichment of virus targets in the clinical samples.

Furthermore, 20 clinical samples with fever or diarrhea of unknown causes were used to evaluate the performance of V8 and N6 in parallel, the fragment length analysis of generated library was conducted with Agilent 2100 Bioanalyzer. The results showed that V8 produced longer fragments distribution and lower length diversity as compared to N6. This phenomenon might be attributed to the difference of V8 and N6 primer sets in the specificity and diversity. V8 has stronger specificity inclining to amplify limited viral genome regions, while N6 has higher diversity tending to generate fragments with a broad range. The sequencing results with ion torrent PGM system revealed an obviously superior performance using V8. By the use of V8 primers, one more case of dengue virus infection was found and validated using real-time PCR in serum sample tests. In the testing of fecal samples, both methods were able to report the causative viruses in several cases, including one case of multiple infections. However, V8 showed apparently higher genome coverage and more reads hit for all virus-positive samples in comparison with N6, except for one case with Norovirus infection, where N6 covered slightly larger region on the Norovirus genome (86.93% vs. 86.53%), although V8 achieved more reads hit. Overall, NGS data implied that V8 approach was superior to N6 in the identification of virus from clinical samples.

Other than virus-targeted amplification, efforts have also been made in pre-treatment of samples to improve virus identification with NGS. It is reported that ultracentrifugation of clinical samples improved virus discovery by removal of free RNA and DNA as well as bacteria^[39]. Pre-degradation of free nucleic acid with cocktail enzyme system has also been reported^[33]. The enrichment strategy in this study can be combined with those reported methods to provide an integrated and comprehensive enrichment solution.

In this study, only serum, swab, and feces samples were used to test the reliability of the VSITA approach followed by NGS. The performance of VSITA in other sample types, such as cerebrospinal fluid (CSF), urine, and other viral infection needs to be further validated in future study.

In conclusion, the VSITA approach designed in this study is demonstrated to possess higher sensitivity and broader genome coverage than traditionally used hexamers in NGS-based identification of viral pathogens directly from clinical samples.

All the authors approved the final manuscript and they have no conflict of interest to declare.

All aspects of the study were performed in accordance with national ethics regulations and approved by the Institutional Review Boards of National Institute for Viral Disease Control and Prevention, Center for Disease Control and Prevention of China.

We acknowledge the Pediatric Research Institute, Children's Hospital of Hebei Province, Shijiazhuang, Hebei, China for providing fecal clinical samples, the National Laboratory for Hemorrhagic Fever, National Institute for Viral Disease Control and Prevention (IVDC), Chinese Center for Disease Control and Prevention (CDC) for providing serum clinical samples.

1. Chuang TW, Lee KM, Lou YC, et al. A Point Mutation in the Exon Junction Complex Factor Y14 Disrupts Its Function in mRNA Cap Binding and Translation Enhancement. J Biol Chem, 2016; 291, 8565-74.

2. Kalia M, Khasa R, Sharma M, et al. Japanese encephalitis virus infects neuronal cells through a clathrin-independent endocytic mechanism. J Virol, 2013; 87, 148-62.

3. Verheyen J, Timmen-Wego M, Laudien R, et al. Detection of adenoviruses and rotaviruses in drinking water sources used in rural areas of Benin, West Africa. Appl Environ Microbiol, 2009; 75, 2798-801.

4. Kageyama T, Kojima S, Shinohara M, et al. Broadly reactive and highly sensitive assay for Norwalk-like viruses based on real-time quantitative reverse transcription-PCR. J Clin Microbiol, 2003; 41, 1548-57.

5. Krafft AE, Russell KL, Hawksworth AW, et al. Evaluation of PCR testing of ethanol-fixed nasal swab specimens as an augmented surveillance strategy for influenza virus and adenovirus identification. J Clin Microbiol, 2005; 43, 1768-75.

6. Pang Z, Li A, Li J, et al. Comprehensive multiplex one-step real-time TaqMan qRT-PCR assays for detection and quantification of hemorrhagic fever viruses. PLoS One, 2014; 9, e95635.

7. Zou S, Han J, Wen L, et al. Human influenza A virus (H5N1) detection by a novel multiplex PCR typing method. J Clin Microbiol, 2007; 45, 1889-92.

8. Poon LL, Chan KH, Smith GJ, et al. Molecular detection of a novel human influenza (H1N1) of pandemic potential by conventional and real-time quantitative RT-PCR assays. Clin Chem, 2009; 55, 1555-8.

9. Chen Y, Cui D, Zheng S, et al. Simultaneous detection of influenza A, influenza B, and respiratory syncytial viruses and subtyping of influenza A H3N2 virus and H1N1 (2009) virus by multiplex real-time PCR. J Clin Microbiol, 2011; 49, 1653-6.

10. Tan EL, Yong LL, Quak SH, et al. Rapid detection of enterovirus 71 by real-time TaqMan RT-PCR. J Clin Virol, 2008; 42, 203-6.

11. Cui A, Xu C, Tan X, et al. The development and application of the two real-time RT-PCR assays to detect the pathogen of HFMD. PLoS One, 2013; 8, e61451.

12. Sugita S, Shimizu N, Watanabe K, et al. Use of multiplex PCR and real-time PCR to detect human herpes virus genome in ocular fluids of patients with uveitis. Br J Ophthalmol, 2008; 92, 928-32.

13. Zhao JR, Bai YJ, Zhang QH, et al. Detection of hepatitis B virus DNA by real-time PCR using TaqMan-MGB probe technology. World J Gastroenterol, 2005; 11, 508-10.

14. Martell M, Gomez J, Esteban JI, et al. High-throughput real-time reverse transcription-PCR quantitation of hepatitis C virus RNA. J Clin Microbiol, 1999; 37, 327-32.