Large-scale single-virus genomics uncovers hidden diversity of river water viruses and diversified gene profiles

Abstract

Environmental viruses (primarily bacteriophages) are widely recognized as playing an important role in ecosystem homeostasis through the infection of host cells. However, the majority of environmental viruses are still unknown as their mosaic structure and frequent mutations in their sequences hinder genome construction in current metagenomics. To enable the large-scale acquisition of environmental viral genomes, we developed a new single-viral genome sequencing platform with microfluidic-generated gel beads. Amplification of individual DNA viral genomes in mass-produced gel beads allows high-throughput genome sequencing compared to conventional single-virus genomics. The sequencing analysis of river water samples yielded 1431 diverse viral single-amplified genomes, whereas viral metagenomics recovered 100 viral metagenome-assembled genomes at the comparable sequence depth. The 99.5% of viral single-amplified genomes were determined novel at the species level, most of which could not be recovered by a metagenomic assembly. The large-scale acquisition of diverse viral genomes identified protein clusters commonly detected in different viral strains, allowing the gene transfer to be tracked. Moreover, comparative genomics within the same viral species revealed that the profiles of various methyltransferase subtypes were diverse, suggesting an enhanced escape from host bacterial internal defense mechanisms. Our use of gel bead-based single-virus genomics will contribute to exploring the nature of viruses by accelerating the accumulation of draft genomes of environmental DNA viruses.

Keywords: DNA viruses; droplet microfluidics; environmental viruses; single-virus genomics; whole genome amplification.

Figure 1

Gel beads enable high-throughput WGA of viruses at the single-particle level. (A) Workflow of single-virus genomics. Viruses are encapsulated into 30-μm of microfluidic droplets with ultra-low melting temperature agarose. After the solidification of agarose, each virus captured in a gel bead proceeded to undergo capsid lysis and WGA. Gel beads with amplified DNA are isolated with FACS and proceeded to the library preparation and high-throughput sequencing. (B) Microscopic image of gel beads after WGA. The amplified viral DNA is stained with SYBR Green I, with a scale bar of 100 μm. (C) Comparison of the number of sequence reads mapped to the Lambda and Charomid reference sequences. Samples were labeled with the reference with the most sequence reads to be mapped. Only 4 of the 96 gel beads contained >1% of the other viral sequence reads, suggesting a low risk of DNA contamination. (D) Genome coverage by the sequence reads derived from each gel bead.

Figure 2

Sequence quality evaluation of vSAGs and vMAGs collected from river water. (A) Negative staining TEM images of viral suspensions were taken at a magnification of 100 000×. Various shapes of virus-like particles were recovered from river water. Most particles are dispersed, whereas some seem to be physically attached. The scale bar is 500 nm. Viruses were collected using two different methods ((I) suction filtration and (II) flocculation with FeCl₃), and single-virus genomics was performed with each viral suspension. Metagenomic DNA was extracted from the viral suspensions collected by method (I). After high-throughput sequencing and the construction of vSAGs and vMAGs, (B) the sequence length, (C) the collected number and the estimated quality by CheckV, and (D) the sequence similarity against reference sequences in the IMG/VR v3 database were evaluated. If the average nucleotide identity (ANI) was >95% and the alignment fraction (AF) was >85%, the viral sequences were judged to be known sequences (dashed square).

Figure 3

Protein-sharing network of 1431 vSAGs and 100 vMAGs. (A) Protein-sharing network analysis of 1431 vSAGs and 100 vMAGs with RefSeq reference database using vConTACT2. Only the sequences sharing at least one PC with other sequences are represented. (B) Relative abundance of vSAGs and vMAGs in the metagenomic raw reads. The vertical axis is the number of mapped metagenomic sequence reads to the representative sequences, expressed as the logarithm of reads per base. (C) Metagenomic sequence reads recruitment patterns (also referred to as diversity curves) for the representative sequences of each VC. Curves represent the percentage of recruited sequence reads at each nucleotide identity value [13].

Figure 4

Identification of AMGs commonly detected from multiple VCs. (A) The distribution profiles of AMGs among vSAGs. Each row represents a single vSAG, and each column represents a single AMG. Straight lines separate each VC, and subclusters within a VC are separated by dashed lines. The AMG categories were determined by DRAM-v. The left side of the figure shows the phylogenetic tree of vSAGs, and the right-side bar plot shows the completeness of each vSAG. The histogram at the top shows the number of VCs containing each AMG, showing 33 types of AMGs from multiple VCs and 37 types of AMGs from a single VC. (B) Proteomic tree created by ViPTree of vSAGs with BTLCP detected. Maximum-likelihood phylogenies of (C) TerL and (D) BTLCP. Branches with bootstrap support by SH-like aLRT higher than 80% and ultrafast bootstrap 95% are indicated by black dots. If those branches were only supported by SH-like aLRT higher than 80%, they are shown in gray. Each color of the outer circle represents a different VC value. Each scale bar means tree scale: 1.

Figure 5

Comparative genomics within the same vOTU (vOTU572). (A) An OGs matrix where each row represents a single OG and columns represent six vSAGs and one vMAG grouped into vOTU572. The top sequence is derived from the vMAG, and the others are derived from vSAGs. The presence or absence of a gene is indicated by color. Gray columns indicate genes that were not annotated, and white columns indicate the genes that were not detected. (B) Alignment of six vSAGs partial genomes. Open reading frames are represented by block arrows. Gray-colored arrows indicate genes in the core region, and other colored arrows indicate genes in the flexible region. Six different MTase subtypes were detected in vOTU572. The lead depth in the flexible region of vSAGs was 258.8 on average and 23.3 at the lowest.