Bacterial Community Profiles of Cr-Contaminated Sites
There was considerable variability in physicochemical and nutritional properties of the contaminated soil, including gradients in Cr contamination (Fig.S1 and Table S1). Based on the Cr pollution level and nutrient conditions, these seven samples can be roughly divided into three groups: a slightly-contaminated group (L1, L2 and Z1 with available Cr concentrations of 0.11, 0.27, and 0.91 mg/kg, respectively), a moderately-contaminated group (L3 and Z2 with Cr concentrations of 6.76 and 7.09 mg/kg, respectively) and a highly-contaminated group (Z3 and Z4 with Cr concentrations of 413.84 and 465.42 mg/kg, respectively). Accordingly, there was a logical decrease in microbial biomass with the increase of Cr contamination and decrease in nutrient availability. The soil bacterial abundance, quantified by 16S rRNA abundance, fell dramatically by several orders of magnitude along the increasing Cr concentration gradient, from 9.78 to 5.97 log (copies)/gdw (per gram dry weight) in LZ, and from 8.89 to 5.96 log (copies)/gdw in ZY (Fig.1A).
Metagenomic analysis of dominant bacterial species and beta diversity revealed a significant difference in the bacterial communities of soils with different levels of Cr pollution. Phylogenetic analysis and heatmap of the 565 dominant species (with a relative abundance of > 0.1%) in LZ and ZY sites showed that the bacterial species from the comparable contamination levels (slight, moderate, and severe) were more closely related, while those from different contamination levels differed observably (Fig.1B), although they were all mainly from Proteobacteria and Actinobacteria. Consistently, redundant analysis (RDA) of bacterial communities showed that bacterial profiles in L1, L2 and Z1 with low available Cr concentration were more similar, and those in L3 and Z2 with moderate pollution levels were clustered together. The bacterial communities in Z3 and Z4 were dramatically different from those in other sites, possibly due to the severe degree of pollution and associated poor nutritional status (Fig.S2).
Virome Patterns along Cr contamination gradients.
In total, 5099 viral clusters (VCs) were identified from the 7257 viral contigs (>5 kb) originating from site LZ and 3021 VCs were identified from the 4561 viral contigs (>5 kb) originated from site ZY (Fig.S3A & B). As with most characterized viromes in various habitats [14, 32], the annotated phages in Cr-contaminated soil predominately belonged to Siphoviridae, Podoviridae and Myoviridae families (Fig.2A). The relative abundance of Siphoviridae, Podoviridae and Myoviridae families among the seven different soil samples was 33.3-93.4%, 0.7-31.9% and 0.5-12.2%, respectively (Fig.2A). Similar to the trend of bacterial abundance, the abundance of viral-like particles (VLP, determined by fluorescence microscopy imaging) decreased from 9.58 to 7.18 log(VLP)/gdw in site LZ and from 9.13 to 7.86 log(VLP)/gdw in site ZY with the increase of Cr contamination (Fig.2B).
Consistent with the trends observed in the bacterial communities, the phage populations among these seven different locations were significantly affected by soil Cr contamination and soil nutrient status. Cluster analysis of the top 20 viruses (with a relative abundance of > 1.0 %) showed that viruses from samples with similar physical-chemical properties were more closed related: L1 and L2 clustered together, Z1, Z2 and L3 clustered together, and Z3 was most closely clustered to Z4 (Fig.S4). RDA corroborated the viral composition patterns with virome of soil with the same level of Cr contamination clustered together (Fig.S5).
Compared with the prokaryotic viruses in the viral database (RefSeq, version 75) [33], the viromes of Cr-contaminated sites were highly novel and variant in terms of genetic profiles (Fig.3). Specifically, approximately 92.3% of viral contig sequences are not related to RefSeq virus sequences. Only 7.7% of the viral contigs had high genetic similarity to viruses recorded in RefSeq (weighted pairwise similarity scores of ≥ 1). Whereas 24.6% of the matched viral contigs in the contaminated soil were linked with viruses found in soil, the values decreased from 38.0% to 20.9% with the increase of pollution concentration in each sample (Fig.S6). Consistent with the similarities in viral profiles, the soil viral contigs with similar levels of Cr contamination were assigned with relevant habitats (Fig.S6). Interestingly, 14.6% of the matched viral contigs were linked with aquatic viruses. Besides, there were less than 10.0% viral contigs linked with viruses from sediment, human body and other habitats.
Variation in host-specific features
The host spectra of viruses were estimated via matching CRISPR spacers in the IMG/VR database to viral contigs (Fig.4). A total of 1044 hits were obtained between 646 phage contigs and 279 bacterial genera. Following the typical parasite-host pattern, the matched 279 host genera were mainly from eight bacterial phyla (Fig.4A & Fig.S7A), which were comprised of the most abundant phyla in the bacterial community (Fig.S7B). For example, the phage contigs linking to Proteobacteria and Actinobacteria, the most abundant phyla (36.4% and 33.0%) in the bacterial community, accounted for 42.4% and 34.3% of the total phyla, respectively. Whereas most soil viruses (731 contigs) seemed to be specific with contigs linked to one specific bacterial genus, a considerable fraction of viruses (363 contigs) may have broad host ranges with phage contigs linked to more than two bacterial genera. The relative abundance of the possibly polyvalent phages (having links to two or more distinct genera), which had higher fitness under high environmental stresses and low bacterial abundance [34, 35], was higher in soils with higher Cr levels (Kruskal-Wallis (W-K) test, p < 0.05). Specifically, the fractions of polyvalent phages in LZ site increased from 26.3% in L1 to 33.3% in L3, and their fractions in ZY site increased from 21.3% in Z1 to 55.2% in Z4 (Fig.4B).
In the low and medium contaminated sites, the relative abundance of Pseudomonas, Cronobacter, Klebsiella and Halomona genera within the bacterial community increased with the available Cr concentrations in soil (Fig.S8). These four bacterial genera are characterized by their biofilm formation capabilities, low membrane permeability and strong efflux pump systems conducive to bacterial survival under heavy metal stress [36, 37]. Accordingly, the relative abundance of phages linking to these four genera increased from 3.6% -13.9% in the low contamination sites (L1, L2 and Z1) to 27.3- 40.3% in the medium contamination sites (L3 and Z2) (Fig.4 & Fig.S8). However, the dominant bacterial genera linked by phage in Z3 and Z4 with high Cr levels (413.8 and 465.4 mg/kg) were significantly different from those in the previous slightly and moderately contaminated sites. The dominant genera of Z3 included Micropruina, Brevibacterium, Acidithiobacillus, and Zoogloea, and that of Z4 included Bacillus, Stenotrophomonas, Comamonas and Lysinibacillus. These bacteria are often detected in heavy metal-contaminated environments because of their strong resilience and adaptivity to high levels of salts and heavy metals [38-40]. Overall, the host-phage linkage corroborated phages as bacterial parasites depend heavily on the appropriate host flora to persist in Cr contaminated soils.
Variation in the relative abundance of lysogenic phage indicators (LPIs)
The annotation of viral genomes was performed by comparing the predicted proteins against the PFAM database (33.1) to explore the lysogenicity of virome. Viral contigs with integrase genes (Table S2) were counted as lysogenic phages [41, 42]. The relative abundance of lysogenic phages increased with Cr levels in LZ site: L1 (7.0%) < L2 (7.5%) < L3 (35.1 %) (Fig.5A). For ZY, the relative abundance of lysogenic phage formed a bell-shaped curve as the Cr level increased. The relative abundance of lysogenic phages increased from Z1 (8.3%) to Z2 (23.1%), while decreased in Z3 (14.7%) and Z4 (13.3%) (Fig.5A). There was a significant increase in the relative abundance of viral integrase genes in bacterial genomes in response to the increased available Cr levels (W-K test, p < 0.01): L1 (0.1%) ≈ L2 (0.1%) < L3 (0.3%) in LZ site, and Z1 (0.1%) < Z2 (0.4%) < Z3 (0.6%) < Z4 (0.7%) in ZY site (Fig.5B). Furthermore, more phages could be chemically induced from bacteria in more contaminated soil, which corroborated the host genome contained more prophage as Cr levels increased (W-K test, p < 0.05). In L1, L2 and Z1, the lysogenic fractions were all lower than 2 VLP/cell, while the lysogenic fractions increased to 5.45 and 4.14 VLP/cell in L3 and Z2, respectively. The corresponding fractions further increased to 19.14 and 22.84 VLP/cell in Z3 and Z4, respectively (Fig.5C).
Community-wide comparative genetic profiles of viromes
The annotation of viral genomes was performed by comparing the predicted proteins against the Non-supervised Orthologous Groups (eggNOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Nonmetric multidimensional scaling (NMDS) analysis of viral genomes revealed that the viral contigs from different locations were more similar in functional composition when the spatially distinct locations had similar Cr levels (Fig.S9). Specifically, the low Cr level samples (i.e., L1, L2, and Z1) clustered together, medium Cr level samples (L3 and Z2) clustered together, while high Cr level samples (i.e., Z3 and Z4) but distinct from the other samples.
Moreover, viruses carried multiple AMGs related to host metabolism, and their relative abundance increased significantly with available Cr concentrations in slight and medium contamination sites (Fig.S10). For example, the relative abundance of genes related to energy production and conversion (Group C) and lipid transport and metabolism (Group I) increased within the virome with available Cr concentrations increased from 0.11 to 6.97 mg/kg. However, at severe pollution locations (Z3 and Z4), the relative abundance of these genes in virome decreased significantly compared to Z2. Similar trends were observed for genes related to inorganic ion transport and metabolism (Group P), secondary metabolites biosynthesis, transport and catabolism (Group Q), intracellular trafficking, secretion, and vesicular transport (Group U) and defense mechanisms (Group V). These genes are well recognized for mediating host transport and ion secretion. Upregulation of these genes is generally beneficial to microbial energy conversion under environmental stress, and to enhance detoxification processes (e.g., efflux and chromate reduction) [43, 44].
Bacterial Cr detoxification mainly involves efflux of intracellular Cr (VI) ions through membrane transporters and enzymatic reduction of Cr (VI) into less toxic Cr (III) [45, 46]. KEGG databases were used to screen functional genes corresponding to microbial heavy metal detoxification (i.e., membrane transporter gene (Table S3) and reductase gene (Table S4) in viromes. As shown in Fig.6A, the relative abundance of reductase genes in viromes increased from 0.03% to 3.8% and membrane transporter increased from 0.1% to 16.5% as the concentrations of available Cr increased from 0.11 to 6.76 mg/kg, which demonstrated that phages can be important reservoirs of bacterial heavy metal resistance genes (MRGs) in contaminated sites. However, the relative abundance of bacterial MRGs in Z3 and Z4 decreased to 0.2% and < 0.1%, respectively. Notably, over 77.0% of the genes encoding membrane transporters and 61.7% of reductases were located on lysogenic phages (Fig.6B), by investigating whether these genes are carried by viral contigs with the integrase genes. This suggests that lysogenic phages may serve as carriers of MRGs in highly contaminated sites.
Case study of lysogenic phages with MRGs.
To corroborate that putative lysogenic phage-encoded AMGs (especially MRGs), four viral contigs (VC-2, VC-16, VC-20, VC-22) with integrase gene were randomly selected from the viromes as representatives and their genomic profiles were systematically analyzed. All VCs carried prokaryotic genes associated with heavy metal detoxification (Fig.6C). Notably, viral contig VC-2 encoded genes including Arsenate reductase, Gutathione S-transferase, Ferric iron ABC transporter, Fe-S oxidoreductase, Heavy metal RND efflux outer membrane protein, Cobalt/zinc/cadmium efflux RND transporter, and Molybdenum ABC transporter. Viral contig VC-20 encoded Dehydrogenase oxidoreductase protein and Permease of the drug/metabolite transporter (DMT) superfamily. For VC-22, there were genes encoding Ferredoxin-NADPH reductases, Sulfonate ABC transporter and Reductase.