Comprehensive Insights Into Composition, Metabolic Potentials, and Interactions Among Archaeal, Bacterial, and Viral Assemblages in Meromictic Lake Shunet in Siberia

Microorganisms are critical to maintaining stratified biogeochemical characteristics in meromictic lakes; however, their community composition and potential roles in nutrient cycling are not thoroughly described. Both metagenomics and metaviromics were used to determine the composition and capacity of archaea, bacteria, and viruses along the water column in the landlocked meromictic Lake Shunet in Siberia. Deep sequencing of 265 Gb and high-quality assembly revealed a near-complete genome corresponding to Nonlabens sp. sh3vir. in a viral sample and 38 bacterial bins (0.2–5.3 Mb each). The mixolimnion (3.0 m) had the most diverse archaeal, bacterial, and viral communities, followed by the monimolimnion (5.5 m) and chemocline (5.0 m). The bacterial and archaeal communities were dominated by Thiocapsa and Methanococcoides, respectively, whereas the viral community was dominated by Siphoviridae. The archaeal and bacterial assemblages and the associated energy metabolism were significantly related to the various depths, in accordance with the stratification of physicochemical parameters. Reconstructed elemental nutrient cycles of the three layers were interconnected, including co-occurrence of denitrification and nitrogen fixation in each layer and involved unique processes due to specific biogeochemical properties at the respective depths. According to the gene annotation, several pre-dominant yet unknown and uncultured bacteria also play potentially important roles in nutrient cycling. Reciprocal BLAST analysis revealed that the viruses were specific to the host archaea and bacteria in the mixolimnion. This study provides insights into the bacterial, archaeal, and viral assemblages and the corresponding capacity potentials in Lake Shunet, one of the three meromictic lakes in central Asia. Lake Shunet was determined to harbor specific and diverse viral, bacterial, and archaeal communities that intimately interacted, revealing patterns shaped by indigenous physicochemical parameters.

Microorganisms are critical to maintaining stratified biogeochemical characteristics in meromictic lakes; however, their community composition and potential roles in nutrient cycling are not thoroughly described. Both metagenomics and metaviromics were used to determine the composition and capacity of archaea, bacteria, and viruses along the water column in the landlocked meromictic Lake Shunet in Siberia. Deep sequencing of 265 Gb and high-quality assembly revealed a near-complete genome corresponding to Nonlabens sp. sh3vir. in a viral sample and 38 bacterial bins (0.2-5.3 Mb each). The mixolimnion (3.0 m) had the most diverse archaeal, bacterial, and viral communities, followed by the monimolimnion (5.5 m) and chemocline (5.0 m). The bacterial and archaeal communities were dominated by Thiocapsa and Methanococcoides, respectively, whereas the viral community was dominated by Siphoviridae. The archaeal and bacterial assemblages and the associated energy metabolism were significantly related to the various depths, in accordance with the stratification of physicochemical parameters. Reconstructed elemental nutrient cycles of the three layers were interconnected, including co-occurrence of denitrification and nitrogen fixation in each layer and involved unique processes due to specific biogeochemical properties at the respective depths. According to the gene annotation, several pre-dominant yet unknown and uncultured bacteria also play potentially important roles in nutrient cycling. Reciprocal BLAST analysis revealed that the viruses were specific to the host archaea and bacteria in the mixolimnion. This study provides insights into

INTRODUCTION
Meromictic lakes have a mixed upper oxic mixolimnion, an interface chemocline, and a lower stagnant anoxic monimolimnion, which does not mix with either of the upper layers (Bowman et al., 2000). Their specific zones having stratified and unique biogeochemical characteristics provide novel opportunities to compare microbial loops and biogeochemical processes among the distinct zones (Bowman et al., 2000;Lauro et al., 2011;Comeau et al., 2012).
Lake Shunet (54 • 25 ′ N, 90 • 13 ′ E), located in the central Altai-Sayan mountain region in southern Siberia, is one of the only three meromictic lakes in the entire Asian region of Russia. This lake has sharp gradients of salinity, dissolved oxygen, hydrogen sulfide (H 2 S), and a dense population of anoxygenic phototrophic (purple sulfur) bacteria in the chemocline zone, which is comparable to that of Lake Mahoney (Canada), known for having the most purple sulfur bacteria of any lake in the world (Overmann et al., 1991;Rogozin et al., 2012). With its stable stratification, Lake Shunet provides a valuable opportunity for comparative studies on biogeochemistry (Kallistova et al., 2006), physicochemical parameters , and microbial communities (Lunina et al., 2007;Rogozin et al., 2009Rogozin et al., , 2010Rogozin et al., , 2012. Microscopic techniques, pigment analysis, culture, FISH-TSA techniques and PCR-DGGE have been used to disentangle the seasonal changes of anoxygenic photosynthetic bacterial community, sulfate reduction and methanogenesis and ecophysiology of phototrophic sulfur bacteria at the chemocline of Lake Shunet (Kallistova et al., 2006;Lunina et al., 2007;Rogozin et al., 2009Rogozin et al., , 2010Rogozin et al., , 2012. However, these approaches are limited in providing a fine-scale, comprehensive insight into the microbial assemblages and their potential functions. Therefore complete microbial community composition still remain unknown along the three zones. Until now, studies on the makeup of microbial populations and their functional properties in meromictic salt lakes by using whole-genome shotgun approaches have only been undertaken in Antarctica; that is, Ace Lake (Ng et al., 2010;Lauro et al., 2011) and Organic Lake . Similar to Ace Lake (68 • S) (Rankin et al., 1999), Lake Shunet (54 • N) is also a saline, high-latitude meromictic lake. Ace Lake is one of the most investigated lakes in Antarctica, in terms of its physicochemical profile, community structure, and functional dynamics (Laybourn-Parry and Bell, 2014). Both lakes are analogous, exhibiting characteristics such as permanent water stratification into three main zones, high H 2 S in the monimolimnion, and a gradient of dissolved oxygen ( Table 1) Laybourn-Parry and Bell, 2014).
Studies on marine viruses have demonstrated their contributions to the microbial community and the regulation of biogeochemical cycling (Suttle, 2007;Breitbart, 2012). Furthermore, the interaction between bacterioplankton populations and those of its phages have been uncovered in freshwater lakes and reservoir ecosystems (Roux et al., 2012;Skvortsov et al., 2016;Ghai et al., 2017), but still few studies have been employed on the ecologically important meromictic salt lakes. Meromictic salt lake ecosystems are known to have specific hydrochemical parameters and harbor particular microorganisms, providing excellent sedimentary archives of possible global paleoclimatic events (Dix et al., 1999), however they have received less attention than other aquatic ecosystems . Therefore, we chose to study Lake Shunet, a high-latitude meromictic salt lake in central Asia. The objectives of this study were to (1) explore the relationships between gene content or limnological parameters and microbial community structures over a depth profile; (2) reconstruct elemental nutrient cycles and investigate the potential role of the prevalent unclassified bacterial and archaeal groups in nutrient cycles; (3) unravel the history of potential interactions among viral, bacterial, and archaeal assemblages on the basis of the detection of clustered regularly interspaced short palindromic repeats (CRISPRs) and correlation of abundance profiles.

Geographic Description of Study Site and Sampling Procedures
Sampling was conducted in Lake Shunet (54 • 25 ′ N, 90 • 13 ′ E), which is 1.2 × 0.4 km in size (area, 0.47 km 2 ; depth, 6.2 m). This saline meromictic lake is located in the northern polar region in the Republic of Khakassia, Siberia  (Parnachev and Degermendzhy, 2002). It is permanently stratified, comprising an upper oxic layer (mixolimnion, 0.0-4.0 m), a transition layer (chemocline, 5.0 ± 0.2 m), and a lower anoxic layer (monimolimnion, 5.5-6.0 m) (Kallistova et al., 2006;Degermendzhy et al., 2010). All samples were collected on July 21, 2010. Approximately 20 L of water was pumped from pre-defined depths corresponding to each stratified layer (i.e., 3.0, 5.0, and 5.5 m), and the water was then collected in sterile containers for transportation. The entire sampling time was ∼2-2.5 h, including transportation of the samples back to the field station at Shira. The samples were filtered immediately at the station to reduce the potential bias in our samples involving transportation. These samples were subsequently filtered through a 10-µm plankton net and placed in a Millipore-Pellicon TFF system (0.22-µm filter membrane) to separate microbes and viruses (Tseng et al., 2013). Microbes were retained on 0.22-µm polycarbonate membrane filters, and viruses were collected from the filtrate water using a 50-kDa filter, which means viruses larger than 0.22 µm were excluded in the analysis. The viral samples were treated with DNase I (New England Biolabs UK Ltd., Herts, UK). The virus-like particles were purified in cesium chloride gradients by ultracentrifugation at a density of 1.35-1.50 g mL −1 (Angly et al., 2006;Tseng et al., 2013) for subsequent observation using TEM.

Analysis of Limnological Parameters
The vertical profiles of temperature, dissolved oxygen (O 2 ), and pH were measured using Hydrolab Data-Sonde 4a (Hydrolab, Austin, TX, U.S.A.) and YSI 6600 (Yellow Springs, OH, U.S.A.) submersible profilers. Conductivity readings at in situ temperatures (C t ) were standardized to specific conductance at 25 • C using the formula K 25 = C t × (1 + 0.0204 × (T-25)) −1 , where T is the in situ temperature in degrees Celsius (Hydrolab, YSI). Salinity was calculated on the basis of the relationship between salinity and conductivity, S (g L −1 ) = 1.117 K 25 -7.9716. The formula for conductivity and salinity (determined as ash content) was derived from water obtained from Lake Shira, another meromictic salt lake located 8 km northwest of Lake Shunet. Conductivity sensors were calibrated against 0.2 M KCl (Hydrolab, YSI) before each sample was collected. To determine sulfide concentrations, subsamples were fixed with zinc acetate and measured using the colorimetric method (Volkov and Zhabina, 1990). Total carbon (C tot ) and total nitrogen (N tot ) were determined by the CN-elemental analyzer FlashEA 1112 NC Soil/MAS 200 (Neolab LLC, U.S.A.). Ferric (Fe) and Sodium (Na) were detected using an inductively coupled argon plasma optical emission spectrometer. Phosphorus (P), NH + 4 , NO − 2 , NO − 3 , SO 2− 4 , and total minerals (Tot miner ) were all measured as described by Kalacheva et al. (2002). Limnological parameters were checked for collinearity using Spearman's rank correlation, and a set of non-collinear parameters was maintained for further analyses.

DNA Extraction, Archaeal/Bacterial Amplicon Libraries, and DNA Sequencing
The membrane filter with microbial retentate was rinsed with sterile double-distilled water and subsequently with 10-mL Milli-Q water. The filter was then transferred into a 50-mL Falcon tube for microbial DNA extraction using the cetyltrimethylammonium bromide method (Wilson, 2001). Viral DNAs were extracted from purified virus-like particles using the formamide/cetyltrimethylammonium bromide method (Angly et al., 2006) and inspected using agarose gel electrophoresis. To ensure that sufficient viral DNA was available for sequencing, viral DNAs were amplified using a Genomiphi kit (GE Healthcare Life Science, Piscataway, NJ, U.S.A.).
A portion of microbial DNAs was used to analyze the bacterial and archaeal community composition on the basis of the V4 and V1-V2 hypervariable regions of the 16S rRNA gene. These amplicon libraries were generated using the universal primer pair of the 16S rRNA gene, namely 571F/751R for archaea and 27F/341R for bacteria. The PCR reaction was carried out in a total volume of 50 µL, containing 2.5U Superrun EX taq TM HS, 5 µL of 10X EX taq Buffer, 200 µM of dNTPs, 0.2 µM of each primer, and 2-5 µg of diluted template DNA (final concentration 100 ng). The PCR program included an initial denaturation at 95 • C for 5 min, followed by 30 cycles of 95 • C for 20 s; 54 • C for 10 s for archaea or 52 • C for 20 s for bacteria; 72 • C for 20 s; and finally 72 • C for 5 min, with cooling at 4 • C. Thereafter, PCR products were separated by 1% agarose gel electrophoresis with 1X TE Buffer and SYBR R Green I. Bands of expected sizes (∼180 bp for archaea and 315 bp for bacteria) were cut from the gel and purified using a QIAEX II Gel Extraction Kit (QIAGEN Inc., Valencia, CA, U.S.A.). Purified DNA fragments were quantified using a NanoDrop spectrophotometer (Thermo Scientific, Vantaa, Finland). In the second round of the PCR process, individual tags were added to 5 ′ ends of the 571F/751R and 27F/341R primers for each sample. The PCR mixture contained 2.5U TaKaRa EX taq TM HS (TaKaRa Bio, Otsu, Japan), 5 µL of 10X EX taq Buffer, 200 µM of dNTPs, 0.4 µM of each tagged primer, and 100 ng V4/ V1-V2 amplicon in a final volume of 50 µL. The PCR program for tag addition had an initial denaturation at 95 • C for 5 min, followed by five cycles of 95 • C for 20 s, 56 • C for 10 s, and 72 • C for 20 s, with the final step of 72 • C for 5 min, and then cooling at 10 • C. The PCR products were purified, and a 200-ng mixture of tagged V4/V1-V2 regions was subject to 454 pyrosequencing using the Roche GS454 FLX Titanium System (Roche 454 Life Sciences, Branford, CT, U.S.A.) at Mission Biotech (Taipei, Taiwan).
After quality trimming of sequences, including chimera checking, and removal of ambiguous nucleotides (N), mismatched primers, and incomplete barcodes, 44,281/13,894 qualified bacterial/archaeal reads were retained. These reads were sorted into subgroups according to specific barcodes and then classified using a Ribosomal Database Project classifier (v2.3) with a bootstrap value of 0.8 for taxonomic assignment of sequences (Wang et al., 2007). Chloroplast reads were removed from subsequent analyses. Reads of each sample were aligned using MUSCLE (http://www.drive5.com/muscle). The distance matrix was calculated using PHYLIP package (v3.69) and clustered in MOTHUR (v.1.14.0) (http://www.mothur.org/) to assign OTUs at a threshold of 97% sequence similarity (Schloss et al., 2009). Notably, singletons were used only in estimating community richness and diversity; otherwise, they were excluded from analysis. In taxa names, NA means "not available at further levels." Comparison studies of Lake Shunet with Sakinaw (Gies et al., 2014), Ursu, and Fara Fund lakes (Andrei et al., 2015) were executed using the same bioinformatics strategy described in Gies et al. (2014); for example, taxonomic assignments for both archaeal and bacterial sequences were blasted against the SILVA database.

Sequencing, Assembly, and Annotation of Metagenomes
The purified microbial and viral DNAs (i.e., metagenomes) were sequenced on an Illumina HiSeq 2000 sequencing platform (San Diego, CA, U.S.A.) at Yourgene Bioscience (Taipei, Taiwan). Metagenomic reads with ambiguous nucleotides (>2) and short lengths (<35 bp) were removed and assembled through a de novo assembly algorithm within the CLC Genomics Workbench using a 40-bp minimum overlap and 99% consensus. ORFs and 16S rRNA were predicted on assembled contigs using MetaGeneMark (http://exon.gatech.edu) and RNAmmer (v1.2) (Lagesen et al., 2007). Taxonomic assignment of predicted 16S rRNA was executed using the Ribosomal Database Project classifier (v2.3) with a bootstrap value of 0.8. ORFs were annotated by searching against the EggNOG (http://eggnog. embl.de) and Kyoto Encyclopedia of Genes and Genomes protein (Kanehisa et al., 2012) databases using BLASTp (e-value ≤ 10 −5 and bit score ≥ 100). The abundance of ORFs was initially calculated on the basis of the coverage of read counts and then normalized by total sample size as relative abundance for further analyses. Data were log(x+1) transformed, and COG functional annotation was evaluated using STAMP (Parks and Beiko, 2010).

Comparative Metagenomic Analyses
For comparative analysis, metagenomic reads of Ace Lake (filtered by 0.8 and 3.0 µm) were assembled using the same assembly workflow (https://www.imicrobe.us/) with default settings. ORF prediction and functional annotation were then performed using the same strategy applied to the Lake Shunet metagenomes (see Sequencing, assembly, and annotation of metagenomes). The ORF abundance was normalized by total sample size as relative abundance for the analysis. NMDS with Bray-Curtis distance metrics was used to explore the similarity in COG profiles between Ace Lake and Lake Shunet. Statistical analyses were performed using the Vegan package (Oksanen et al., 2011) in R (https://www.r-project.org/).

Metagenomic Binning and Reconstruction of Carbon, Nitrogen, and Sulfur Cycles
Metagenomic contigs (>1,000 bp) were binned using a twotiered binning algorithm (Saeed et al., 2012). The binning method used has previously been validated on both simulated and realworld data, similar to the data sets explored in this study (Saeed et al., 2012). Each bin was subjected to a conservative BLAST/MEGAN analysis to confirm that bin was meaningful. The binning method takes advantage of unsupervised machine learning to group related sequences into phylogenetic clusters (a technique which relies less on training data, which can skew results, and more on extracting common patterns in the data itself). The two-tiered nature of the algorithm (a course grouping, followed by a more refined grouping) was demonstrated to group sequences at a higher resolution and with better accuracy than other comparative methods. Moreover, the unsupervised characteristics of the binning method make it highly aligned with metagenomic studies, such as this one, where a complete set of a priori information about the underlying communities is not always available. Therefore, the selected binning method was an ideal candidate to aid the analysis of raw sequence data, and that the results obtained using the binning method are representative of the underlying microbial communities that we expected to find. However, it is worthy to note that the taxonomy assignment result would have some variation using different binning methods. ORFs were predicted using the contigs contained within each microbial bin, and flavobacteria from the viral samples were annotated using a BLASTp search against the KEGG protein database (e-value ≤ 10 −5 ) to reconstruct carbon (C), nitrogen (N), and sulfur (S) cycles in Lake Shunet. The relative abundance of each bin was initially calculated according to the read counts of coverage of ORFs and subsequently normalized by total sample size.

Bioinformatics Analyses of Viral Communities and Measurement of Diameter of Virus-Like Particles
Viral richness and the Shannon-Wiener Index (H ′ ) were estimated using PHACCS (v1.1.3, http://sourceforce.net/ projects/phaccs), on the basis of contig spectra generated by Circonspect (v0.2.5, http://sourceforge.net/projects/circonspect) in combination with the Cap3 assembler. Viral community composition was defined according to a homology search. Viral ORFs were searched against the NCBI RefSeq viral and microbial protein collection (ftp://ftp.ncbi.nlm.nih.gov/refseq/release) using BLASTp (e-value ≤ 10 −5 , single best match) to infer taxonomic information. Only 1.06, 0.88, and 0.82% of ORFs in the respective viral metagenomes were assigned to viruses ( Figure S1). The final taxonomy of each contig was determined through the voting of constituted ORFs. However, many viral ORFs were not found to have any corresponding sequences in the databases; the composition prediction might be biased due to limited information. The relative abundance of each viral taxon was measured as described by Tseng et al. (2013). The function of viral ORFs was annotated using GhostKOALA of the KEGG database.
The diameter of virus-like particles was measured on the basis of TEM images using NIS-Elements AR 4.0 of Nikon microscopy and was calculated as follows: Diameter of virus − like particle = nm scale bar pixels scale bar × pixels of virus − like particle

Analysis of Flavobacterial Genomes in Viral Samples
Several contigs derived from viral metagenomes (the largest being ∼4.1 Mb) had hits against flavobacteria with a cutoff evalue of 10 −5 and bit score of ≥100. Analyses involving 16S rRNA genes, reciprocal BLAST, and genome completeness confirmed that these contigs belonged to flavobacteria. Furthermore, a reciprocal BLASTp was used to examine the similarity of flavobacteria from microbial and viral samples, with a cutoff e-value of 10 −5 and bit score of ≥100. Hidden Markov models using TIGRFAM/Pfam libraries (Haft et al., 2003;Finn et al., 2010) and CheckM (Parks et al., 2015) were used to estimate genome completeness. The fragment recruitment analysis was conducted using MUMmer (http://mummer. ourceforge.net/manual/) with default settings. Moreover, rRNA genes were predicted with WebMGA (Huang et al., 2009), and BLASTn searches against the NCBI nucleotide database (nt) were performed to identify taxonomically related species. According to BLASTn results, Nonlabens dokdonensis DSW-6 was the best hit (97% similarity), followed by Psychroflexus torquis ATCC 700755 (92% similarity) and Flavobacteria bacterium BBFL7 (92% similarity). The CLC Genomics Workbench (CLC bio) software was used to map contigs assigned to flavobacteria onto a draft genome. Contigs measuring >10 kb were used for single-copy gene analysis to estimate genome completeness. A genome map was drawn with Circos (Krzywinski et al., 2009) and annotated according to reference genomes (Nonlabens dokdonensis DSW-6, Psychroflexus torquis ATCC 700755, and Flavobacteria bacterium BBFL7).

Analyses of Interactions Between Microbial and Viral Assemblages
Appropriate computational approaches can predict bacteriophage-host relationships (Edwards et al., 2016). We used signal categories such as abundance profiles, CRISPRs, and exact matches. Congruence between ordinations of viral and bacterial/archaeal communities was quantified using the procrustean superimposition method (Procrustes) and examined for statistical significance using a Monte Carlo procedure (999 permutations). The m 12 values, representing a goodness-of-fit statistic that measures congruence between two ordination configurations, were transformed to Procrustes correlation (r) by calculating the square root of their complements using the Vegan package (Oksanen et al., 2011). Pearson correlation coefficients between viral and bacterial/archaeal diversity indices and the richness estimate were calculated using the Psych package (Revelle, 2011). CRISPRs were predicted by PILER-CR (v1.06, http://www.drive5.com/pilercr) with default settings. To predict potential interactions between hosts and viruses, reciprocal BLAST analysis between microbial and viral metagenomics was initially conducted, followed by ORF clustering, using CD-HIT with a threshold of 90% sequence similarity (Li and Godzik, 2006). Moreover, contigs of microbial bins were searched against viral contigs to infer potential host-virus interactions (BLASTn, e-value ≤ 8 × 10 −21 , bits ≥ 100).

Statistical Analyses
The following statistical analyses were undertaken with R software (Team, 2011). Functional diversity was calculated using the R package ShotgunFunctionalizeR (Kristiansson et al., 2009). Count data (including community structure and functional profiles) were transformed into relative abundance by dividing by sample size. Limnological parameters were log(x+1) transformed for all analyses. The clustering heat map was generated using a gplots package (Bolker et al., 2012). Differences among the three depths in terms of contents of genes involved in energy metabolism were examined using the Wilcoxon test (Bauer, 1972), whereas Student's t-test was used for testing differences in viral diversity, richness estimates, and viral-like particle sizes. Non-metric Multidimensional Scaling (NMDS) based on Bray-Curtis distances was performed using the Vegan package to visualize the distribution of the bacterial/archaeal community composition or the 50 most abundant functional genes involved in metabolism, followed by post-hoc regression of individual explanatory variables on the ordination scores. Goodness-of-fit values and their significance were calculated using 999 random permutations. Correspondence analysis (Perrière and Thioulouse, 2003) using the Vegan package was applied to compare the variances of functional genes involved in energy metabolism and the 10 most abundant unclassified bacterial/archaeal OTUs. Only the top 50 genes based on the derived relative abundance were shown on the ordination plot.

Phylogenetic Analysis
Phylogenetic trees of the functional gene psbA were constructed using the neighbor-joining method from 500 bootstrap iterations using MEGA 6.0 (Tamura et al., 2013). In addition to the four partial psbA genes obtained from virome data in this study, we downloaded genes of various Synechococcus and Prochlorococcus phages from the UniProt and NCBI databases. Phylogenetic trees were constructed separately, because of the lack of common regions among the four partial psbA genes.

Physicochemical Profile of Lake Water
The gradient of each physicochemical parameter varied according to depth (Table 1). Both pH (8.1-6.7) and temperature (15.5-7.5 • C) gradually decreased with depth. Furthermore, O 2 , H 2 S, and salinity changed sharply with depth; specifically, the O 2 level was reduced from 10.9 mg L −1 to nearly zero in the chemocline, whereas the salinity (∼71 g L −1 ) and H 2 S (∼400 mg L −1 ) levels were much higher in the monimolimnion (i.e., bottom water).

Specific Bacterial and Archaeal Community Structure and High Diversity Along the Depths in Lake Shunet
Three representative depths for the mixolimnion (3.0 m), chemocline (5.0 m), and monimolimnion (5.5 m) layers were selected for a microbial community survey, according to hydrographical data. Archaeal and bacterial 16S ribosomal RNA genes were subjected to massively parallel pyrosequencing to characterize the community composition.
After tagging and quality trimming of sequences, qualified sequences were submitted to the RDA pipeline to define the operational taxonomic unit (OTU) at 97% nucleotide identity. Following the removal of singletons, a total of 293 archaeal and 125 bacterial OTUs were identified. According to Shannon's diversity index (H ′ ), the 5.5-and 3.0-m depths had the most diverse bacteria (H ′ = 5.46) and archaea (H ′ = 6.24), respectively. Rarefaction curves showed that only the curves for samples obtained at the 5.0-m depth approached an asymptote, suggesting that additional sequencing efforts at the other two depths would detect greater diversity.
The archaeal and bacterial community structures comprised 3 phyla, 5 classes, and 5 genera and 9 phyla, 15 classes, and 26 genera (relative abundance > 1%), respectively, with the corresponding prevailing phyla being Euryarchaeota and Proteobacteria. Most archaeal sequences were affiliated with a Euryarchaeota in both the upper and lower layers and with Methanomicrobia in the chemocline. The bacterial assemblage of the upper layer comprised pre-dominantly Cyanobacteria, whereas Gammaproteobacteria were abundant in the other lower layers. At the genus level, an unknown genus of Euryarchaeota and Thiocapsa were the most abundant archaea and bacteria, respectively (Figures 1A,B). Every depth harbored unique archaea and bacteria: Nitrososphaera and Methanospirillum were only present specifically in the mixolimnion (Figure 1A), and specific depth-dependent bacterial genera were also identified, including Loktanella in the mixolimnion, Enterobacter in the chemocline, and Candidatus Cloacamonas in the monimolimnion ( Figure 1B). Notably, pre-dominant unclassified bacteria and archaea, as well as those at unknown genus levels, were observed across the three depths (Figures 1A,B), which echoed the relative abundance of genomic bins, revealing that unknown bacteria were abundant specifically in both upper and lower layers ( Table 2).
Bacterial assemblages were dominated by Cyanobacteria in the mixolimnion and Gammaproteobacteria in both the chemocline and monimolimnion in Lake Shunet. Shunet differs from Sakinaw Lake, in which Actinobacteria, Gammaproteobacteria, and Chloroflexi are pre-dominant along the layers (Gies et al., 2014); Ursu and Fara Fund Lakes, in which Gammaproteobacteria or Deltaproteobacteria pre-dominate along these layers (Andrei et al., 2015). The bacterial community compositions of these four meromictic lakes, each have their own dominant genera. Lake Shunet, for example, has dominant bacteria that the other lakes do not, including Thiocapsa, Candidate division OD1, and Halomonas ( Figure 1C). In terms of archaeal compositions, Lake Shunet and Sakinaw Lake were more similar to each other than the Ursu and Fara Fund Lakes; both groups of lakes harbored common archaeal groups ( Figure 1D). These results reflect the stratified nature and uniqueness of meromictic lakes.

Functional Profiles Varied Among Depths
Using a deep sequencing strategy to uncover and comprehend potential metabolisms of microbial communities in Lake Shunet, we generated sequences with sizes ranging between 26 and 72 Gb for each of the six samples (Table S1), as well as microbial and viral sequences with total approximate sizes of 132 and 133 Gb, respectively. Consequently, total bases of contigs were between 121 Mbp and 14 Gbp in length for the six metagenomes. The largest microbial and viral contigs were 0.47 and 0.09 Mbp, respectively (Table S1). A total of 1,302,473 microbial and 1,379,639 viral open reading frames (ORFs) were predicted from the contigs and were annotated in 4356 COGs families, of which 71.9% were shared by the three depths, and 9.1, 3.7, and 2.5% were specific to the respective depths. A functional diversity index was calculated on the basis of metagenomic data by using ShotgunFunctionalize R; the index of Lake Shunet (Shannon index = 7.67-7.81) was higher than that of the North Pacific Ocean (Shannon index = 6.46-6.93) (DeLong et al., 2006;Kristiansson et al., 2009).
The three layers harbored significantly different contents of genes involved in energy metabolism (P < 0.05; Figure 2). In terms of COG functional profiles, most of the genes were annotated as playing roles in replication, recombination and repair, amino acid transport and metabolism, cell wall/membrane/envelope biogenesis, and energy production and conversion. STAMP analysis conducted on the basis of the COG profiles revealed that the observed microbial metabolism, including energy production and conversion, significantly differed (P < 0.05) among the three depths ( Figure S2).

Relationships Between Metabolic Functions and Archaeal/Bacterial Community Structure and Limnological Parameters
Nonmetric multidimensional scaling (NMDS), based on Bray-Curtis distances with post-hoc regression, was used FIGURE 1 | Community structure of (A) archaea and (B) bacteria based on 16S rRNA amplicons (>1% relative abundance) at the genus level across the three depths of Lake Shunet. Heat map and functional clustering based on relative abundance (>1%) of (C) archaeal and (D) bacterial assemblages (against SILVIA database) of meromictic Shunet, Sakinaw, Ursu, and Fara Fund lakes along the depths at the genus level. sh, Lake Shunet; sak, Sakinaw Lake; ursu, Ursu Lake; farafund, Fara Fund Lake. Bacterial community compositions of the four meromictic lakes respectively cluster while Lake Shunet and Sakinaw Lake harbored common archaeal groups. to explore the correlation between the 50 most abundant genes involved in KEGG orthology (KO) profiles of metabolism and bacterial/archaeal assemblages (>0.1% relative abundance). Five archaeal groups, including Methanococcoides, Methanosphaera, and Methanolobus, and bacterial groups, including Desulfobacterium, Desulfuromonadaceae, and Thiocapsa, had a significant goodness of fit (P ≤ 0.05) regarding the relative importance of the 50 most abundant genes involved in metabolism (Figures 3A,B). Specifically, Methanococcoides were correlated with K00163 (pyruvate dehydrogenase E1 component), K00366 (ferredoxin-nitrite reductase), K00266 (glutamate synthase (NADPH/NADH) small chain), and K03891 (ubiquinol-cytochrome c reductase cytochrome b subunit), whereas K10203 (elongation of very long chain fatty acids protein 6) was enriched by Thiocapsa. We conducted NMDS on the basis of bacterial and archaeal community structures, followed by post-hoc regression of individual environmental variables; the results revealed that parameters including H 2 S, O 2 , NH + 4 , and NO − 2 , as well as 11 other parameters, significantly explained the variance of both community structures with depth (Figures 3C,D). Specifically, Methanococcoides and Methanolobus were positively associated with increases in total nitrogen (N tot ) and temperature, respectively ( Figure 3C). The bacterial clusters including Desulfonatronum, Halothiobacillus, Desulfobacteraceae, Desulfobacterium, and Candidatus_Cloacamonas were enriched as H 2 S increased, whereas Cyanobacteria, Marinobacter, Alcaligenaceae:NA, and Planctomycetaceae:NA were correlated with changes in O 2 (Figure 3D).

Reconstruction of Nutrient Cycles Through the Depths
Many long-sized contigs obtained by deep sequencing considerably facilitated binning metagenomic data, thus enabling the prediction of potential elemental nutrient networks among the microbial groups. A total of 14, 6, and 18 microbial bins (namely taxonomically related groups) were identified in the respective metagenomes of the three depths ( Table 2). The elemental nutrient cycles of the three layers were reconstructed according to the gene annotation in individual microbial genomic bins (Figure 4). Serial numbers were assigned to each bin in the form of "depth-number, " as shown in the following descriptions.
The constructed C, N, and S cycles annotated on the basis of the 3.0-m metagenome were also annotated in 3.0-m genomic bins; however, a few pathways involved in N cycles annotated on the basis of the 5.0 and 5.5-m metagenomes were not found in any genomic bins at the depths, such as nitrogen fixation, denitrification, and nitrate reduction. This means that some microbes playing a role in N cycles were still not defined.

Virus-Like Particle Size, Viral Genome Size, Diversity, and Function Annotation
The sizes of virus-like particles were measured using transmission electron microscopy (TEM; Figure 5A); the results revealed that the sizes decreased significantly across the depths, with average sizes of ∼ 32, 30, and 28 nm at the three sampling depths, respectively (P < 0.05; Figure 5B). A similar pattern was also observed in viral contig length, which was longer (28.9 kb on average) at 3.0 m than at 5.0 and 5.5 m (24.4 and 19.9 kbp, respectively).
The viral community structure differed among the three sampling depths (Figure 6). Fourteen abundant viral families (>0.1% relative abundance) were detected through the water column, including several typical bacteriophage groups, such as Siphoviridae, Microviridae, Myoviridae, and Podoviridae, and three viral groups (Nanoviridae, Circoviridae, and Phycodnaviridae). The viral assemblage at 3.0 m was dominated by Siphoviridae (35.6%), followed by Myoviridae, whereas the viral assemblages at 5.0 and 5.5 m were dominated by Microviridae (34.4%) and Nanoviridae (27.4%), followed by Nanoviridae and Siphoviridae, respectively. Furthermore, many unknown viruses were detected at the three depths, with the greatest enrichment (24.4%) being observed at 5.5 m (Figure 6). Viral diversity determined using PHAACS differed significantly among the three depths ( Table 3). The 3.0-m layer had the most genotype-rich (15,877 genotypes) and diverse (H ′ = 9.15) viral assemblage, followed by the 5.5-m layer (10,000 genotypes and H ′ = 8.06) and 5.0-m layer (2,498 genotypes and H ′ = 6.70). The highest evenness estimate was derived at 3.0 m, indicating that all viral genotypes were close to being equally abundant, whereas some viral genotypes at 5.0 and 5.5 m were pre-dominant. The lowest evenness was observed in the viral community at 5.0 m, and the most abundant genotype constituted 5.51% of the community.
In total, 9, 22, and 13 unique viral sequences were found in microbial bin sequences at the respective depths ( Table S2). The detected viral taxa, including Cronobacter phages Vb_CsaM_GAP31, Pseudomonas phage YuA, and Prochlorococcus phage P-SSM7, did not show host specificity to any bacteria. The microbial groups observed at the 3.0-m depth seemed to be attacked only by viruses obtained from this depth, except for Flavobacteria (3-2), which were attacked by Cronobacter phage obtained from the 5.0-m depth. However, the microbial groups observed at the 5.0-and 5.5-m depths seemed to be infected randomly by viruses from both depths.

Flavobacteria in Viral Samples at 3.0 m
Flavobacteria in the viral samples at 3.0 m were verified to be distinct from those detected in the metagenomes at 3.0 m, according to two results: (1) Reciprocal BLAST analysis indicated that the flavobacteria from the viral samples (160,318 ORFs) and those from the microbial metagenomes (78,840 ORFs) shared only 2043 ORFs (0.01 and 0.03%, respectively); and (2) 16S rRNA genes predicted from both flavobacteria shared only 94% identity. Over 90% of the genome completeness of the flavobacteria in the viral samples was estimated using bacterial single-copy genes ( Table S3) and up to 99.8% completeness   using CheckM. Fragment recruitment analysis demonstrated that the flavobacteria had high similarity and coverage toward the genome of Nonlabens dokdonensis DSW-6 ( Figure S4), suggesting that the identified flavobacteria in the viral samples were closely related to Nonlabens. The draft genome size of Nonlabens sp. sh3vir. was ∼4.2 Mb, comprising 113 contigs (Figure 8), with complete gene suites for assimilatory sulfate reduction and thiosulfate oxidation as well as genes for denitrification (e.g., nirK and nosZ), nitrite reduction (e.g., nirBD), organic carbon oxidation, and carbon fixation.

DISCUSSION
We applied both 16S rRNA amplicon pyrosequencing and approaches to comprehensively explore microbial communities (bacteria, archaea, and viruses) and their potential metabolic capacity in Lake Shunet. Approximately 132 Gb of microbial and 133 Gb of viral metagenomic data were generated, which provided a much higher resolution of the microbial community structure than that in similar studies conducted in meromictic lakes Yau et al., 2013). These sequence data facilitated clustering the sequences into individual microbial bins, reconstructing a draft genome of Nonlabens sp. sh3vir. discovered in the viral samples, and generating near-complete viral genomes.

Archaeal and Bacterial Diversity Was High in Lake Shunet
Despite deep sequencing, rarefaction curves of the sequence reads plotted against OTU numbers did not reach a clear asymptote; this thus suggests that additional sequencing efforts would detect even greater diversity in this lake. Nevertheless, the observed OTU richness and diversity estimates are generally higher than those observed in other aquatic ecosystems, including two hypersaline meromictic lakes, namely the Ursu Lake (archaea H ′ = 3.34, bacteria H ′ = 5.45) and Fara Fund Lake (archaea H ′ = 4.06, bacteria H ′ = 3.74), Yellowstone Lake (Chao1 = 1074) Kan et al., 2011), freshwater reservoirs (bacteria H ′ = 5.34) (Tseng et al., 2013), and marine and mesohaline sites (archaea H ′ = 4.2) (Hugoni et al., 2015); this is probably due to the use of different primers or sequencing depths.

Variations of Archaeal and Bacterial Communities Were Significantly Explained by Physicochemical Parameters
Comparisons in archaeal and bacterial community compositions among Lake Shunet, Sakinaw Lake, Ursu Lake, and Fara Fund Lake suggested a unique microbial loop in Lake Shunet. However, we cannot exclude a potential error made from hypervariable regions of the 16S rRNA and the PCR bias caused by primer efficiency among studies. It might be because the environment surrounding Lake Shunet is different from the other three lakes. The observed archaeal and bacterial assemblages were dissimilar along the depths, and the variations were significantly explained by indigenous environmental parameters; these findings are in accordance with those of a study on Lake Kivu (Inceoglu et al., 2015). The sharp salinity gradient between mixolimnion and monimolimnion resulted in a lack of homogeneity of limnological parameters through stratification and likely resulted in various niches with distinct features that constituted some special microbial assemblages and metabolic capabilities (Lindstrom, 2001;Degermendzhy et al., 2010;Wear et al., 2013). Many bacterial genera, including Cyanobacteria, Pseudoalteromonas, Alcanivorax, and Rhodobacteraceae:NA, were increased by O 2 , temperature, pH, and NH + 4 in the mixolimnion; in this layer, key genes involved in oxygenic photosynthesis, nitrogen fixation, and nitrite reduction were detected in the bin sequences of these bacteria. At 5.5 m, H 2 S not only increased Desulfonatronum, Desulfobacteraceae, Desulfobacterium, and Halothiobacillus, which are either sulfatereducing or sulfur-oxidizing bacteria (Castro et al., 2000;Shi et al., 2011), but also enriched Candidatus Cloacamonas, which is probably a syntrophic bacterium involved in anaerobic digestion (Pelletier et al., 2008) and producing H 2 S during the breakdown of organic matter. In the archaeal community composition, Nitrososphaera and Methanospirillum were significantly explained by NO − 3 . Nitrososphaera is a typical ammonia-oxidizing archaeon that mediates nitrification (Spang et al., 2012), whereas Methanospirillum is a common syntrophic benzoate-degrading partner in benzoate degradation (Dolfing and Tiedje, 1988) and the genes involved in this degradation. The occurrence of benzoate degradation and the availability of nitrate as an electron acceptor possibly facilitated benzoate mineralization (Berry et al., 1987). This potential syntrophic relationship may have enriched both archaea at 3.0 m.

Genes Involved in Energy Metabolism Associated With Archaeal and Bacterial Assemblages
The potential KO profiles of energy metabolism, as well as archaeal and bacterial community compositions, differed significantly among the three sampling depths, because microbial functional profiles are usually associated with the community (Gilbert et al., 2010) or specific microbial groups (Debroas et al., 2009). For example, in the NMDS ordination, the relative abundance of Thiocapsa (a bacterial genus) and Methanococcoides (an archaeal genus) significantly explained the relative importance of potential KO involved in metabolism; specifically, Thiocapsa enriched K00395 (adenylylsulfate reductase, subunit B) involved in sulfur metabolism, and Methanococcoides enriched genes, including K00266 (glutamate synthase (NADPH/NADH) small chain), involved in nitrogen metabolism. Thiocapsa sp. dominated the bacterial community structure in Lake Shunet; they used hydrogen sulfide, thiosulfate, and elemental sulfur as electron donors in the anoxic chemocline zone (Caumette et al., 2004;Rogozin et al., 2012;Baatar et al., 2016). Reconstruction of elemental nutrient cycles in Lake Shunet also indicated an important role of Thiocapsa in sulfate reduction, consistent with a previous observation (Caumette et al., 2004). Allen et al. (2009) discovered that Methanococcoides burtonii isolated from Ace Lake could assimilate ammonia using a two-step glutamine synthetase (Mbur_1975) and glutamate synthase (Mbur_0092) pathway.

Unclassified Archaea and Bacteria Have an Important Role in Nutrient Cycles in Lake Shunet
Abundant untapped archaea and bacteria were detected at 5.0 and 5.5 m, respectively. This detection is similar to that in Sakinaw Lake, which had prevailing unassigned archaea in the chemocline. According to the elemental nutrient cycles of the lake that were reconstructed on the basis of binning results, several unknown or uncultured bacteria in high relative abundance (based on coverage of read counts) contributed to the elemental nutrient cycles at both 3.0 and 5.5 m. Furthermore, underexplored microbial groups were shown in various adaptive needs including bin 5.5-2, which harbored genes (betP/proP/opuA) involved in osmoadaptive mechanisms. The unknown archaea and bacteria were pre-dominantly diverse and played roles in nutrient cycling in this lake ecosystem.
To explore the potential roles of unclassified bacterial and archaeal OTUs (Table S4), we applied correspondence analysis to visualize the relationships between the 10 most abundant bacterial and archaeal OTUs and genes involved in energy metabolism. For example, gene K14083 (trimethylaminecorrinoid protein Co-methyltransferase, mttB) involved in methane metabolism was correlated with the abundance of two unknown archaeal OTUs (arcOTU_15 and arcOTU_72) and unknown bacterial OTUs (including bacOTU_14 and bacOTU_17; Figure S5). Notably, K00266 (glutamate synthase small chain) and K00262 (glutamate dehydrogenase) were potentially enriched by the specific unclassified archaea (arcOTU_31). This unknown archaeon might also assimilate ammonia using pathways similar to those in M. burtonii (Allen et al., 2009), suggesting that the unknown archaeon (different from all known archaeal phyla) might provide an ecological function for nitrogen assimilation in this ecosystem.

Reconstruction of Nutrient Cycles in Lake Shunet
The observed nutrient cycles differed along the three sampling depths. For example, aerobic carbon fixation was detected only in the oxic mixolimnion, and no genes involved in sulfite oxidation were found in the chemocline, whereas nitrite oxidation was detected only in the chemocline layer. These differences were attributed to variations in the microbial community composition in the three layers, because individual microbial groups evolved to address disparate core redox reactions (Falkowski et al., 2008).
On the basis of 16S rRNA genes and ORF read counts, Chroococcales was the most abundant bacterial group in the oxic mixolimnion and the only biotic source of O 2 production. Similar to Ace Lake, both the mixolimnion and monimolimnion of Lake Shunet were involved in CO production; this was due to the incomplete oxidation of organic compounds, indicating that CO oxidation might be an important pathway through which Rhodobacteraceae and unknown aerobic/anaerobic groups can generate energy, because numerous microbes can utilize CO for growth (Mörsdorf et al., 1992;Oelgeschläger and Rother, 2008).
Because both the upper and lower layers of the lake harbored diverse and abundant bacteria and archaea in the microbial loop and might compete for organic carbon, some bacteria might obtain energy from CO. The Wood-Ljungdahl pathway of acetyl-CoA biosynthesis in this lake is in the oxic zone, which conflicts with the findings of previous studies that have suggested that this pathway is commonly found in anaerobes (methanogens and acetogens). We suggest that suspended particles in the oxic layer could create macroscale oxic-anoxic interfaces for anaerobic metabolism.
Although most of the light is absorbed by dense purple sulfur bacteria in the chemocline, we detected genes involved in light-dependent metabolism, including those in monimolimnetic assemblage. We suppose that dominant photoautotrophic Thiocapsa and Cyanobacteria were detected in the monimolimnion because they attached to the sinking particle and some weak light penetrated into the bottom layer. This echoes the results of a previous study that measured bacteriochlorophyll a in the monimolimnion of Lake Shunet (Rogozin et al., 2009). Lake Shunet was reported to have methane concentrations of up to 32 µM (Kallistova et al., 2006); however, no genes involved in methane oxidation or methanogenesis were annotated in the current study, in contrast to previous research that has suggested that methane production commonly occurs in lake systems (Bastviken et al., 2004;Pasche et al., 2011). Our results also contrast with those reported for Ace Lake, where methane concentrations were up to 22.4 µM and genes for methane production were detected . We also observed prevalent methanogenic archaea in Lake Shunet; therefore, we speculate that methane was not an important carbon source in the anoxic zone or that conditions at the time of sampling were not favorable for methane production. Moreover, the lack of methane oxidation genes could be due to zero consumption from methanotrophic bacteria, which were not identified in the lake.
Many bins in high relative abundance potentially showed the capability of nitrate reduction and nitrogen fixation, leading to the production of ammonium, which is consistent with high ammonium concentrations detected in the mixolimnion in Lake Shunet. The N 2 source for nitrogen fixation was from the oxic zone. Complete processes of denitrification and nitrogen fixation co-occurred throughout the depths, supporting that nitrogen fixation might mainly provide bioavailable nitrogen to compensate for the loss of nitrogen by denitrification (Halm et al., 2009). In a previous study, abundant denitrifying bacteria were observed at equal levels (3.2-5.3 × 10 5 copy mL −1 ) from the surface to the bottom of the water in the Arabian Sea; however, denitrification activity, as probed using a stable isotope tracer method, was detectable only in the middle surface-the oxygen minimum zone (Ward et al., 2009). Therefore, we inferred that the presence of denitrifying bacteria did not always activate the genes involved in denitrification. Additionally, particles such as suspended sediment in the oxic layer could create macroscale oxic-anoxic interfaces for denitrification (Jia et al., 2016).
Regarding sulfur cycling, Thiocapsa (bins 5-6 and 5.5-18) were the dominant purple sulfur bacteria in the lake that also utilized thiosulfate (S 2 O 2− 3 ) instead of H 2 S as the electron donor and CO 2 as the carbon source under anoxic conditions (Caumette et al., 2004). Annotation of genome bins revealed that not only Thiocapsa but also other bacteria, including Hyphomonas neptunium-like bacteria (bin 3-11), Enterobacteriaceae (bin 5-2), and Desulfobacteraceae (bin 5.5-12), throughout the lake harbored genes involved in reducing sulfate to H 2 S as a terminal product (Figure 4). According to the relative abundance of bins, H 2 S seemed to be produced mainly in the chemocline. The H 2 S concentration was nearly zero in both the mixolimnion and chemocline but increased sharply to 400 mg L −1 in the monimolimnion. Some H 2 S possibly leached from the sediment and was generated by Thiocapsa and Desulfobacteraceae; in addition, because of gravity and increasing salinity, H 2 S accumulated at the bottom of the lake from the upper layers.
Three potential chemolithoautotrophic groups, namely Hyphomonas neptunium-like bacteria (bin 3-11), Enterobacteriaceae (bin 5-5), and Halomonas (bin 5.5-17), were identified and found to be potentially capable of deriving carbon from carbon fixation and energy from nitrate and sulfate reduction for growth. Specifically, only bin 5-5 harbored genes for the second step of nitrification (not detected in Ace Lake metagenomes). The three potential chemolithoautotrophs all belonged to nonphototrophic proteobacteria, and this finding is consistent with the discovery of chemolithoautotrophy among ubiquitous bacterial lineages (Swan et al., 2011). Nonlabens sp. Sh3vir.: Tiny Flavobacteria Discovered in a Viral Sample Notably, we detected numerous flavobacterial contigs at 3.0 m when analyzing viral metagenomic data. Flavobacteria are ubiquitous, strictly aerobic bacteria, comprising diverse members (Alonso et al., 2007). According to a comparison of 16S rRNA gene sequences, these flavobacteria were different from the other flavobacteria at 3.0 m, with low similarity (only ∼88%), and they had the highest phylogenetic similarity to Nonlabens dokdonensis (97% identity), of which cells were previously identified as Donghaeana dokdonensis, and the second highest similarity to Psychroflexus torquis ATCC 700755 (92% identity) (Yoon et al., 2006;Yi and Chun, 2012).
The reason for the detection of Nonlabens sp. sh3vir. in the viral samples warrants exploration. The cell morphology is generally rod-shaped or filamentous with a size range of 1.2-50 µm × 0.2-1.5 µm (Thomas-Jinu and Goodwin, 2004;Yoon et al., 2006;Yoshizawa et al., 2012); Psychroflexus are 0.2 µm in width (Bowman et al., 1998). However, the flavobacteria detected in the viral samples (<0.22 µm) in this study were Nonlabens sp. sh3vir., which has apparently not been described. Although all known flavobacteria are larger than 0.22 µm, we speculated that this Nonlabens sp. sh3vir. was likely smaller than 0.2 µm in wide and therefore could have passed through the 0.22-µm filter. Nonlabens dokdonensis is strictly aerobic, which corresponds to its exclusive detection in the oxic zone 3.0 m. In addition, some of its phenotypic characteristics, including the requirement of NaCl for growth and the capability of producing H 2 S and reducing nitrate (Yoon et al., 2006), were consistent with the environmental conditions of the surface water of Lake Shunet.
In the draft genome of the tiny flavobacteria, more than 70% of the genes were similar to the reference genome of Nonlabens dokdonensis DSW-6 and ∼3% of the genes were different, suggesting that the flavobacteria were yet to be identified Nonlabens sp. sh3vir. with a small size (<0.2 µm wide). Compared with microbial genomic bins, this genome is nearly complete; thus, we can provide more comprehensive insights into its ecological role in this lake. Moreover, this Nonlabens sp. sh3vir. is the first to be detected in an inland lake ecosystem, according to our review of the literature. The other Nonlabens sp. sh3virp. were found only in oceanic environments.
So far, Plasmaviridae have been reported only in Antarctic lakes (Rosario and Breitbart, 2011). Certain viral types are possibly enriched by local environmental conditions through selective pressure (Angly et al., 2006;Winter et al., 2013). With the exception of bacteriophages, terrestrial eukaryotic viruses were detected in Lake Shunet-specifically, abundant Nanoviridae at 5.5 m. Nanoviridae are DNA viral pathogens of plants (Lukert et al., 1995). Blast results indicated that such plants were mostly faba bean necrotic yellows, subterranean clover stunt, and milk vetch dwarf. Because Lake Shunet has an inflow from a stream , terrestrial viruses could flow into the lake during the snow melt period and subsequently accumulate at the bottom of the lake, probably due to their attachment to the sinking particles. Viral reference genomes are still few, and viral communities are only presented herein based on a minority of the annotatable data. However, through the use of a nonquantitative metavironmic method, the viral assemblages in Lake Shunet can be initially described and can serve as the reference baseline or further study of meromictic lakes in Asia.

Intimate Interactions of Viral, Bacterial, and Archaeal Assemblages
The viral richness and diversity of Lake Shunet significantly differed among the depths and were associated with bacterial and archaeal richness and diversity. Therefore, we concluded that bacterial and archaeal diversity was also influenced by local viral assemblages; this is because viruses have been deemed important vehicles, shuttling genetic materials and subsequently influencing speciation of microorganisms in aquatic ecosystems (Weinbauer and Rassoulzadegan, 2003). Additionally, our results indicate that some viral taxa were specific to certain bacterial bins, echoing the notion that viruses are strain-specific predators. Consequently, high archaeal and bacterial diversity could have resulted in the high viral diversity observed (Rohwer et al., 2009). Procrustes correlation analysis revealed a significant congruence among viral, archaeal, and bacterial communities in Lake Shunet, consistent with the interactions among archaeal, bacterial, and viral populations. Microbial bins at 5.0 and 5.5 m, but not at 3.0 m, seemed to be randomly infected by viruses from both depths, we predicted that viruses attached to particles and sank into the bottom layer.
A CRISPR is an antiviral defense system common in microbial genomes (Marraffini and Sontheimer, 2010). Infection traces viral sequences incorporated in the CRISPR arrayswere detected in the microbial metagenomes of Lake Shunet, suggesting previous interactions between viral and microbial communities (Berg Miller et al., 2012). For example, CRISPRs detected at both 5.0 and 5.5 m were randomly matched to viral sequences at both depths, whereas viral sequences at 3.0 m seemed to specifically match CRISPRs at 3.0 m only. The observed infection patterns were in accordance with the patterns characterized in interactions between bins and viruses; this finding supports the assertion that both the chemocline and monimolimnion might harbor similar viral assemblages, potentially infected diverse taxonomic groups, whereas 3.0 m depth contained viruses that were more strain-specific (Rohwer et al., 2009;Atanasova et al., 2012).

Lake Shunet Is a Unique Meromictic Lake
Lake Shunet has several characteristics and physicochemical profile patterns that resemble Ace Lake in Antarctica (Laybourn-Parry and Bell, 2014). To determine the specificity of Lake Shunet in terms of microbial communities and functions, we initially performed an NMDS analysis based on embedded COG profiles obtained from microbial metagenomes of both lakes ( Figure S6). In the NMDS ordination plot, Lake Shunet and Ace Lake were discriminated, indicating the existence of dissimilar microbial communities characterizing local biogeochemical processes. In addition, among meromictic lakes, Lake Shunet harbored unique prevalent bacterial and archaeal assemblages. Taken together, Lake Shunet is a unique meromictic lake in terms of microbial communities and potential functions; thus, it can serve as a baseline reference in studies of meromictic lake ecosystems, particularly in central Asia.