Saline lakes on the Qinghai-Tibet Plateau harbor unique viral assemblages mediating microbial environmental adaption.

The highest plateau on Earth, Qinghai-Tibet Plateau, contains thousands of lakes with broad salinity and diverse and unique microbial communities. However, little is known about their co-occurring viruses. Herein, we identify 4,560 viral Operational Taxonomic Units (vOTUs) from six viromes of three saline lakes on Qinghai-Tibet Plateau, with less than 1% that could be classified. Most of the predicted vOTUs were associated with the dominant bacterial and archaeal phyla. Virus-encoded auxiliary metabolic genes suggest that viruses influence microbial metabolisms of carbon, nitrogen, sulfur, and lipid; the antibiotic resistance mediation; and their salinity adaption. The six viromes clustered together with the ice core viromes and bathypelagic ocean viromes and might represent a new viral habitat. This study has revealed the unique characteristics and potential ecological roles of DNA viromes in the lakes of the highest plateau and established a foundation for the recognition of the viral roles in plateau lake ecosystems.

Introduction

High-altitude plateau areas, which among the least-explored terrestrial biospheres on Earth, are characterized by persistent low temperatures and very sensitive to global warming. The Qinghai-Tibet Plateau is the highest plateau in the world, with an average altitude over 4,000 m. It is the primary water source for south and east Asia (Immerzeel et al., 2010) and contains more than 32,000 lakes, with a total area of 44,993 km2, which comprise the world's largest alpine lake group (Zhang et al., 2011, 2014). In addition, the Qinghai-Tibet Plateau contains 36,793 glaciers, with a total area and volume of 49,873 km2 and 4,561 km3, respectively (Immerzeel et al., 2010; Yao, 2010). In recent decades, under the influence of global climate change, the glaciers of the Qinghai-Tibet Plateau have started to melt and retreat. The melting glacier waters have then increased the total area of lakes in this area (Liu et al., 2019a). Most of these lakes are remote from human disturbance. Together, they have broad environmental gradients, including in salinity (from 0.1 to over 400‰), pH (5.4–10.2), and altitude (Yang et al., 2016; Zheng and Liu, 2009), and provide a unique, natural, and extreme laboratory to study microbial variation and adaption to environmental changes.

Recently, the prokaryotic and eukaryotic communities of soil, sand, lakes, and wetlands on high plateaux, such as the Qinghai-Tibet Plateau, Yunnan Plateau, and Colorado Plateau, have been investigated using high-throughput sequencing techniques (Lee et al., 2016b; Li et al., 2017; Liu et al., 2019a, 2016; Yang et al., 2019b; Zhang et al., 2015, 2020). The composition and function of the microorganisms in these high-altitude environments were found to be significantly different from those of low-altitude environments and were also different between lakes of different environments (Ji et al., 2019; Yang et al., 2016). Microbial diversity in the high-altitude lakes decreased significantly with dilution from melting glaciers (Liu et al., 2019a), whereas the network complexity seemed to increase in the hypersaline lakes (Ji et al., 2019). In addition, abundant microbes in wetland soils of Qinghai-Tibet Plateau have better potential adaptability to environmental change than rarer ones (Wan et al., 2021). However, information on their corresponding viromes in these high-altitude plateau lakes remains rare, especially over 3,000 m.

Viruses are the most abundant and diverse “biological entities” on Earth (Suttle, 2005). They significantly influence the microbial community structure, alter metabolic output through virus-encoded auxiliary metabolic genes (vAMGs), mediate horizontal gene transfer, drive global biogeochemical cycles, and contribute to carbon sequestration through viral lysing and aggregation (Suttle, 2005, 2007). During the past decade, the importance of viruses in global biogeochemical cycles and the composition of the viral community at a genome level has been explored using metagenomics and high-throughput sequencing, covering ocean, soil, and low-altitude lake environments (Brum and Sullivan, 2015; Chow and Suttle, 2015; Coutinho et al., 2020; Gregory et al., 2019b; Jin et al., 2019). Recently, 2,314,329 uncultivated viral genomes were identified with their relationship to habitat type from 21,075 public datasets in IMG/VR v3 (Roux et al., 2021). Most of the understanding of these viruses has been reported from marine environments. Our knowledge of high-altitude plateaux viromes remains limited (Eissler et al., 2020; Zhong et al., 2020). In particular, until now, there has been no information on viromes at a genome level from lakes on the Qinghai-Tibet Plateau.

Here we use the metagenomic method to examine the genetic composition of viruses in three Qinghai-Tibet Plateau high-altitude lakes, Qinghai Lake (QhL), Da Qaidam Lake (DQL), and Xiao Qaidam Lake (XQL), which span a salinity gradient from 13 to 246. This study aims to characterize the viral assemblages by assessing viral diversity, virus-host interaction, and virus-encoded auxiliary metabolic genes (vAMGs) and explore their adaption to different salinities.

Results

Overview of six viromes of three lakes on Qinghai-Tibet Plateau

Six viral samples were collected from the surface waters of three lakes, with distinct environmental parameters, on the Qinghai-Tibet Plateau; this includes four samples around QhL, one from XQL and one from DQL (Figure 1). All samples were alkaline with pH values from 7.8 to 9.1. Salinities ranged from 246 in DQL, 79 in XQL, and from 13 to 17‰ in QhL (Table S1). In DQL, sodium and chloride ions accounted for about three-quarters of the total dissolved solids.

Figure 1

Sampling stations and community structure of classified viral operational taxonomic units (vOTUs) of the six viromes of three lakes on Qinghai-Tibet Plateau

The pie charts show relative abundance on virus taxonomic composition, assessed for all vOTUs at the family level (except haloviruses) based on the results of CAT. The proportion of annotated vOTUs in all predicted vOTUs of each sample is shown below the names of stations.

Six viral metagenomic libraries of three lakes on Qinghai-Tibet Plateau were obtained after Illumina sequencing and contained about 72G bases raw data (a total of 239,539,878 raw paired-end reads [150bp each-end]). After quality control, these six viromes contained about 210 million paired clean reads. A total of 120,515 contigs longer than 1.5 kilobases were assembled. Combined with the prediction outputs of VirSorter, VirFinder, and CAT, a total of 4,835 viral contigs were detected (Gregory et al., 2019b; Ren et al., 2017; Roux et al., 2015; von Meijenfeldt et al., 2019). After redundant contigs were removed through alignment, 4,560 viral Operational Taxonomic Units (vOTUs) were classified (Table S2). Based on the results of VIBRANT, a fraction of the vOTUs were predicted to be temperate phages (Kieft et al., 2020). The proportion of lytic vOTUs in the four QhL samples (92.4%–98.2%) were lower than in XQL (98.7%) and DQL (99.2%, Figure S1), which was consistent with the salinity pattern.

Viral community structure

Shannon index values of the six viromes were calculated based on the relative abundance of vOTUs (Table S1); these varied from the lowest value of 4.03 at QhL1 to the highest of 6.88 at QhL10. However, the alpha diversities of the prokaryotic communities were different from those of the viruses. The sample from the hypersaline DQL had the lowest Shannon index value (0.80) and XQL sample had the highest (3.37) (Table S1, X. Li, Q. Liu, C. Zhang, X. Zhou, C. Gu, X. Yu, M. Wang, H. Shao, J. Li, Y. Jiang, unpublished data). The most abundant vOTUs differed between the three lakes (Figure S2) but the four viromes of QhL were more similar to each other than the other lakes.

At a family level, vOTUs were annotated using CAT, based on the lowest common ancestor (LCA) algorithm, and represented an average relative abundance of 10.9% (von Meijenfeldt et al., 2019). Siphoviridae (33.8%–74.0%) and Podoviridae (17.9%–55.9%) were the abundant families, although DQL mainly consisted of haloviruses (71%, Figure 1). This observation is in accordance with the dominance of Halobacteria in the prokaryotic community of DQL (Li Xianrong et al., unpublished data). In addition, Microviridae, a nonenveloped ssDNA virus, accounted for a moderate proportion of viruses (21.2% in QhL3, 1.0% in QhL10, and 11.8% in DQL; Figure 1) (Fard, 2016). Nucleocytoplasmic large DNA viruses (NCLDVs, mainly Mimiviridae and Phycodnaviridae) were detected in QhL3, QhL10, and XQL by CAT (Figure 1) and also detected in QhL10 and XQL by NCLDVs marker gene (major capsid protein, MCP) using hmmsearch (Figure S4B) (Mistry et al., 2013; Schulz et al., 2020; von Meijenfeldt et al., 2019). In addition, Lavidaviridae, virophages related to NCLDVs, was detected either by CAT (QhL3, QhL10, QhL11, and XQL) or by hmmsearch against virophage marker gene MCP (QhL1, QhL3, and QhL10) (Mistry et al., 2013; Paez-Espino et al., 2019; von Meijenfeldt et al., 2019). The relative abundances of NCLDVs and virophages were very low and less than 1% (Figure 1).

The representative sequences of vOTUs of six viromes were searched against NCBI RefSeq genomes of viruses using BLASTN to classify the vOTUs at species level. However, no hit met the threshold (≥95% average nucleotide identity [ANI], ≥ 85% alignment fraction, ≤ 1 × 10−3 e-value). In addition, the vOTUs were searched against the IMG/VR v3 dataset, which is the largest available public viral dataset consisting of over 2,314,329 uncultivated viral genomes viral sequences, and against the Global Ocean Viromes (GOV 2.0) dataset with over 190,000 viral populations, with the same threshold to test the distinctiveness of the six viromes (Gregory et al., 2019b; Roux et al., 2021). However, again there was no hit against GOV 2.0, and only 0.7% of the vOTUs were aligned with the IMG/VR v3 dataset, most of which were linked to vOTUs from nonmarine saline and alkaline, marine, and freshwater environments (Figure 2C).

Figure 2

Taxonomic diversity of the six viromes of three lakes on Qinghai-Tibet Plateau

(A) Gene-sharing network of viral sequence based on sequences (≥10k) of the six viromes and IMG/VR v3 (seawater, freshwater, nonmarine saline, and alpine and terrestrial habitats). Nodes represent viral genomes, and edges indicate similarity based on shared protein clusters.

(B) Venn diagram of shared viral clusters among the five different environmental habitats.

(C) Assignment of vOTUs of the six viromes of three lakes on Qinghai-Tibet Plateau against IMG/VR v3.

To investigate the relationship between vOTUs of six viromes and viral sequences from other habitats, 261 viral contigs longer than 10k were selected and compared with over 6,000 reference and predicted complete viral genomes of seawater, freshwater, nonmarine saline and alkaline, and terrestrial environment in IMG/VR v3 dataset using vContact 2 (Figure 2A) (Bin Jang et al., 2019). Seventy-four viral clusters (VCs, with a very high probability of the same genus in a VC) were detected in the six viromes, whereas 14 VCs (18.9%) were not found in the other four habitats of the IMG/VR v3 dataset (Figure 2B) (Roux et al., 2021).

Phylogenetic analysis of viral contigs and complete and near-complete viral genomes

A total of 24 complete and four near-complete viral metagenome-assembled genomes (MAGs) were assembled from the six viromes, according to the thresholds of the GOV dataset (Tara Oceans Coordinators et al., 2016). Eighteen complete viral genomes were shorter than 10 kbp. The marker gene of the ssDNA Microviridae, the MCP gene, was recognized in 13 complete viral genomes using DIAMOND BLASTP against all Microviridae MCP and major coat protein downloaded from NCBI (Buchfink et al., 2015). The phylogenetic tree was constructed based on 1,418 MCP gene sequences of Microviridae (and major coat protein for Bullavirinae) (Figure 3A). Most MCP sequences of Microviridae were divided into four superclades, and the 13 sequences from the six viromes were all classified into superclade one. This superclade included all reported viruses in the Bullavirinae subfamily and did not include any virus in the Gokushovirinae subfamily, suggesting the close taxonomic relationship between the new-found microviruses of the six viromes and Bullavirinae. These results extended our knowledge of the diversity of Microviridae.

Figure 3

Phylogenetic and comparative genomics-based synteny analysis of microviruses

(A) The phylogenetic tree of maximum likelihood was generated from the major capsid protein. Only the nodes with the bootstrap test value higher than 95 were shown. The star indicated the new-found microviruses detected in the six viromes of three lakes on Qinghai-Tibet Plateau.

(B) Comparative genomics-based synteny analysis of 13 microviruses from the viromes and one reference genome sequence and relative abundance of microviruses of the viromes.

To further understand the genomic similarity of superclade one of Microviridae, 24 microvirus genomes were chosen to compare the ANI using OrthoANI (Lee et al., 2016a). A heatmap of ANI showed that most genomes in the same clade were homologous (Figure S9). DQL_55, QhL10_5914, and QhL3_4981 were very similar with an ANI higher than 99%, which is the same for DQL_44, QhL10_5533, and QhL3_4669, suggesting that these genomes might represent two species of Microviridae. Therefore, three microviruses (DQL_55 group, DQL_44 group, and QhL3_4806) were classified into clade 1–2. In addition, because of their high ANI values, three sequences (DQL_40, QhL10_5371, and QhL3_4528) might represent the same species, and two other sequences (DQL_151 and QhL3_7004) could be an additional species. Furthermore, the virome from the XQL was different from those of the other two lakes, and no complete microvirus genome was detected. Alignment between the 13 complete genomes and the reference (Antarctic microvirus TYR_006_V_SP_13) is shown in Figure 3B (Sommers et al., 2019). All 14 genomes had three conserved genes, including MCP, DNA pilot protein, and replication initiation protein. The results of alignment were consistent with the ANI results (Figure S9).

To further investigate the diversity of the two main groups of dsDNA viruses, Caudovirales and NCLDVs, phylogenetic trees, based on the gene sequences of Caudovirales terminase large subunit and NCLDVs MCP, are shown in Figure S4. A total of 281 phage terminase sequences was detected from the six viromes of salt lakes. The phylogenetic tree shows that most of putative phage terminase in the high diversity were closely related to Siphoviridae and Myoviridae (Figure S4A). Also, many viral clades consisting of terminase have only been detected in these six viromes, which suggests that these might represent an endemic Qinghai-Tibet Plateau Caudovirales taxa.

By comparison, only 14 NCLDVs MCP sequences were detected in QhL and XQL and none in hypersaline DQL. The phylogenetic tree shows that most NCLDVs of the saline lake viromes were close to several algal viruses belonging to Phycodnaviridae, such as Micromonas sp. Rcc1109 virus mpv1, Phaeocystis globosa virus 12T, Bathycoccos sp. RCC1105 virus BpV1, and Yellowstone lake phycodnavirus 3 (Figure S4B).

Host prediction of viruses

After comparing sequence similarity, tRNA sequences, CRISPR spacers, and oligonucleotide frequencies between prokaryotic metagenome-assembled genomes and vOTUs, putative hosts were predicted for 358 vOTUs of the six viromes (Figure 4), accounting for a relative abundance of 13.5% (Li et al., 2021). Most of the predicted vOTUs had narrow host ranges, with only 11 potentially exhibiting a broader host range across several classes. Interestingly, 25 vOTUs in the DQL were linked to the dominant archaea, Halobacteria, but no vOTUs were predicted to link to archaeal host in the other two lakes. Predicted prokaryotic hosts spanned one archaeal and 17 bacterial class-level taxa, with Alphaproteobacteria (34 vOTUs, accounting for 35.5% relative abundance of predicted vOTUs) being the most frequently predicted (Figure 4). The most abundant vOTU of the six viromes (vOTU_1374) was also linked to Alphaproteobacteria, which was the dominant prokaryotic class in QhL. Most of the vOTUs were linked to Proteobacteria and Bacteroidetes, which were the dominant prokaryotes in QhL and XQL at phylum level.

Figure 4

Host prediction of viral operational taxonomic units (vOTUs) of the six viromes of three lakes on Qinghai-Tibet Plateau against prokaryotic metagenome-assembled genomes

Undirected virus-host prediction network with nodes that represent vOTUs (gray small nodes) and their putative hosts (big nodes), colored on the basis of taxonomy at the class level. Edges represent in silico virus-host linkages and are colored on the basis of host. If a node of vOTU is connected to different host nodes, it means this vOTU is predicted to link to different host classes, indicating potential broad host range.

Focusing on the 1,244 most abundant vOTUs with a cumulative relative abundance of 80% (Figure S2A), putative hosts were predicted for 123 with a relative abundance of 15.4%. Fifteen of 123 vOTUs were classified into Podoviridae and haloviruses. Most of the 123 vOTUs were predicted to link to Bacteroidetes (41), Gammaproteobacteria (23), and Halobacteria (23). In addition, vOTU_13 with a broad putative host range and linking to Alphaproteobacteria and Betaproteobacteria, coexisted in DQL and QhL. vOTU_285 with a putative host (Bacteroidetes) was detected in XQL and QhL.

Functional analysis of protein clusters

A total of 27,321 putative open reading frames (ORFs) in the six viromes of saline lakes were detected using metaProdigal (Hyatt et al., 2010). According to the NCBI COG functional categories of genes annotated by eggNOG-Mapper (Figure S5), genes with an unknown function were the most abundant annotated ORFs (Huerta-Cepas et al., 2017). In addition, genes dealing with information storage and processing and genes related to replication, recombination, and repair were also abundant categories.

To analyze the functional diversity of genes of viruses in the six viromes, all-verses-all BLAST was run by DIAMOND against all predicted genes (Buchfink et al., 2015; Zhao et al., 2013). This method detected reciprocal best-BLAST hits and received putative orthologs information, which was used to classify protein clusters (PCs). A total of 22,367 PCs were classified (Figure 5A). However, only 3.67% of PCs were shared among different lakes and no PC was detected in all three lakes, which suggests that there were considerable differences in virus functional genes between the different lakes. The diversity of viral PCs in the six viromes was the highest in the XQL virome, which is about 3.6 times higher than in the DQL virome (Figure 5A).

Figure 5

Functional protein clusters of viral contigs of the six viromes of three lakes on Qinghai-Tibet Plateau

(A) Venn graph of viral protein clusters of three different lakes.

(B) KEGG functional categories of different viromes of three saline lakes.

(C) Number of protein clusters of different categories of virus-encoded auxiliary metabolic genes.

To better understand the functional genes in the six viromes, the relative abundance of each functional gene in each sample detected by GhostKOALA, except the group “not included in pathway or brite”, is shown in Figure 5B (Kanehisa et al., 2016). In brief, only a small proportion of the PCs (less than 15%) were annotated as genes with classified functions. The three categories of “BRITE hierarchies,” “metabolism,” and “genetic information processing” were all abundant (Figure 5B). The category “BRITE hierarchies”, which contains noncoding RNAs and protein families of metabolism, genetic information processes, and signaling and cellular processes, was the most abundant category, except in QhL3.

Virus-encoded auxiliary metabolic genes

Several kinds of vAMGs were annotated after comparing with different databases; these mainly included osmoregulation, carbohydrate-active, nitrogen metabolism, sulfur metabolism, lipid metabolism, and antibiotic resistance genes (Figure 5C). Carbohydrate active and lipid metabolism genes were the most diverse vAMGs. For carbohydrate-active genes, no polysaccharide lyases or carbohydrate esterases were detected in the six viromes. Glycoside hydrolases, which hydrolyze or rearrange glycosidic bonds and glycosyltransferases to form glycosidic bonds, were more diverse in QhL than the other two lakes. Interestingly, three adjacent glycoside hydrolase genes were found in the viral contig DQL-84, suggesting the potential viral mediation of carbohydrate metabolism of the host cells. For nitrogen metabolism, three tandem glutamine synthetase genes were detected in the viral contig QhL10-122, suggesting the potential for viral mediation of glutamine metabolism in the host cells. Moreover, no nitrogen metabolism genes were found in the DQL virome and no sulfur metabolism genes were detected in the viromes from QhL. Furthermore, three associated sox genes, L-cysteine S-thiosulfotransferase (soxA), S-sulfosulfanyl-L-cysteine sulfohydrolase (soxB), and sulfur-oxidizing protein (soxZ), were detected in the viral contig XQL-17618, which shows its potential mediation ability involving sulfur oxidation in the host cells. Many kinds of lipid metabolism genes were detected, especially for fatty acid biosynthesis and glycerophospholipid metabolism (Figure 5C). These genes might contribute to the formation of the components of the host's cell membrane or the virus envelope.

Interestingly, some antibiotic resistance genes, relating to different mechanisms of resistance to antibiotics, were detected in XQL and QhL, such as antibiotic efflux, antibiotic inactivation, and antibiotic target alteration (Figure 5C). In addition, many viral osmoregulation genes were detected in the six viromes, which might be beneficial for osmotic pressure regulation in the host. However, the lower the salinity of the lake, the more kinds of glucan exporter associated genes were detected. None were detected in the hypersaline DQL virome, whereas one was found in the XQL viromes and three in QhL.

Ribosomal genes and potential CRISPR-like sequence

Twenty different ribosomal genes were encoded by viruses in QhL and XQL but none were detected in DQL (Figure S6). These genes, such as L1 and S1, are different parts of different ribosomal components, but they did not consist of an integrated component.

A putative CRISPR-like sequence was detected from one near-complete viral genome (XQL_135) (Figure S7A). In the predicted ORF20 of XQL_135, a CRISPR-like sequence consisting of nine similar repeats and eight spacers was detected. However, not only the repeats but also the spacers had no hits with the IMG/M dataset. Using BLASTP hits against NCBI nr database, XQL_135 was found to be most similar to Lake Baikal phage Baikal-20-5-C28. Both viral genomes had a low GC content and a similar gene order (Cabello-Yeves et al., 2017). There was no ortholog fragment between the ORF20 of XQL_135 and Baikal-20-5-C28, whereas some sequences before and after this ORF contained reciprocal orthologs with Baikal-20-5-C28. According to the tetranucleotide analysis, the region of CRISPR-like sequence in ORF20 of XQL_135 had lower tetranucleotide adaptability compared with the whole sequence (Figure S7B), whereas the corresponding region in Baikal-20-5-C28 had a higher self-adaptability (Figure S7C). XQL_135 contained four T4-like genes (terminase large subunit T4 headful, portal protein Gp20, RNA polymerase sigma factor Gp55, and DNA polymerase clamp loader subunit Gp62), suggesting its close relationship with genus “T4-like viruses” of the Myoviridae family. No tRNA sequences were detected in XQL_135, although six were found in Baikal-20-5-C28. In addition, five carbohydrate-associated vAMGs were detected in XQL_135, whereas vAMGs associated with the iron-sulfur cluster were detected in Baikal-20-5-C28 (Figure S7).

Relationship with environmental factors

Canonical correspondence analysis (CCA) was conducted to examine the relationship between the viral communities and environmental factors (Figure 6B). The results show that the four QhL viromes were clustered together, whereas viromes of the other two lakes were separate; this is similar to the result of PCoA analysis (Figure 6A). After forward selection, four environmental factors, elevation, latitude, and concentrations of magnesium and calcium, were found to be related to variability in the viral community structure (P < 0.05). The QhL viromes were mainly correlated with elevation, whereas the DQL virome correlated with magnesium concentration, and the XQL virome was mainly correlated with calcium concentration (Figure 6B).

Figure 6

Comparisons among the six viromes of three lakes on Qinghai-Tibet Plateau and other environmental viromes

(A) Principal co-ordinate analysis (PCoA) of the six viromes of six viromes.

(B) Canonical correspondence analysis (CCA) of the relationship between viral populations and environmental factors.

(C) Nonmetric multidimensional scaling (NMDS) analysis of the six viromes of three lakes on Qinghai-Tibet Plateau and other environmental viromes. According to habitats, viromes are divided into four groups, bathypelagic (BATHY), other marine (Epi and Meso), Qinghai-Tibetan Plateau (QTP), and other terrestrial (LAND). Group BATHY is inside Group QTP.

Pearson's correlation analysis was used to examine the relationships between viral community structure, environmental factors, and prokaryotic alpha diversity (Shannon Index) (Figure S3). Viral alpha diversity indices were not significantly correlated with any factor (P > 0.05), although the correlation coefficient between prokaryotic and viral alpha diversity was relatively high (R = 0.775).

Comparison between the six viromes and other environmental viromes

To compare the six viromes of three lakes on Qinghai-Tibet Plateau with other environmental viromes, a nonmetric multidimensional scaling (NMDS) analysis was conducted, which showed that the 40 virome samples could be divided into three groups (Figure 6C). The first group consisted of viromes from the Arctic, Antarctic, temperate, and tropical epi- (0–200 meters below surface, mbs) and mesopelagic (200–1000 mbs) zones in GOV 2.0 (Gregory et al., 2019b), whereas the second group contained nonmarine viromes and XQL viromes. The third group included the high-altitude aquatic viromes from the Qinghai-Tibet Plateau (four viromes from QhL and the viromes from DQL), two viromes from northwestern Qinghai-Tibet Plateau glaciers, and ten viromes from the bathypelagic zone (1000–4000 mbs) of the ocean (Gregory et al., 2019b; Zhong et al., 2020). In addition, the DQL virome is similar to a virome from the alpine, hypersaline Lake Salar de Uyuni in Bolivia, South America (Figure 6C) (Eissler et al., 2020; Ramos-Barbero et al., 2019).

Discussion

Although metagenomic analyses greatly expanded our understanding of viromes, our knowledge of high-altitude saline lake viral communities is still sparse. In this study, 4,560 vOTUs with 24 complete and four near-complete viral MAGs from six viromes of three saline lakes on the Qinghai-Tibet Plateau were detected. Almost all vOTUs of the six viromes were not cultured, based on a comparison of results with the NCBI viral RefSeq genomes dataset, GOV 2.0, and IMG/VR v3 dataset (Gregory et al., 2019b; Roux et al., 2021). These vOTUs contained many vAMGs associated with biogeochemical cycles, ARGs, and the osmotic pressure regulation of host cells. Multivariate analysis showed that the viromes from the Qinghai-Tibet Plateau were significantly different from viromes of land and epi- and mesopelagic oceans. Overall, the Qinghai-Tibet Plateau probably harbors many unique viruses that had not been found in other natural environments because of the unique alpine and saline-alkaline aquatic environments.

Qinghai-Tibet Plateau harbor unique alpine and saline-alkaline viral assemblages

No identified hit was found for the six viromes of three lakes on Qinghai-Tibet Plateau against the NCBI viral RefSeq genomes dataset and GOV 2.0. On average, only 0.7% vOTU sequences could be found in the IMG/VR v3 dataset, which included diverse habitats, such as freshwater, different marine zones, and nonmarine saline and alpine waters (Figure 2B) (Roux et al., 2021). These results suggest that these viromes were distinct from currently reported viromes from other environments and reflect the inadequate research and absence of datasets of viromes from high-altitude lakes, especially from the Qinghai-Tibet Plateau. The proportion of lysogenic vOTUs in the saline lakes was less than 7% in the six viromes (Figure S1) and decreased with increasing salinity, suggesting that the viruses inhabiting these alpine, saline-alkaline lakes prefer a lytic life strategy.

Only a few viral sequences of the six viromes could be assigned to a taxonomic rank, which is similar to previous viral metagenomic studies (Liang et al., 2019; Roux et al., 2016; Yang et al., 2019a). Caudovirales were abundant in QhL and XQL, which had low to moderate salinities (Ackermann, 1998). This distribution is similar to other environments, including seawater, freshwater lakes, soil, and marine sediments (Fancello et al., 2013; Jasna et al., 2018; Kim et al., 2017; Parmar et al., 2018; Segobola et al., 2018; Skvortsov et al., 2016; Yang et al., 2019a; Zhang et al., 2021a; Zheng et al., 2021). However, viruses of DQL, a hypersaline lake, were unique and dominated by haloviruses (Figure 1). Hypersaline environments are a reservoir of viruses with diverse shapes and have the highest reported abundance in water (up to 1010 mL−1) (Boujelben et al., 2012; Santos et al., 2012). However, only a limited number of tailed haloviruses have been isolated from hypersaline environments, which is not in accordance with the reports based on transmission electron microscope observations (Atanasova et al., 2015; Santos et al., 2012). This paradox suggests that previous reports may be incomplete and biased because only isolated haloviruses were analyzed in hypersaline viromes. Viruses within Siphoviridae, with long and noncontractile tails, were the most abundant viral family in the saline lakes (Figure 1) and many areas of seawater (Jasna et al., 2018; Suttle, 2005; Zhang et al., 2021a). Siphoviruses tend to infect host cells with a lysogenic life strategy and are not often the most dominant viral family in natural aquatic environments (Parmar et al., 2018; Suttle, 2005). In contrast, myoviruses, with contractile, long tails and a relatively broader host range, are more commonly observed as the dominant family in aquatic environment (Parmar et al., 2018; Suttle, 2005; Zhang et al., 2021a). However, Myoviridae was relatively rare in the six viromes (Figure 1); this suggests the unique environmental characteristics of lakes on the Qinghai-Tibet Plateau, such as the presence of a broad salinity gradient and high altitude, might shape the host communities. So the presence of viruses with a broad host range would not be conducive to precise control of the host community, which would be supported by the results of host's prediction of viruses (Figure 4) (Ji et al., 2019; Liu et al., 2016, 2019b; Wang et al., 2011; Yang et al., 2016).

Some NCLDVs were detected in XQL and QhL10, which had moderate salinities (Figure 1); most of these were similar to phycodnaviruses that infect microalgae (Figure S4B) (Schulz et al., 2020). This occurrence pattern is similar to marine viral communities (Yang et al., 2019a). The NCLDV supergroup is commonly abundant and diverse in oceans but can also be found in freshwater and associated with most major eukaryotic lineages; it might be the most important contributor to the control of primary production in hypersaline environments despite the low biomass (Schulz et al., 2020; Uritskiy et al., 2019). However, only a limited number of NCLDVs have been isolated from hypersaline environments (Atanasova et al., 2015). This study has expanded the number of NCLDVs found in hypersaline environments.

In addition, without multiple displacement amplification (MDA) some short putative ssDNA complete genomes were detected and 13 complete microvirus genomes were selected for phylogenetic and comparative genomic analysis. The ssDNA viruses are abundant in environments, such as sediments and the gut (Creasy et al., 2018; Gregory et al., 2019a; Yoshida et al., 2018; Yu et al., 2018). ssDNA viruses that have a small genome, only approximately 2–15 kb, may encode only a single structural protein and a single protein involved in their DNA replication (O'Carroll and Rein, 2016). The Microviridae contains two subfamilies, Gokushovirinae and Bullavirinae (Cherwa and Fane, 2011). Unlike dsDNA Caudovirales, no visible tails could be detected on microviruses (Cherwa and Fane, 2011). By comparison with MCP genes of microviruses in GenBank, four major superclades were clustered on the phylogenetic tree (Figure 3A). New microviruses in the six viromes and the isolated bullaviruses were classified into superclade one, without any gokushoviruses, suggesting a close relationship between these new-found microviruses and the Bullavirinae subfamily. Three annotated coding sequences, MCP, DNA pilot protein, and replication initiation protein, were detected in all selected microviruses (Figure 3B); this suggests that these microviruses probably had a common ancestor. Our study defined and expanded the diversity of Microviridae to a global scale.

From the host prediction of viruses, abundant vOTUs were linked to the dominant bacterial and archaeal populations. Three hundred fifty-eight of 4,560 vOTUs could be linked to putative hosts, which represents 13.5% of the viral abundance in the six viromes of three lakes. This abundance is lower than in other environments, such as cold, deep-sea seep sediments, suggesting the uniqueness of the viral and microbial communities in the Qinghai-Tibet Plateau lakes (Li et al., 2021). Although most predicted vOTUs were linked to Bacteroidetes (90), Gammaproteobacteria (73), and other bacteria (51), the most abundant vOTUs (35.5% in relative abundance) were linked to Alphaproteobacteria, including the most abundant vOTU (vOTU_1374). Alphaproteobacteria, which contains many important and abundant microbes, can utilize terrestrial-derived dissolved organic matter, especially that related to CHO and CHOS mixtures (Kang et al., 2013, p. 116; Yang et al., 2020; Zhao et al., 2013, p. 11). These results suggest that these abundant viruses linked to Alphaproteobacteria might have an important role in Qinghai-Tibetan Plateau lakes’ biogeochemical cycles.

In addition, 25 vOTUs linked to the dominant archaea, Halobacteria, were detected in DQL, the hypersaline lake (Figure 4). This observation is consistent with the community structure of viruses and prokaryotes in DQL (Figure 1, Li Xianrong et al., unpublished data), suggesting that DQL harbors special archaeal viruses.

Potential virus-mediated biogeochemical cycling, antibiotic resistance, and osmotic pressure regulation of host cells

In this study, 21 different categories of vAMGs related to carbon, nitrogen, sulfur, and lipid metabolism were detected. As viruses can carry and express some vAMGs to directly mediate the metabolism of the host cells, such as genes involved in photosynthesis, phosphate removal, and nitrogen metabolism, and indirectly affect biogeochemical cycles (Brum and Sullivan, 2015; Zimmerman et al., 2020), the presence of vAMGs suggests that viruses might mediate the metabolism of host cells and influence biogeochemical cycles in this area. Glycoside hydrolases, which hydrolyze or rearrange glycosidic bonds and glycosyltransferases to form glycosidic bonds, were more diverse in QhL than the other two lakes. The presence of these amino acids suggests that the viruses in the Qinghai-Tibet Plateau lakes might participate in carbon cycling through the mediation of carbohydrate synthesis and decomposition metabolism of the host cells by gene expression of vAMGs. In these lakes, viruses might mediate the carbon cycle, through lysis of the host cells, i.e. the “viral shunt,” increasing aggregation and vertical carbon flux to deeper waters by viral lysates, named the “viral shuttle,” and manipulation of the hosts' metabolic pathways by the expression of virus-encoded carbohydrate-associated genes (Sullivan et al., 2017; Suttle, 2007; Zimmerman et al., 2020).

Interestingly, viral antibiotic resistance genes (ARGs) were detected in the six viromes, indicating that viruses might help their host survive environmental antibiotics introduced by human activities (Moon et al., 2020; Yang et al., 2019b). Previous studies have shown that viruses can be key vehicles of ARG transfer in bacteria. This might be a significant threat to global health, as viral ARGs are common in urban water, oceans, and animal feces (Debroas and Siguret, 2019). ARGs have been observed in the wetlands across the Qinghai-Tibet Plateau but with lower abundance than in coastal estuaries and sites associated with human activity (Yang et al., 2019b). It is suggested that viruses with these genes might help the host cell survive the relatively high concentrations of antibiotics. The antibiotic concentrations in QhL (1.14–17.3 ng/L) were higher than in input rivers (Li et al., 2020; Yang et al., 2018). These results suggest that the viruses on the pristine Qinghai-Tibet Plateau are an important gene pool for ARGs, and so special attention should be paid to the usage of antibiotics in this area.

Viruses might also provide assistance to host cells for the osmotic pressure regulation, which is important for bacterial adaption to saline environments (Singh et al., 2013). Different virus-encoded osmoregulation genes were detected in saline lakes (Figure 6). These genes could help their host adapt to changing osmotic pressure in two different ways. The first approach, “salt-in-cytoplasm mechanism,” is to accumulate salt to match the environmental osmolarity (Kunte, 2006; Ma et al., 2010). The second approach, “organic-osmolyte mechanism,” requires highly water-soluble and uncharged organic compounds to maintain an osmotic equilibrium with the surrounding environment (Kunte, 2006; Ma et al., 2010). Interestingly, three and one vAMG encoding glucan exporter ATP-binding protein was detected in QhL, with low osmotic pressure, and in XQL, with moderate osmotic pressure, respectively, suggesting that viruses might help host cells to export glucan from the cytoplasm to periplasm to adapt to the low osmotic pressure (Lu et al., 2020).

Qinghai-Tibet Plateau harbors viral-encoded ribosomal protein genes and CRISPR-like sequence

Twenty virus-encoded elements of ribosomal proteins, which might control host transcription, were found (Figure S6). These genes, which have been detected in viral genomes, might replace some components of host ribosomes to moderate host translation and might give a functional fitness advantage during phage evolution (Mizuno et al., 2019).

One putative CRISPR-like sequence, consisting of nine similar repeats and eight spacers, was detected from the near-complete viral genome XQL_135 (ORF20), which might belong to the “T4-like virus” genus of the Myoviridae family; it was most similar to the Lake Baikal phage Baikal-20-5-C28 (Figure S7) (Cabello-Yeves et al., 2017; Chen et al., 2019). However, neither repeats nor spacers had hits with the IMG/M dataset and Baikal-20-5-C28, whereas some sequences before and after the ORF20 contained reciprocal orthologs with Baikal-20-5-C28. This indicates that the putative CRISPR-like sequence might be one of the inserted or removed parts of their common ancestor. According to the tetranucleotide analysis, the region of the CRISPR-like sequence in ORF20 of XQL_135 had lower tetranucleotide adaptability compared with the whole sequence, whereas the corresponding region in Baikal-20-5-C28 had a higher self-adaptability, suggesting that this region in XQL_135 might have been inserted later (Duhaime et al., 2011).

Qinghai-Tibet Plateau might represent a viral ecological zone

Marine viral communities are separated into five distinct viral ecological zones (VEZs), according to GOV 2.0 (Gregory et al., 2019b). However, few comparisons have been made between the five marine VEZs and terrestrial viromes, especially from alpine and saline-alkaline lakes. Here we compared the similarity of the six viromes, two viromes from glacier ice of the Qinghai-Tibet Plateau, eight other nonmarine viromes, and 24 viromes of GOV 2.0 (Adriaenssens et al., 2017; Coutinho et al., 2020; Eissler et al., 2020; Gregory et al., 2019b; Skvortsov et al., 2016; Zhong et al., 2020). According to the NMDS analysis and hierarchical clustering, the viromes could be divided into three groups. The first group consisted of viromes from the Arctic, Antarctic, temperate and tropical epipelagic zones, and mesopelagic zones in GOV 2.0, whereas the second group contained other nonmarine viromes and XQL viromes. However, the other seven Qinghai-Tibet Plateau viromes and the bathypelagic viromes clustered together forming the third group. As these represent the highest and lowest parts of the Earth, it is very interesting that viruses from these two areas are clustered together. Given the great environmental differences between the Qinghai-Tibet Plateau and the bathypelagic zone, it is speculated that viruses of these two environments might either have similar origins or possess adaption strategies to the cold and high-/low-pressure environments, and these two environments together might constitute a VEZ for viruses. However, the virome from XQL, which was assigned to another group, was closed to other QTP viromes (Figure 6C) and might represent a transitional state from the VEZ to other aquatic environments.

A study of marine viral communities has found that temperature is one of the most important environmental factors influencing viruses (Gregory et al., 2019b). The main factor influencing the planktonic viral communities in estuaries, however, was salinity (Zhang et al., 2021a). Four environmental factors, elevation, latitude, and concentrations of magnesium and calcium, were related to the variations of viral community structures. Elevation and latitude contributed to the variation of temperature, which might be an important factor affecting viral diversity. Alternatively, the concentrations of calcium and magnesium ions varied between lakes. Adsorption, as the first step of viral infection cycles, is influenced by divalent cations, and so calcium ions might inhibit the combination of phage QHHSV-1 with its host Halomonas ventosae QH52-2 (Fu et al., 2016). Thus, magnesium and calcium ions might further influence viruses through adsorption and therefore the beginning of the infection cycle (Shao and Wang, 2008).

Conclusions

Lakes on the Qinghai-Tibet Plateau harbor diverse and unique microbes, yet their associated viruses and influences on lake microbiomes have not yet been identified. In this study, for the first time, six viromes from three lakes with altitude above 3,000 m on the Qinghai-Tibet Plateau are revealed. The unique high-altitude plateau, which contains lakes spanning a broad salinity gradient, results in a unique and distinctive viral community and distribution. Almost all vOTUs of alpine saline lakes were different from those of other environments and most had predicted hosts that were linked to bacteria, whereas 25 linked to the dominant archaea, Halobacteria, which was only detected in the DQL, the hypersaline lake. Phylogenetic analysis expanded our knowledge of the Caudovirales, NCLDVs, and redefined the ssDNA Microviridae on a global scale. The life strategies of the viruses inhabiting the high-altitude plateau seem to indicate that they have adopted a lytic life cycle and potentially help their host adapt to osmotic pressure through the expression of virus-encoded osmoregulation genes. It also suggests that the six viromes of three lakes on Qinghai-Tibet Plateau might be delineated separately from land and epi- and mesopelagic ocean viromes and form a VEZ with bathypelagic ocean viromes. These data begin to fill the information gap on viruses from lakes on the highest Qinghai-Tibet Plateau and shed light on their potential influences on the associated microbiomes. Future studies expanding viral spatial-temporal dynamics of more lakes with deeper sequencing, covering a larger area and longer time span on the Qinghai-Tibet Plateau, will provide a better understanding of viral diversity and co-evolution and interactions with host cells.

Limitations of the study

In this study, only six viromes of three lakes were analyzed. Overall, little is known about viral communities in lakes on the whole plateau. In addition, viral communities in these lakes at different depths of different seasons are not known. Future studies about viral spatial-temporal dynamics of more lakes would help us understand viral diversity on Qinghai-Tibet Plateau comprehensively.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Min Wang (mingwang@ouc.edu.cn)

Materials availability

This study did not generate new unique reagents.

Method details

Experimental model and subject details

Six samples were collected form QhL (4 surface stations), XQL (1 surface) and DQL (1 surface), three saline lakes in the Qinghai Province of China, during May 2019 (Figure 1). The methods used for measuring environmental factors and prokaryotic alpha diversity based on 16s amplicons were as before (Xiong et al., 2012). Descriptions of all six samples is shown in Table S1.

In brief, six lake water samples, 50L each, were collected for the following study.

Virus concentration, DNA extraction and metagenomic sequencing

Lake water samples were collected and filtered through 3.0 μm and 0.2 μm polycarbonate membrane filters (Millipore, MA, USA). 50 L of the filtrate was added to 2.5 mL of 10 g/L Fe Stock (FeCl3) Solution. After shaking and settling, the mixed water was filtered through 0.8 μm polycarbonate membrane filters (Millipore, MA, USA). The 0.8 μm polycarbonate membrane filters were stored at – 80°C (John et al., 2011).

The viral particles on the filters were resuspended in fresh 0.1 M EDTA-0.2 M MgCl2 buffer (pH 6.0) at 4°C. Viruses were concentrated by centrifugal ultrafiltration into 400 μL. The QIAamp DNA Mini Kit (QIAGEN) was used to extract DNA from the concentrated solution following manufacturer handling instructions. Library construction and next generation sequencing of the viral DNA was carried out by Novogene (Nanjing, China) using Illumina NovaSeq 6000 (pair-end sequencing, 2 × 150 bp). In this study, the RNA component of the viromes was not considered as only DNA was extracted.

Quality control of reads, assembly and identification

The clean reads were picked from the raw reads by Novogene. Then, the high-quality paired-end reads were filtered using cutadapt and Perl scripts using the following conditions: (1) without N; (2) no more than 20% bases with a quality score less than 20; (3) no more than 30% bases with a quality score less than 30 (Martin, 2017; Yang et al., 2019a).

High-quality reads in each sample were assembled using metaSPAdes (v 3.12.0). Contigs which are less than 500 bp were removed (Nurk et al., 2017). MetaQUAST (v 5.0.2) was used to evaluate these metagenomic assemblies (Mikheenko et al., 2016). Then, contigs ≥ 1.5 kb were piped through VirFinder and VirSorter (v 1.0.6) and contigs that sorted with a VirFinder score ≥ 0.7 and p < 0.05 or VirSorter categories 1-6 were selected for further viral identification (Ren et al., 2017; Roux et al., 2015). The contigs with the following conditions were considered as viral contigs: (1) contigs sorted by VirFinder with a score ≥ 0.9 and p < 0.05; (2) contigs sorted by VirSorter in categories 1 and 2; (3) contigs sorted by VirFinder with a score ≥ 0.7 and p < 0.05 and VirSorter categories 1-6; (4) the remaining contigs were run through CAT (v 5.0.3) and contigs with CAT results of superkingdom not classified and the viruses were selected (Gregory et al., 2019b; von Meijenfeldt et al., 2019).

Dereplication, calculating relative abundances and taxonomic profiling

Viral contigs of each sample were piped through CD-hit-est (v 4.8.1) separately (Fu et al., 2012). Viral contigs were considered to be the same viral population if they shared ≥ 95% nucleotide identity across ≥ 80% of short viral contigs and the longest contig in each population was considered as a representative sequence. In order to calculate the relative abundances of viral populations in each sample, clean reads were mapped to representatives of the viral population using bowtie2 and counted using SAMtools (v 1.9) (Langmead and Salzberg, 2012; Li et al., 2009). Then, relative abundances of viral populations were calculated as transcripts per million (TPM) reads mapped (Li et al., 2010).

The drafts of taxonomic information of vOTUs were generated by CAT (von Meijenfeldt et al., 2019). As the CAT software found some NCLDVs and Lavidaviridae, all predicted genes of vOTUs were searched against NCLDV and virophage major capsid with Hidden Markov Models using hmmsearch (v 3.1b2, -E 1e-10) separately to test the results (Paez-Espino et al., 2019; Schulz et al., 2020).

VIBRANT (v 1.2.1) was used to detect the lifestyle of viruses with the default parameters (except -virome) to filter some contigs which were hard to detect (Kieft et al., 2020).

Identification of shared vOTUs

The viral populations in all six samples were piped through CD-hit-est in the same way as dereplication. Viral populations in different samples were considered to have shared vOTUs if they shared ≥ 95% nucleotide identity across ≥ 80%, while other viral populations were considered as individual vOTUs.

Virus–host prediction

Microbiome sequencing was performed by Shanghai Biozeron Biotechnology Co., Ltd. (Shanghai, China.) using 0.2 μm polycarbonate membrane filters (Li Xianrong et al., unpublished data). Reads from the quality control (by Trimmomatic v 0.36) were assembled using metaSPAdes (Bolger et al., 2014; Nurk et al., 2017). Contigs were binned (parameters: --metabat2 --maxbin2 --concoct) and refined (parameters: -c 70 -x 5) using metaWRAP (v 1.3.2) (Uritskiy et al., 2018). Taxonomy of each prokaryotic MAG was assigned by CAT (von Meijenfeldt et al., 2019).

Four different in silico methods were used to predict virus-host linkages, nucleotide sequence homology, oligonucleotide frequency, tRNA match and CRISPR spacer match (Li et al., 2021). (1) Sequences of vOTUs were compared with prokaryotic MAGs using BLASTN (≥70% minimum nucleotide identity, ≥75% coverage of viral sequences, ≥50 bit score and ≤0.001 e-value). (2) Sequences of vOTUs and prokaryotic MAGs were piped into VirHostMatcher (v1.0) to detect linkages (d2∗ values ≤0.2) (Ahlgren et al., 2017). (3) ARAGORN (v 1.2.38) was used to identify tRNAs from sequences of vOTUs and prokaryotic MAGs (parameters: -t) (Laslett, 2004). Each Match required ≥90% length with identity ≥90% of the sequences by reciprocal BLASTN (Coutinho et al., 2017; Zhao et al., 2013, p. 11). (4) CRISPR arrays were assembled from prokaryotic quality-controlled reads by crass (v 1.0.1) and CRISPR spacers were compared with viral contigs using BLASTN with ≤1 mismatch over the whole spacer (Li et al., 2021; Skennerton et al., 2013). For each match, BLASTN was used to compare the repeat and prokaryotic MAGs with the same parameters.

Comparison between the six viromes and other environmental viromes

To examine the relationship between six viromes of three lakes on Qinghai-Tibet Plateau and other environmental viral communities, 34 viral metagenomes from different environmental habitats were selected and download from NCBI Sequence Read Archive (SRA) data (https://www.ncbi.nlm.nih.gov/sra/) and CyVerse Data Commons (https://datacommons.cyverse.org/) (Table S3, Figure S8). Frequencies of di-, tri- and tetranucleotides, calculated by Python script.

The vOTUs were BLASTN (BLAST+ v 2.9.0) searched against the NCBI virus reference genomes with 95% identity, 85% of alignment length and an evalue of 0.001 were used to modify the drafts (Roux et al., 2019). In addition, vOTUs were blasted to the Integrated Microbial Genome/Virus (IMG/VR) system v.3 dataset with the same threshold in order to obtain further information (Figure 2C) (Roux et al., 2021). The vOTUs of the viromes were also BLASTN against GOV 2.0 contigs, over 10 k and circular, with the same threshold.

All vOTUs more than 10 k were piped into vConTACT 2 to with 113 reference genomes of freshwater, 3679 of seawater, 69 terrestrial viral reference genomes and 2187 complete, high quality viral sequences with IMG/VR v3 (Roux et al., 2021). All open reading frames (ORFs) were predicted by metaProdigal (v 2.6.3) (Hyatt et al., 2010). All proteins were compared using all-verses-all DIAMOND (v 0.9.29.130) BLASTP (e-value ≤1 × 10−5, query coverage ≥50%, identity ≥25%) (Buchfink et al., 2015; Zhang et al., 2021b). vConTACT 2 was then used to calculate the network (--db 'None' --pcs-mode MCL --vcs-mode ClusterONE) (Bin Jang et al., 2019). The network is illustrated using Gephi (v 0.9.2) (Figure 2A) (Bastian et al., 2009). The number of VCs were plotted using the VennDiagram in R (Figure 2B) (Chen and Boutros, 2011).

Open reading frames prediction and functional analyses

ORFs were recognized from vOTUs in each sample using metaProdigal with default parameters (Hyatt et al., 2010). To remove the contamination of microbial fragments, the non-viral ORFs in the contigs with provirus predicted by VIBRANT were removed (Kieft et al., 2020). Relative abundance of ORFs was calculated with the relative abundance of vOTUs.

ORFs were compared to the NCBI nonredundant protein database to evaluate the percentage of known genes using DIAMOND (v 0.9.29.130) BLASTP (Buchfink et al., 2015). In order to obtain the proportion of various functional genes, ORFs were piped through eggNOG V5.0 (http://eggnog5.embl.de/#/app/home) to identify COG functional classifications with default parameters, through PfamScan (e value < 1 × 10−5; bit score > 40) to identify domains and hidden Markov models (HMMs) and through KEGG GhostKOALA (https://www.kegg.jp/ghostkoala/) and KEGG KAAS (https://www.genome.jp/tools/kaas/) to predict KEGG Orthology pathways (Huerta-Cepas et al., 2017, 2019; Kanehisa et al., 2016; Mistry et al., 2013; Moriya et al., 2007; Tara Oceans Coordinators et al., 2016). ORFs were BLASTP against Carbohydrate-active enzymes (CAZy) database (http://www.cazy.org/) to find glycol-metabolism related genes using DIAMOND (e value < 1 × 10−5; > 60% identity with > 60% length of the query amino acid sequence) (Cao et al., 2020; Lombard et al., 2014). RGI (v 5.1.1) was run locally to detect antibiotic resistance genes excluding loose hits but nudging more than 95% identity loose hits to strict based on CARD reference data (Alcock et al., 2019). For further analysis of function, all ORFs were RPSBLAST against Conserved Domain Database (CDD) to annotate (Lu et al., 2020). All ORFs were all-verses-all BLAST to detect reciprocal BLAST hits (e value < 1 × 10−5; > 60% identity with > 60% length of the query amino acid sequence) (Cao et al., 2020; Zhao et al., 2013, p. 11). Query sequence and subject sequence of each reciprocal BLAST hit were considered to belong to the same protein cluster with the same function. Heatmaps of predicted ORFs function were plotted by TBtools (Chen et al., 2020).

Analysis of complete and near-complete viral genomes

All circular contigs detected by VirSorter and VIBRANT were viewed as complete viral genomes and all contigs more than 50 k with terminase detected by PfamScan and eggNOG-mapper were considered as near-complete viral genomes (Huerta-Cepas et al., 2017; Kieft et al., 2020; Roux et al., 2015; Tara Oceans Coordinators et al., 2016). These sequences were piped into Rapid Annotation using Subsystem Technology (RAST) (https://rast.nmpdr.org/) to annotate genome features (Overbeek et al., 2014).

MetaCRT (modified from CRT1.2) was used to detect CRISPR-like components of these high quality draft genomes (Rho et al., 2012). Only XQL-135 had a CRISPR-like sequence and not only repeaters but also spacers had no identical hit when BLASTN against RefSeq, CRISPR databases using CRISPRs web server and IMG/M (Chen et al., 2019; Grissa et al., 2007). All genes of XQL-135 predicted by prodigal were aligned against NCBI-nr database using DIAMOND to obtain the best hit of each gene (Buchfink et al., 2015). TBLASTX was run to detect homology parts plotted by Easyfig with Lake Baikal phage Baikal-20-5m-C28 (Cabello-Yeves et al., 2017; Sullivan et al., 2011). Tetranucleotide frequencies of two sequences were calculated by Python script with 10 kbp window size and 1 kbp step size and were normalized by z-scoring (Duhaime et al., 2011). The Pearson’s correlation coefficient was calculated by script from either the array of each fragment and the whole genome.

In order to obtain more credible predicted ORFs of shorter complete genomes, all complete genomes shorter than 10k were upload to GeneMarkS (http://topaz.gatech.edu/GeneMark/genemarks.cgi, sequence type: virus) (Besemer, 2001). These predicted protein sequences were aligned against all Microviridae MCPs and major coat proteins downloaded from NCBI using DIAMOND (e value < 1e-10) and two protein sequences at the beginning and end of the same genomes were overlapped if both of them had hits. 13 genomes of predicted microviruses were overlapped and adjusted to positive sequences starting at predicted MCP genes using the Sequence Manipulation Suite Version 2 (Stothard, 2000). All downloaded sequences dereplicated on the amino acid sequence combined with 13 new-found MCP genes were aligned with MAFFT (v 7.471) (auto mode) and trimmed using trimAL (v 1.4.rev15) to remove position with gaps over 90% (Capella-Gutierrez et al., 2009; Katoh and Standley, 2016; Schulz et al., 2020). The maximum-likelihood phylogenetic tree was computed by IQ-TREE (v 2.0.3) after choosing model (LG+F+R10) (Minh et al., 2020). The phylogenetic tree was visualized with iTOL v5 (https://itol.embl.de/) (Letunic and Bork, 2021).

24 genomes of microviruses were chosen to compare ANI using OrthoANI, which was plotted by TBtools (Chen et al., 2020; Lee et al., 2016a).

Phylogenetic analysis of Caudovirales and NCLDVs

All terminase sequences over 300 amino acids of Caudovirales were downloaded from NCBI virus and dereplicated. Terminase genes of the six viromes of three lakes on Qinghai-Tibet Plateau were detected by BLASTP using DIAMOND (-e 1 × 10−5) against downloaded sequences (Buchfink et al., 2015). Sequences with less than 300 amino acids were removed.

The MCP sequences of the six viromes were detected by hmmsearch (-E 1e-10) against MCP genes of NCLDVs dataset (Schulz et al., 2020). Sequences less than 100 amino acids were removed.

Dereplicated sequences combined with reference sequences were aligned with MAFFT (v 7.471) (auto mode) and trimmed using trimAL (v 1.4.rev15) to remove position with gaps over 90% (Capella-Gutierrez et al., 2009; Katoh and Standley, 2016; Schulz et al., 2020). The maximum-likelihood phylogenetic tree was computed by IQ-TREE (v 2.0.3) after choosing model (terminase: LG+F+R10; MCP: LG+F+R8) (Minh et al., 2020). The phylogenetic tree was visualized with iTOL v5 (Letunic and Bork, 2021).

Quantification and statistical analysis

The α- diversity (Shannon’s H) and β- diversity (Bray-Curtis dissimilarity) statistics were performed using vegan (v 2.5-6) in Rstudio (v 1.2.5001; R Stats v 3.6.1), based on the relative abundances of vOTUs (Dixon, 2003). PCoA based on Bray-Curtis dissimilarity matrices was conducted using cmdscale (k = 3, eig=T) in Rstudio. Figures were generated by ggplot2 (v 3.2.1) (Figure 6A). CCA was conducted by Canoco 5 with iterative selection of environmental factors (Figure 6B) (Jiangshan, 2014). The Pearson’s correlation between environmental factors and viral alpha diversity was calculated using IBM SPSS Statistics v 25.

Oligonucleotide frequencies of 40 viromes were piped into Rstudio and analyzed to compute euclidian distance and hierarchical cluster using vegan and pvclust (v 2.2-0) (method.hclust = “average”, method.dist = “cor”, nboot = 1000) libraries respectively (Gregory et al., 2019b; Roux et al., 2014; Suzuki and Shimodaira, 2006). The boundaries of groups were based on approximately unbiased (AU) p-value threshold 0.70 using function pvrect. The non-metric multidimensional scaling was based on frequencies using the metaMDS function (distance = “euclidean”, k = 2, trymax = 1000) and plotted using ggplot2 (Figure 6C).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Critical commercial assays
QIAamp DNA Mini Kit	QIAGEN	Cat. 51304
Deposited data
Global Ocean Viromes 2.0 raw reads	Gregory et al. (2019a, 2019b)	European Nucleotide Archive (ENA) and Integrated Microbial Genomes (IMG) with Joint Genome Institute - see Table S3 for details
Qinghai-Tibet Plateau lake viromes raw reads	This paper	Genome Sequence Archive: PRJCA005626
Tibet Glacier viromes QC reads	Zhong et al. (2020)	https://datacommons.cyverse.org/browse/iplant/home/shared/iVirus/Tibet_Glacier_viromes_2017/QC_Reads
Antarctic soil viromes raw reads	Adriaenssens et al. (2017)	Biosample: SAMN05767751, SAMN05767740
High-altitude wetland viromes raw reads	Eissler et al. (2020)	Bioproject: PRJNA613071
Lake Baikal virome raw reads	Coutinho et al. (2020)	Biosample: SAMN15248698
Salt lake Salar de Uyuni virome raw reads	Ramos-Barbero et al. (2019)	https://de.cyverse.org/dl/d/CEAB6764-643B-455C-BC1D-942B2B3544C5/ZN5.fa
Software and algorithms
cutadapt 2.4	Martin (2017)	https://github.com/marcelm/cutadapt
metaSPAdes 3.12.0	Nurk et al. (2017)	https://github.com/ablab/spades
MetaQUAST 5.0.2	Mikheenko et al. (2016)	http://bioinf.spbau.ru/metaquast
VirFinder	Ren et al. (2017)	https://github.com/jessieren/VirFinder
VirSorter 1.0.6	Roux et al. (2015)	https://github.com/simroux/VirSorter
CAT 5.0.3	von Meijenfeldt et al. (2019)	https://github.com/dutilh/CAT
prodigal 2.6.3	Hyatt et al. (2010)	https://github.com/hyattpd/Prodigal
CD-hit 4.8.1	Fu et al. (2012)	https://github.com/weizhongli/cdhit
bowtie2	Langmead and Salzberg (2012)	https://github.com/BenLangmead/bowtie2
SAMtools 1.9	Li et al. (2009)	https://github.com/samtools/
hmmer 3.1b2	Mistry et al. (2013)	http://www.hmmer.org/
VIBRANT 1.2.1	Kieft et al. (2020)	https://github.com/AnantharamanLab/VIBRANT
R 3.6.1	https://www.r-project.org/	https://www.r-project.org/
Rstudio 1.3.1093	https://www.rstudio.com/	https://www.rstudio.com/
Canoco 5	Jiangshan (2014)	http://www.canoco5.com/
IBM SPSS Statistics 25	https://www.ibm.com/products/spss-statistics	https://www.ibm.com/products/spss-statistics
Trimmomatic 0.36	Bolger et al. (2014)	https://github.com/usadellab/Trimmomatic
metaWRAP 1.3.2	Uritskiy et al. (2018)	https://github.com/bxlab/metaWRAP
BLAST 2.9.0+	ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/	ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/
VirHostMatcher 1.0	https://github.com/jessieren/VirHostMatcher	https://github.com/jessieren/VirHostMatcher
ARAGORN 1.2.38	Laslett (2004)	http://www.ansikte.se/ARAGORN/
crass 1.0.1	Skennerton et al. (2013)	https://github.com/ctSkennerton/crass
vConTACT 2	Bin Jang et al. (2019)	https://bitbucket.org/MAVERICLab/vcontact2
DIAMOND 0.9.29.130	Buchfink et al. (2015)	https://github.com/bbuchfink/diamond
Gephi 0.9.2	Bastian et al. (2009)	https://gephi.org/
eggNOG 5.0	Huerta-Cepas et al. (2019)	http://eggnog5.embl.de/#/app/home
PfamScan 1.6	ftp://ftp.ebi.ac.uk/pub/databases/Pfam/Tools/	ftp://ftp.ebi.ac.uk/pub/databases/Pfam/Tools/
KEGG GhostKOALA	Kanehisa et al. (2016)	https://www.kegg.jp/ghostkoala/
KEGG KAAS	Moriya et al. (2007)	https://www.genome.jp/tools/kaas/
RGI 5.1.1	Alcock et al. (2019)	https://github.com/arpcard/rgi
TBtools 1.09854	Chen et al. (2020)	https://github.com/CJ-Chen/TBtools
RAST	Overbeek et al. (2014)	https://rast.nmpdr.org
MetaCRT	Rho et al. (2012)	https://omics.informatics.indiana.edu/CRISPR/
Easyfig 2.2.5	Sullivan et al. (2011)	http://mjsull.github.io/Easyfig/
GenaMarkS	Besemer (2001)	http://topaz.gatech.edu/GeneMark/genemarks.cgi
Sequence Manipulation Suite Version 2	Stothard (2000)	http://www.bioinformatics.org/sms2/
MAFFT 7.471	Katoh and Standley (2016)	https://mafft.cbrc.jp/alignment/software/
trimAL 1.4.rev15	Capella-Gutierrez et al. (2009)	http://trimal.cgenomics.org/
IQ-TREE 2.0.3	Minh et al. (2020)	http://www.iqtree.org/
iTOL v5	Letunic and Bork (2021)	https://itol.embl.de/
OrthoANI[ 0.5.0	Lee et al. (2016a, 2016b)	https://github.com/althonos/orthoani
Vegan (R package)	Dixon (2003)	https://cran.r-project.org/web/packages/vegan/index.html
ggplot2 (R package)	https://cran.r-project.org/web/packages/ggplot2/index.html	https://cran.r-project.org/web/packages/ggplot2/index.html
pvclust (R package)	Suzuki and Shimodaira (2006)	https://cran.r-project.org/web/packages/pvclust/index.html
VennDiagram (R package)	https://cran.r-project.org/web/packages/VennDiagram/index.html	https://cran.r-project.org/web/packages/VennDiagram/index.html
Other
Illumina sequencing	Illumina	NovaSeq 6000