Skip to main content

Enhanced mutualistic symbiosis between soil phages and bacteria with elevated chromium-induced environmental stress



Microbe–virus interactions have broad implications on the composition, function, and evolution of microbiomes. Elucidating the effects of environmental stresses on these interactions is critical to identify the ecological function of viral communities and understand microbiome environmental adaptation. Heavy metal-contaminated soils represent a relevant ecosystem to study the interplay between microbes, viruses, and environmental stressors.


Metagenomic analysis revealed that Cr pollution adversely altered the abundance, diversity, and composition of viral and bacterial communities. Host–phage linkage based on CRISPR indicated that, in soils with high Cr contamination, the abundance of phages associated with heavy metal-tolerant hosts increased, as did the relative abundance of phages with broad host ranges (identified as host–phage linkages across genera), which would facilitate transfection and broader distribution of heavy metal resistance genes in the bacterial community. Examining variations along the pollutant gradient, enhanced mutualistic phage–bacterium interactions were observed in the face of greater environmental stresses. Specifically, the fractions of lysogens in bacterial communities (identified by integrase genes within bacterial genomes and prophage induction assay by mitomycin-C) were positively correlated with Cr contamination levels. Furthermore, viral genomic analysis demonstrated that lysogenic phages under higher Cr-induced stresses carried more auxiliary metabolic genes regulating microbial heavy metal detoxification.


With the intensification of Cr-induced environmental stresses, the composition, replication strategy, and ecological function of the phage community all evolve alongside the bacterial community to adapt to extreme habitats. These result in a transformation of the phage–bacterium interaction from parasitism to mutualism in extreme environments and underscore the influential role of phages in bacterial adaptation to pollution-related stress and in related biogeochemical processes.

Video Abstract


Phages, bacterial viruses that constitute the most abundant and genetically most diverse biological entities in the biosphere, are critical components of microbial communities [1, 2]. The rapid development of viral metagenomics has advanced our understanding of phage abundance and composition in various ecosystems [3,4,5,6,7], allowing us to identify their functional importance more comprehensively. The sheer abundance and diversity of phages, together with their multifaceted impacts on microbial hosts, strongly indicate that phages can influence entire macroscopic ecosystems by affecting the native microbial communities [8,9,10,11,12]. However, previous virome studies mainly focused on phage species and genetic profiles from different ecosystems, while phage–bacterium interactions have been heretofore relatively understudied [5, 13]. To better understand the ecological and evolutionary roles of phages in ecosystems, it is important to elucidate phage–bacterium interactions and their response to environmental stresses.

Phages are obligate intracellular bacterial parasites who rely on the presence of appropriate hosts to persist and propagate [14, 15]. The classical view of phage–bacteria parasitism holds that phages propagate at the expense of their bacterial hosts and redirect host resources to produce virions [16]. Beside this parasitic pattern often observed under favorable conditions, mutualistic phage–bacteria interactions also occur under unfavorable circumstances (e.g., solar radiation or nutrient limitation) [17, 18]. Specifically, lysogenic phages can promote horizontal gene transfer by lysogenic conversion or transferring packaged bacterial genes to new hosts, which may expand the host metabolic profile and enhance microbial environmental adaptability [11, 19,20,21,22]. In return, the phages that are integrated into the host genomes in dormant forms (i.e., prophages) can take advantage of microbial environmental resistance to avoid direct exposure to harsh environments [23]. Typically, the extent of phage–bacteria mutualism is positively correlated with the abundance of lysogenic phages in a virome, lysogenized bacteria in a microbiome, and auxiliary metabolic genes (AMGs) carried by lysogenic phages [19].

Heavy metal contamination of soil poses a serious threat to public health and ecological security due to its widespread occurrence, persistence, and ecotoxicity [24, 25]. Heavy metal-contaminated soils generally exhibit a decreasing geochemical gradient from the contamination source [26, 27]. Such soils are recognized as ecosystems under complex stresses and are both relevant and convenient to study the interplay between phages, bacteria, and environmental stressors. However, the viromes in heavy metal-contaminated soils remain poorly characterized, and the response of phage-bacterium interactions to heavy metal-induced environmental stress is largely unknown. This critical knowledge gap should be addressed to inform microbiome adaptation to hostile environments and the role of phages in related biogeochemical processes. Advancing understanding of phage–bacterium interactions under different heavy metal stress may also help better assess the ecotoxicity of heavy metals and potential of microbiome engineering for contaminated soil remediation.

To this end, the objective of this study was to explore the variability of soil viromes and phage–bacterium interactions along gradients of environmental stress. Chromium (Cr)-contaminated sites were chosen to represent heavy metal soil contamination, due to the high frequency of occurrence and well-recognized ecotoxicity [28, 29]. Soil samples were taken along Cr gradients in one lightly polluted site (Luzhou, LZ for short) and one heavily polluted site (Zhangye, ZY for short). Metagenomic analyses for total environmental DNA and viral DNA were separately conducted to inform indigenous bacterial and viral profiles. Bacteria–phage linkages, commonality of lysogenicity within the virome, abundance of phage-carried AMGs and the lysogenized fraction of bacteria were investigated to demonstrate how variations in the virome and its functions were related to the indigenous bacteria under Cr-generated environmental stress. This metagenomic study reveals an enhanced mutualistic symbiosis between phages and bacteria with elevated Cr-induced environmental stress and highlights the roles of phages in mediating microbiome adaptation to hostile environments.


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 6.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.75 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).

Fig. 1
figure 1

Profile of bacterial communities in Cr-contaminated soils. A The biomass of bacteria in Cr-contaminated soil. The petal lengths of the rose represent the value of log (copies/g soil) quantified by fluorogenic quantitative PCR reactions with 16S rRNA gene, which generally decreased with the increase of Cr contamination in both sites. B Redundant analysis (RDA) of bacteria in different soil samples. The results showed nitrate, ammonium, total nitrogen, organic matter and total phosphorus had the closest correlation with the composition of bacteria. The bacterial communities of each sample clustered most closely with samples under similar geochemical conditions. Specifically, L1, L2 and Z1 with low available Cr concentration were more similar, and L3 and Z2 with slightly higher pollution levels were clustered together. Z3 and Z4 were dramatically different from those in other sites

Metagenomic analysis of species richness, dominant bacterial species and beta diversity revealed a significant difference in the bacterial communities of soils with different levels of Cr pollution. Specifically, the Chao1 and ACE index showed that the increase in soil Cr concentration adversely affected the species richness of microflora (Fig. S2). Moreover, 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. S3), 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. 1B).

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.S4A & B). As with most characterized viromes in various habitats [13, 30], 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).

Fig. 2
figure 2

Profile of viral communities in Cr-contaminated soils. A The composition of the viral community at the family level from seven different locations at two Cr-contaminated sites shows that viral populations between these 7 different locations exhibited distinct profiles. B The abundance of viral particles in Cr-contaminated soil counted by fluorescence microscope. The petal lengths of the rose represent the abundance of viral particles per gram of soil (log (VLP/g soil)), which decreased with the increase of Cr contamination in both sites

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 closely related: L1, L2 and Z1clustered together; Z1, Z2, and L3 clustered together; and Z3 was most closely clustered to Z4 (Fig. S5). RDA corroborated the viral composition patterns with virome of soil with the same level of Cr contamination clustered together (Fig. S6).

Compared with the prokaryotic viruses in the viral database (RefSeq, version 75) [31], 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. S7). Consistent with the similarities in viral profiles, the soil viral contigs with similar levels of Cr contamination were assigned with relevant habitats (Fig. S7). Interestingly, 14.6% of the matched viral contigs were linked with aquatic viruses, and there were less than 10.0% viral contigs linked with viruses from sediments, the human body, and other habitats.

Fig. 3
figure 3

Relating viral contigs from Cr-contaminated soils to known viral sequence space. Viral contigs recovered from Cr-contaminated soils were clustered with all RefSeq (v 75) viral genomes or genome fragments with genetic connectivity to these data. Each node colors indicate sample and reference viral sequences’ habitat of origin. The edges (lines) between the nodes represent statistically weighted pairwise similarity scores of ≥ 1. The results show that viromes at different sites along a Cr-contaminated soil gradient have distinct specificity. And only a few of recovered viral contigs were associated with viruses recorded in NCBI databases (RefSeq, v 75) from soil, water, human, and sediment 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.S8A), which were comprised of the most abundant phyla in the bacterial community (Fig. S8B). 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 (409 contigs) seemed to be specific with contigs linked to one specific bacterial genus, a considerable fraction of viruses (237 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 [32, 33], was higher in soils with higher Cr levels (Kruskal-Wallis (W-K) test, p < 0.05). Specifically, the fractions of polyvalent phages in the LZ site increased from 26.3% in L1 to 33.3% in L3, and their fractions in the ZY site increased from 21.3% in Z1 to 55.2% in Z4 (Fig. 4B). In Cr-contaminated soil, the top ten bacterial genera infected by polyvalent phages were Pseudomonas, Salmonella, Cronobacter, Enterobacter, Escherichia, Klebsiella, and Shigella from Proteobacteria phylum, and Salinispora, Micromonospora, and Actinomyces from Actinobacteria phylum.

Fig. 4
figure 4

Predicted virus–host linkages in Cr-contaminated soils. A A total of 1044 hits were obtained between 646 phage contigs and 279 bacterial genera, and these host genera were mainly from 8 bacterial phyla. From left to right and the first columns show the distribution of bacterial hosts at the phylum level; the second column depicts viral contigs linking to hosts in the seven soil samples; the third column shows the distribution of the bacterial host assigned with CRISPR approaches at the genus level. The height of the column section represents the abundance of viruses and hosts. The results revealed that these viral host species predicted in the soil polluted by high Cr concentration mostly belong to the bacterial genera which have been reported to have strong adaptivity to better survive in heavy metal-polluted habitats. B The relative abundance of polyvalent phages within the viral communities of Cr-contaminated soils. Polyvalence was found to be consistently higher in soil samples with higher available Cr concentrations

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. S9). 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 [34, 35]. 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.S9). 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 [36,37,38]. Overall, the host–phage linkage corroborated phages as bacterial parasites depending 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 [39, 40]. The relative abundance of lysogenic phages increased with Cr levels in the LZ site: L1 (7.0%) < L2 (7.5%) < L3 (35.1 %) (Fig. 5A). For ZY, the relative abundance of the 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.15%), 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 the LZ site, and Z1 (0.1%) < Z2 (0.4%) < Z3 (0.6%) < Z4 (0.7%) in the ZY site (Fig. 5B). Furthermore, more phages could be chemically induced from bacteria in more contaminated soil, which corroborated the host genome that 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.140 and 22.84 VLP/cell in Z3 and Z4, respectively (Fig. 5C).

Fig. 5
figure 5

Lysogen occurrence within the soil virome along Cr concentration gradients. A Co-occurrence networks of viral contigs in seven soil samples. Nodes (circles) represent viral contigs. Orange represent the contigs containing lysogenic phage indicators; dark blue nodes are all other viral contigs. Shared edges (lines) between nodes indicate a correlation, and the more edges in a node indicate a closer genetic relationship between contigs. The values to the right of each sample name are the relative abundances of lysogenic phages within the corresponding sample. As the distance between sampling location and source of contamination grew, the relative abundance of lysogenic phages decreased in the lightly polluted LZ site. In the heavily polluted ZY site, an increasing trend from sample Z1 to Z2 was observed, which then showed a decreasing trend between Z2 and Z4. The relative abundance of phage integrase genes found in bacteria (B) and inducible lysogenic fractions (C) at the sampling locations consistently increased in tandem with Cr concentrations

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. S10). 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 slightly and moderately contaminated sites (Fig. S11). 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 with Z2. Similar trends were observed for genes related to inorganic ion transport and metabolism (Group P); secondary metabolite 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) [41, 42].

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) [43, 44]. 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.

Fig. 6
figure 6

Genomic analysis of viral genomes in Cr-contaminated soils. A The relative abundance of viral genes coding membrane transporter (dark blue) and reductase (gray) were annotated by KEGG datasets (Table S3 and Table S4). There was a significant increase in the relative abundance of these genes associated with reduction and transport of heavy metals in L3 and Z2. B The distribution ratio of genes coding reductase (left) and membrane transporter (right) on lysogenic phages (blue) and other viruses (green), which were are calculated by determining whether these genes are carried by viral contigs with the integrase genes. It shows that over 60% of them were on lysogenic phages, so lysogenic phages are an important reservoir for MRGs. C Genome map of four viral contigs containing both the phage integrase gene and the metal detoxification gene. Key of gene families and their colors are found in the figure, and arrow length within the gene maps correspond to the size of the open reading frame (ORF). Detailed function descriptions of the seven viral contigs are listed in Table S5

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, Glutathione 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.


Phages, being the most common biological entity on earth, are a great untapped resource with possible applications in medical, industrial, and agricultural fields [45]. To identify the ecological and evolutionary role of phages in the biosphere and reveal the tactics of microbial adaptivity to hostile environments [46, 47], we systematically investigated variations of the viral community profile and phage–bacterium interactions along horizontal gradients of Cr-induced environmental stresses at two contaminated sites. The abundance and composition of viral communities varied greatly alongside their host communities across contamination and geochemical gradients. As Cr-induced environmental stresses increased, phage communities exhibited an enhanced beneficial relationship with the bacterial communities that is reflected by elevated relative abundance of lysogenic phages in the viromes and higher lysogenic fractions of bacteria. Furthermore, viral genomic analysis showed that phages under higher Cr-induced stresses carried more AMGs that can contribute to microbial heavy metal detoxification and survival in hostile environments. These findings inform variations in interactions between phages and microbes under adverse conditions and advance our understanding of the contribution of phages in biogeochemical cycles.

As is common for heavy metal contamination sites, increased levels of Cr were associated with decreased nutrient availability at both sites LZ and ZY (Table S1). Accordingly, the toxicity of Cr and the accompanying change in trophic conditions significantly altered both soil microbial and viral communities (Figs. 12). Notably, viral and bacterial communities from soil samples with similar Cr contamination levels rather than spatial proximity had more similarities in composition and abundance of functional genes (Fig. 1B, Fig. S6, & Fig. S10). Interestingly, the viromes in Cr-contaminated soils were mostly not linked to known viruses recorded in the RefSeq database (Fig. 3), possibly due to the unique physical and chemical properties of Cr-contaminated sites compared with other characterized habitats (uncontaminated soil, surface water, human gut, and sediment) [48]. This suggests that viromes of metal-contaminated soils (and possibly other extreme environments) might represent an untapped gene reservoir of biological agents with unknown potential for medical and environmental applications.

Phages as obligate intracellular parasites rely on the presence of appropriate hosts to persist and propagate [14, 15], which is corroborated by a high degree of consistency between the predicted phage host and the dominant bacteria in soils (Fig. 4A, Fig.S8 & Fig.S9). Relative abundance of bacterial genera (e.g., Micropruina, Brevibacterium, and Bacillus) with strong heavy metal resistance and survival capability under adverse conditions [49] increased in heavily Cr-polluted environments; CRISPR-based host linkage data revealed a tandem increase in the relative abundance of phages dependent on these resistant hosts. From an evolutionary perspective, bacteria with stronger resistance will have higher biomass and more stable intracellular environments under adverse conditions, providing better survival and reproduction conditions for their corresponding phages [50]. Furthermore, broad host range phages which can use multiple hosts to propagate are enriched in nutrient-poor environments or in conditions of low host density [33]; this research corroborates this trend for the first time within a Cr-contaminated area showing that the polyvalent strategy is viable to enhance phage fitness under severe environmental stress (Fig. 4B).

Furthermore, the capacity and tendency of phages towards a lysogenic lifestyle were generally enriched when faced with elevated levels of heavy metal pollution (Fig. 5). This corroborated the theory that lysogeny is a highly refined relationship whose development may be closely related to growth in unfavorable or even hostile conditions [19, 30, 51]. Both computational and experimental analyses demonstrated that the closer to the source of contamination, the higher abundance of prophages within the genomes of bacterial members of the community (Fig. 5B & C). When faced with hostile environments, phages capable of entering lysogeny have shown a higher inclination to integrate into their host genome as a prophage, establishing a more mutualistic relationship [16].

The relative abundance of lysogenic phages within the free virome in highly contaminated samples (Z2, Z3, Z4, and L3) increased significantly above that of low contamination samples (Z1, L1, and L2) (Fig. 5A). Interestingly, when examining the ratio of lysogenic to lytic phages in the free virome, the extremely contaminated samples (Z3 and Z4), had lower relative abundances of lysogenic phages, decreasing to nearly one half of that for moderately contaminated Z2 at the same sampling site. To explain this drop, we posit that with further deterioration of the surrounding environment, lysogenic phages may transform into prophages and lurk in the host body [52], which was corroborated by the increase of prophage in bacterial genome (Fig. 5B & C). Therefore, the number of free but lysogenic capable phages in the ultra-harsh environment would decrease. These findings inform the potential preference of lysogenic viruses in more extreme environments, suggesting that the lysogeny should be favorable for phages under heavy metal pollution stress. However, it should be noted that the host spectra estimated via the CRISPR-based approach may need further verification by other host-prediction models [53] or combined experimental and computational approaches [54].

Previous studies reported that phage-carried AMGs can contribute to the genetic and phenotypic evolution of their hosts and ultimately provide the hosts with diversified competitiveness and environmental adaptability [16, 52]. In our study, many MRGs encoding machinery for the detoxification of heavy metals from bacteria were detected in the virome (Fig. 6A & C). Our analysis also highlights a consistent trend between the relative abundance of MRGs and lysogenic phages within the virome (Figs. 5A & 6A). That is mainly because most MRGs are distributed on lysogenic phages (Fig. 6B). Therefore, when lysogenic phages integrate into host genomes in prophage form, the relative abundance of MRGs in the free virus particles also decreases. Lysogenic phages carrying AMGs (e.g., MRGs) suggest the potential of lysogenic phages to transfer beneficial genes in the microbial community and thus improve the survival ability of hosts in extreme environments [55]. Therefore, there is strong tendency toward mutualism relationship between phage and their hosts based on lysogenic integration with the further aggravation of Cr pollution in soil. In this relationship, lysogenized bacteria provide a more stable environment for phages [20], in return for which they may benefit from the expression of phage-carried AMGs [56,57,58]. Furthermore, the discovery of MRGs in the virome highlights the influence of viral communities on biogeochemical process in contaminated sites, which may inspire phage-based strategies (e.g., virome transplant) for microbiome engineering [60]. Given that phages are recognized as significant vectors of antibiotic resistance genes and pathogenic genes in various environments [11], the importance of bacteria–phage interactions in mediating bacterial resistance and possible pathogenesis may be underestimated, and their influence on public health warrants further investigation.


This study of virome and metagenome interactions in soils subjected to various levels of stress associated with Cr contamination revealed that phage–bacterium interactions change from antagonistic (parasitic) at low stress towards a more mutualistic relationship with increasing Cr concentrations. This mutualistic symbiosis is mainly based on the lysogenic integration of lysogenic phages. That is, the lysogenic phages carrying MGRs can improve their host’s ability to resist toxicity of heavy metals by integrating themselves into the host genome, which in turn facilitates expression of resistance genes and reproduction of related phages. This demonstrates that viral populations and their functions have a profound influence on the adaptation of prokaryotic communities and their ability to withstand adverse geochemical conditions. In addition, efforts to remediate and restore heavy metal-contaminated sites could benefit from advanced understanding of interactions between bacterial and phages and their manipulations to enhance redox processes that reduce metal toxicity and mobility.


Site description and sampling

Seven soil samples were collected from two Cr residue-contaminated sites, located in Zhangye City (ZY, N 38° 25′ 49″ E 100° 48′ 41″) with relatively high chromium pollution and Luzhou City (LZ, N 28° 54′ 24″, E 105° 31′ 26″) with light chromium pollution, in September 2019. Soil samples are collected at intervals of 50 m, from the Cr-slag dump center towards the edge of the site. According to the site area, 3 soil samples (L1, L2, and L3) were collected in LZ, and 4 soil samples (Z1, Z2, Z3, and Z4) were collected in ZY. Each soil sample was a manually homogenized composite of five subsamples collected at the surface layer (5–10 cm) and sieved through a size 20 mesh. The details of sampling locations are shown in Fig.S1. The physical and chemical properties (i.e., pH, organic matter, total phosphorus, total sulfur, total potassium, available phosphorus, nitrogen, nitrate nitrogen) of soil samples were determined by standard methods [61]. The total Cr, available Cr (i.e., the total amount of water-soluble and exchangeable Cr), and hexavalent Cr concentrations were measured by ICP-MS (Agilent 7700x) according to the USEPA Method 7196A [62]. The concentrations of available Cr in the seven samples varied from 0.11 to 465 mg/kg (Table S1).

Quantification of bacteria and virus in soil

Fluorogenic quantitative PCR reactions was conducted with 16S rRNA gene V4–V5 primer pairs to determine the biomass of bacteria in soil samples from ZY and LZ (Primer sequences (5’-3’) were 16S-515F: “GTGCCAGCMGCCGCGG” and 16S-907R: “CCGTCAATTCMTTTRAGTTT”). To determine the biomass of viruses in soil samples from ZY and LZ, the inverted fluorescence microscope (Nikon, eclipse Ti-S) was used to observe soil virus particles stained with SYBR Gold fluorescent dyes as previously described [63], and Fig. S12 is a typical image, in which viruses are the most numerous tiny dots, while bacteria are larger with definite shapes and edges.

Total DNA extraction, metagenomic sequencing, and analysis

FastDNA Spin Kit for Soil (MP Bio, Irvine, CA) was used to extract total DNA from the soil samples according to the manufacturer’s instructions. DNA quantification was performed with a Qubit 3.0 fluorometer using the dsDNA system (Invitrogen, Waltham, MA). According to the Illumina TruSeq DNA Sample Preparation Guide, by using NEBNext Ultra II DNA Prep Kit (New England Biolabs, Ipswich, MA) and NovaSeq 6000 S4 Reagent Kit (300 cycles) (Illumina, San Diego, CA), genomic libraries were constructed and sequenced on Illumina Hiseq 4000 platform (2 × 150 bp, paired-end reads) at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). Finally, 16-GB sequence data was obtained for each of the samples. The quality control of raw reads was conducted by Cutadapt (v1.17) and FASTP (0.20.0) [64], and high-quality reads were used to assemble sequences with MEGAHIT v1.0. Minimap2 [65], Samtools [66], and blast2lca [67] were used to complete the construction of non-redundant sequence sets. After predicting the open reading frame (ORF) in sequences using MetaGeneAnnotator [68], redundant sequences were removed based on 90% sequence similarity by using CD-HIT [69]. Soap.coverage ( was used to calculate contig abundance and gene TPM (transcripts per million) value using a script ( The functional annotation of predicted ORFs was performed through assigning predicted ORF sequences to the PFAM 33.1 [70], KEGG database [71], and eggNOG v5.0.0 with a cutoff of E value < 0.001 [72]. Taxonomic classification at the species level was annotated via the NCBI-NT database (v2016-6-19) using Kaiju with the default run mode Maximum Exact Matches [73].

Virus particle elution and verification

To separate viral nucleic acid from bacterial nucleic acid, several measures were taken to extract and enrich particles from soil following previous protocols (Fig. S13) [74]. Briefly, fresh soil sample was suspended in 4 °C potassium citrate buffer (10 g/L C6H5K3O7, 1.92 g/L Na2HPO4·12H2O, and 0.24 g/L KH2PO4, pH = 7). After bath sonication in an ice–water mixture for 3 min at 47 kHz, the suspension was centrifuged (7000 rpm, 4 °C, 10 min) to obtain the virus-containing supernatant. This was subsequently followed by centrifugation and filtration (0.2 μm) to remove impurities larger than 0.2 μm [73]. It should be noted that this filtration process may exclude viruses larger than 0.2 μm, leading them not to be included in further analysis of the viral population. The virus particles in the filtrate were initially concentrated by tangential flow filter cartridge (TFF, Sartorius Vivaflow 50 30000MWCO PES) followed by precipitation by polyethylene glycol 8000. The obtained virus pellets were treated with DNase and RNase to remove free DNA and RNA before viral DNA extraction [75]. In addition, virus particles were verified by transmission electron microscope (Zeiss LSM710) with the phosphotungstic acid counterstaining method (Fig. S14) [76]. It should be mentioned that the enrichment method of virome in this study is not perfect in the acquisition of giant phages.

Viral DNA extraction and metagenomic sequencing

Viral DNA was extracted using the AllPrep PowerViral DNA/RNA Kit (MP Bio) following the manufacturer’s instructions. The extracted DNA was subjected to whole genome amplification (KAPA HiFi HotStart ReadyMix) to meet the metagenomic sequencing requirements. Libraries were constructed before sequencing on Illumina Hiseq4000 platform as described in “Total DNA extraction, metagenomic sequencing and analysis”. The protocol of quality control and assembly of virome is also the same as above. It should be mentioned that the enrichment method of viral nucleic acid in this study is not perfect in the acquisition of RNA viral information.

Viral contig identification, taxonomic classification, and functional gene annotation

To eliminate potential bacterial contamination, “Virsorter + vHMM” method was used for the identification and annotation of viral contig > 5 kb in length [14, 77]. According to previous reports, VirSorter can provide near-perfect identification (> 95% recall and 100% precision) [78], and vHMM demonstrates a specificity of 99.6% for viral detection with a 37.5% recall rate [14]. Specifics are as follows: firstly, Virsorter was used with default parameter values, and the categories 1, 2, 4, and 5 of Virsorter predictions were considered as viral contigs. Secondly, the protein sequences of the remaining contigs were scanned with vHMM pipeline under 3 criteria (refer to (a) the contigs had at least 5 hits to viral protein families, and the total number of genes covered with KO terms on the contig was < 20%; the total number of genes covered with Pfams ≤ 40%; and the number of genes covered with viral protein families > 10%; (b) the contigs were selected as viral contigs when the number of viral protein families on the contig were equal or higher than the number of Pfams (bacterial genes); (c) the contigs for which the number of viral protein families was equal or higher than 60% of the total genes were classified as viral contigs. Viral contigs must satisfy at least one of the above three criteria. Then, vcontact2 was used to annotate the taxa of these viral contigs [19]. The functional annotation of non-redundant proteins of viral contigs was performed through assignment of predicted proteins to the PFAM, KEGG, and EggNOG databases with a cutoff of E value < 0.001 [72].

CRISPR-based phage–bacteria linkage

Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) were used to predict phage–bacteria linkage based on the sequence similarity between the spacers in microbial CRISPR regions and the protospacers in viral genomes [19]. The sequences of all viral contigs were aligned to the spacer sequences deposited in the IMG/VR database ( Specifically, all viral contig and spacer sequences from the shared genomic regions were queried by BLASTn (E value 10−10 and 100% nucleotide identity) to obtain the information about the potential host of the virus and finally link the virus to its bacterial host.

Prophage induction assay by mitomycin-C

Induction of prophages was performed as previously described [51]. Briefly, 1 g of each of the original soil samples were suspended in 9 ml sterilized deionized water, and shaken for 20 min at the maximum speed, 20,000 rpm. Then, these resuspensions were inoculated into sterile water-soluble organic carbon (WSOC) solution both with and without mitomycin-C (final concentration of 1 μg/ml). These cultures were incubated at room temperature, shaking at 200 rpm for 24 h under in the dark. To harvest the induced and free phages, the cultures were passed through a 0.2-μm sterile membrane. The filtrates were then treated with glutaraldehyde (final concentration was 2.5%) as a fixative and stored at 4 °C for 24 h prior to staining.

Staining and fluorescence counting were performed as previously described by Anand, et al. [63]. After counting, the number of phages induced from a single bacterium was calculated using the following formula:

$$ \mathrm{X}=\left({\mathrm{V}}_{\mathrm{i}}-{\mathrm{V}}_{\mathrm{ck}}\right)/B $$

Here, Vi is the abundance of viral-like particles in the treatment with mitomycin-C, Vc is the abundance of viral-like particles in control treatments, and B is the bacterial abundance in the corresponding soil sample (by fluorogenic quantitative PCR reactions).

Data analysis

The data analysis and visualization process involved in this study were conducted by Excel 2013 and R ×64 3.6.2. In order to clarify the correlation of viral contigs recovered in Cr-contaminated soil and other known viruses, the gene similarity score between the acquired virus contigs and known bacterial and archaeal viruses (NCBI RefSeq, version 75) was calculated according to Trubl et al. [31] (Genome pairs with a similarity score of 1 were considered remarkably similar). Then Cytoscape3.7.1 ( was used to visualize the co-occurrence network, using an edge-weighted spring-embedded model. The Maximum-likelihood phylogenetic tree and heat maps of dominant bacteria (relative abundance greater than 0.1%) were produced using MEGA7 with the JTT model and iTOL ( The RDAs of viruses and bacteria were performed using R packages ggbiplot to analyze the viral and bacterial composition and diversity of different samples. In addition, functional genes annotated in KEGG databases corresponding to microbial heavy metal detoxification (i.e., membrane transporter gene (Table S3) and reductase gene (Table S4)) in the metagenome and virome were isolated to determine the detoxification potential of the microbial community. The differences between genetic profiles of viral communities were assessed by NMDS analysis utilizing R package vegan [79]. Based on PFAM databases, phage integrase family genes (as listed in Table S2) in the metagenome and virome were screened out to investigate changes in the relative abundance of lysogenic phages [29]. The statistical significance of differences at different Cr-contamination levels (i.e., the lysogenicity fraction of bacterial community, the relative abundance of MRGs, lysogenic phages, and polyvalent phages among soils with different pollution levels) were determined at the 95% confidence level (p < 0.05) using the non-parametric Kruskal-Wallis test (K-W test) [30].

Availability of data and materials

All raw sequence data generated in this research has been deposited in the Genome Sequence Archive (GSA) in BIG Data Center, Beijing Institute of Genomics, Chinese Academy of Sciences, under accession numbers CRA003088 and CRA003796, and are publicly accessible at





Auxiliary metabolic gene


Heavy metal resistance gene


Viral cluster


National Center for Biotechnology Information


Evolutionary genealogy of genes: Non-supervised Orthologous Groups


Kyoto Encyclopedia of Genes and Genomes


Redundant analysis


Nonmetric multidimensional scaling


Polymerase chain reaction


Water-soluble organic carbon


  1. Martha RJC, Andrew DM, Andrey VL, Shaun H. Phages in nature. Bacteriophage. 2011;1(1):31–45.

    Article  Google Scholar 

  2. Wang YL, Jiang XT, Liu L, Li B, Zhang T. High-resolution temporal and spatial patterns of virome in wastewater treatment systems. Environ Sci Technol. 2018;52(18):10337–46.

    Article  CAS  PubMed  Google Scholar 

  3. Petrovich ML, Zilberman A, Kaplan A, Eliraz GR, Wang YB, Langenfeld K, et al. Microbial and viral communities and their antibiotic resistance genes throughout a hospital wastewater treatment system. Front Microbiol. 2020;11:13.

    Article  Google Scholar 

  4. Cao B, Xu H, Mao CB. Phage as a template to grow bone mineral nanocrystals. Methods Mol Biol. 2014;1108:123–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Suttle CA. Viruses in the sea. Nature. 2005;437(7057):356–61.

    Article  CAS  PubMed  Google Scholar 

  6. Motlagh AM, Bhattacharjee AS, Coutinho FH, Dutilh BE, Casjens SR, Goel RK. Insights of phage-host interaction in hypersaline ecosystem through metagenomics analyses. Front Microbiol. 2017;8:352.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Manrique P, Dills M, Young MJ. The human gut phage community and its implications for health and disease. Viruses-Basel. 2017;9(6):141.

    Article  CAS  Google Scholar 

  8. Pan D, Watson R, Wang D, Tan ZH, Snow DD, Weber K. Correlation between viral production and carbon mineralization under nitrate-reducing conditions in aquifer sediment. Isme J. 2014;8(8):1691–703.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Nair RR, Vasse M, Wielgoss S, Sun L, Yu YTN, Velicer GJ. Bacterial predator-prey coevolution accelerates genome evolution and selects on virulence-associated prey defences. Nat Commun. 2019;10(1):4301.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Kotay SM, Datta T, Choi JD, Goel R. Biocontrol of biomass bulking caused by Haliscomenobacter hydrossis using a newly isolated lytic bacteriophage. Water Res. 2011;45(2):694–704.

    Article  CAS  PubMed  Google Scholar 

  11. Calero-Caceres W, Ye M, Balcazar JL. Bacteriophages as environmental reservoirs of antibiotic resistance. Trends Microbiol. 2019;27(7):570–7.

    Article  CAS  PubMed  Google Scholar 

  12. Bhattacharjee AS, Choi J, Motlagh AM, Mukherji ST, Goel R. Bacteriophage therapy for membrane biofouling in membrane bioreactors and antibiotic-resistant bacterial biofilms. Biotechnol Bioeng. 2015;112(8):1644–54.

    Article  CAS  PubMed  Google Scholar 

  13. Shkoporov AN, Hill C. Bacteriophages of the human gut: the “known unknown” of the microbiome. Cell Host Microbe. 2019;25(2):195–209.

    Article  CAS  PubMed  Google Scholar 

  14. Paez-Espino D, Eloe-Fadrosh EA, Pavlopoulos GA, Thomas AD, Huntemann M, Mikhailova N, et al. Uncovering Earth’s virome. Nature. 2016;536(7617):425–30.

    Article  CAS  PubMed  Google Scholar 

  15. Rodriguez-Brito B, Li LL, Wegley L, Furlan M, Angly F, Breitbart M, et al. Viral and microbial community dynamics in four aquatic environments. Isme J. 2010;4(6):739–51.

    Article  PubMed  Google Scholar 

  16. Warwick-Dugdale J, Buchholz HH, Allen MJ, Temperton B. Host-hijacking and planktonic piracy: how phages command the microbial high seas. Virol J. 2019;16:13.

    Article  Google Scholar 

  17. Scanlan PD, Hall AR, Blackshields G, Friman VP, Davis MR, Goldberg J, et al. Coevolution with bacteriophages drives genome-wide host evolution and constrains the acquisition of abiotic-beneficial mutations. Mol Biol Evol. 2015;32(6):1425–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Roossinck MJ. The good viruses: viral mutualistic symbioses. Nat Rev Microbiol. 2011;9(2):99–108.

    Article  CAS  PubMed  Google Scholar 

  19. Emerson JB, Roux S, Brum JR, Bolduc B, Woodcroft BJ, Jang HB, et al. Host-linked soil viral ecology along a permafrost thaw gradient. Nat Microbiol. 2018;3(8):870–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Howard-Varona C, Hargreaves KR, Abedon ST, Sullivan MB. Lysogeny in nature: mechanisms impact and ecology of temperate phages. Isme J. 2017;11(7):1511–20.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Zhang JY, Gao Q, Zhang QT, Wang TX, Yue HW, Wu LW, et al. Bacteriophage-prokaryote dynamics and interaction within anaerobic digestion processes across time and space. Microbiome. 2017;5:10.

    Article  Google Scholar 

  22. Walker CB, Stolyar S, Chivian D, Pinel N, Gabster JA, Dehal PS, et al. Contribution of mobile genetic elements to Desulfovibrio vulgaris genome plasticity. Environ Microbiol. 2009;11(9):2244–52.

    Article  CAS  PubMed  Google Scholar 

  23. Maslov S, Sneppen K. Well-temperate phage: optimal bet-hedging against local environmental collapses. Sci Rep-UK. 2015;5(1):10523.

    Article  CAS  Google Scholar 

  24. Xu JW, Liu C, Hsu PC, Zhao J, Wu T, Tang J, et al. Remediation of heavy hetal contaminated soil by asymmetrical alternating current electrochemistry. Nat Commun. 2019;10:2440.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Mielke HW, Gonzales CR, Powell ET, Laidlaw MAS, Berry KJ, Laidlaw MAS, et al. The concurrent decline of soil lead and children's blood lead in New Orleans. Proc Natl Acad Sci USA. 2019;116(44):22058–64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Wu G, Kang HB, Zhang XY, Shao HB, Chu LY, Ruan CJ. A critical review on the bio-removal of hazardous heavy metals from contaminated soils: issues progress eco-environmental concerns and opportunities. J Hazard Mater. 2010;174(1-3):1–8.

    Article  CAS  PubMed  Google Scholar 

  27. Liu XY, Bai ZK, Shi HD, Zhou W, Liu XC. Heavy metal pollution of soils from coal mines in China. Nat Hazards. 2019;99(2):1163–77.

    Article  Google Scholar 

  28. Zou YD, Wang XX, Khan A, Wang PY, Liu YH, Alsaedi A, et al. Environmental remediation and application of nanoscale zero-valent iron and its composites for the removal of heavy metal ions: a review. Environ Sci Technol. 2016;50(14):7290–304.

    Article  CAS  PubMed  Google Scholar 

  29. Ma Y, Zhong H, He ZG. Cr (VI) reductase activity locates in the cytoplasm of Aeribacillus pallidus BK1 a novel Cr(VI)-reducing thermophile isolated from Tengchong geothermal region China. Chem Eng J. 2019;371:524–34.

    Article  CAS  Google Scholar 

  30. Gao SM, Schippers A, Chen N, Yuan Y, Zhang MM, Li Q, et al. Depth-related variability in viral communities in highly stratified sulfidic mine tailings. Microbiome. 2020;8(1):89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Trubl G, Jang HB, Roux S, Emerson JB, Solonenko N, Vik DR, et al. Soil viruses are underexplored players in ecosystem carbon processing. Msystems. 2018;3(5):21.

    Article  Google Scholar 

  32. Yu PF, Mathieu J, Yang Y, Alvarez PJJ. Suppression of enteric bacteria by bacteriophages: importance of phage polyvalence in the presence of soil bacteria. Environ Sci Technol. 2017;51(9):5270–8.

    Article  CAS  PubMed  Google Scholar 

  33. Yu PF, Mathieu J, Lu GW, Gabiatti N, Alvarez PJ. Control of antibiotic-resistant bacteria in activated sludge using polyvalent phages in conjunction with a production host. Environ Sci Tech Let. 2017;4(4):137–42.

    Article  CAS  Google Scholar 

  34. Yang JX, Wei W, Pi SS, Ma F, Li A, Wu D, et al. Competitive adsorption of heavy metals by extracellular polymeric substances extracted from Klebsiella sp J1. Bioresource Technol. 2015;196:533–9.

    Article  CAS  Google Scholar 

  35. Dogan NM, Kantar C, Gulcan S, Dodge CJ, Yimaz BC, Mazmanci MA. Chromium (VI) bioremoval by Pseudomonas bacteria: role of microbial exudates for natural attenuation and biotreatment of Cr (VI) contamination. Environ Sci Technol. 2011;45(6):2278–85.

    Article  CAS  PubMed  Google Scholar 

  36. Ng TW, Cai QH, Wong CK, Chow AT, Wong PK. Simultaneous chromate reduction and azo dye decolourization by Brevibacterium casei: azo dye as electron donor for chromate reduction. J Hazard Mater. 2010;182(1-3):792–800.

    Article  CAS  PubMed  Google Scholar 

  37. Anderson C, Jakobsson AM, Pedersen K. Influence of in situ biofilm coverage on the radionuclide adsorption capacity of subsurface granite. Environ Sci Technol. 2007;41(3):830–6.

    Article  CAS  PubMed  Google Scholar 

  38. Kumar R, Singh R, Kumar N, Bishnoi K, Bishnoi NR. Response surface methodology approach for optimization of biosorption process for removal of Cr (VI), Ni (II) and Zn (II) ions by immobilized bacterial biomass sp Bacillus brevis. Chem Eng J. 2009;146(3):401–7.

    Article  CAS  Google Scholar 

  39. Lamont I, Richardson H, Carter DR, Egan JB. Genes for the establishment and maintenance of lysogeny by the temperate coliphage 186. J Bacteriol. 1993;175(16):5286–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Villafane R, Black J. Identification of 4 genes involved in the lysogenic pathway of the Salmonella newington bacterial-virus epsilon (34). Arch Virol. 1994;135(1-2):179–83.

    Article  CAS  PubMed  Google Scholar 

  41. Yin H, He BY, Peng H, Ye JH, Yang F, Zhang N. Removal of Cr (VI) and Ni (II) from aqueous solution by fused yeast: study of cations release and biosorption mechanism. J Hazard Mater. 2008;158(2-3):568–76.

    Article  CAS  PubMed  Google Scholar 

  42. Nies DH. Efflux-mediated heavy metal resistance in prokaryotes. Fems Microbiol Rev. 2003;27(2-3):313–39.

    Article  CAS  PubMed  Google Scholar 

  43. Miao Y, Liao RH, Zhang XX, Wang Y, Wang Z, Shi P, et al. Metagenomic insights into Cr (VI) effect on microbial communities and functional genes of an expanded granular sludge bed reactor treating high-nitrate wastewater. Water Res. 2015;76:43–52.

    Article  CAS  PubMed  Google Scholar 

  44. Malaviya P, Singh A. Metagenomic insights into Cr (VI) effect on microbial communities and functional genes of an expanded granular sludge bed reactor treating high-nitrate wastewater. Water Res. 2015;76:43-52.

  45. Martinez-Garcia M, Santos F, Moreno-Paz M, Parro V, Anton JL. Unveiling viral-host interactions within the ‘microbial dark matter’. Nat Commun. 2014;5:8.

    Article  Google Scholar 

  46. Forterre P. The virocell concept and environmental microbiology. Isme J. 2013;7(2):233–6.

    Article  CAS  PubMed  Google Scholar 

  47. Breitbart M, Bonnain C, Malki K, Sawaya NA. Phage puppet masters of the marine microbial realm. Nat Microbiol. 2018;3(7):754–66.

    Article  CAS  PubMed  Google Scholar 

  48. Dhal B, Thatoi HN, Das NN, Pandey BD. Chemical and microbial remediation of hexavalent chromium from contaminated soil and mining/metallurgical solid waste: A review. J Hazard Mater. 2013;250:272–91.

    Article  PubMed  Google Scholar 

  49. Liu YR, Delgado-Baquerizo M, Bi L, Zhu J, He JZ. Consistent responses of soil microbial taxonomic and functional attributes to mercury pollution across China. Microbiome. 2018;6:12.

    Article  Google Scholar 

  50. Shapiro OH, Kushmaro A, Brenner A. Bacteriophage predation regulates microbial abundance and diversity in a full-scale bioreactor treating industrial wastewater. Isme J. 2010;4(3):327–36.

    Article  PubMed  Google Scholar 

  51. Liang XL, Zhang YY, Wommack KE, Wilhelm SW, DeBruyn JM, Sherfy AC, et al. Lysogenic reproductive strategies of viral communities vary with soil depth and are correlated with bacterial diversity. Soil Biol Biochem. 2020;144:10.

    Article  Google Scholar 

  52. Touchon M, de Sousa JAM, Rocha EPC. Rocha Embracing the enemy: the diversification of microbial gene repertoires by phage-mediated horizontal gene transfer. Curr Opin Microbiol. 2017;38:66–73.

    Article  CAS  PubMed  Google Scholar 

  53. Laslett D, Canback B. ARAGORN, a program to detect tRNA genes and tmRNA genes in nucleotide sequences. Nucleic Acids Res. 2004;32(1):11–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Yu PF, Mathieu J, Li MY, Dai ZY, Alvarez PJJ. Isolation of polyvalent bacteriophages by sequential multiple-host approaches. Appl Environ Microb. 2016;82(3):808–15.

    Article  CAS  Google Scholar 

  55. Feiner R, Argov T, Rabinovich L, Sigal N, Borovok I, Herskovits AA. A new perspective on lysogeny: prophages as active regulatory switches of bacteria. Nat Rev Microbiol. 2015;13(10):641–50.

    Article  CAS  PubMed  Google Scholar 

  56. Wang XX, Kim Y, Ma Q, Hong SH, Pokusaeva K, Sturino JM, et al. Cryptic prophages help bacteria cope with adverse environments. Nat Commun. 2010;1:9.

    Article  CAS  Google Scholar 

  57. Moon K, Jeon JH, Kang IN, Park KS, Lee KY, Cha CJ, et al. Freshwater viral metagenome reveals novel and functional phage-borne antibiotic resistance genes. Microbiome. 2020;8(1):5.

    Article  Google Scholar 

  58. Oliver KM, Degnan PH, Hunter MS, Moran NA. Bacteriophages encode factors required for protection in a symbiotic mutualism. Science. 2009;325(5943):992–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Lucas PPB, Aymé S, Witold K, Marie-Christine B, Lars HH, João CS, Laurent P. Impact of phages on soil bacterial communities and nitrogen availability under different assembly scenarios. Microbiome. 2020;8(1):52.

  60. Rowell DL. Soil science: methods & applications. Harlow: Longman Scientific & Technical; 1994.

    Google Scholar 

  61. EPA US. Alkaline digestion for hexavalent chromium (EPA Method 3060a). U.S: EPA; 1996.

    Google Scholar 

  62. Patel A, Noble RT, Steele JA, Schwalbach MS, Hewson I, Fuhrman JA. Virus and prokaryote enumeration from planktonic aquatic environments by epifluorescence microscopy with SYBR Green I. Nat Protoc. 2007;2(2):69–276.

    Article  Google Scholar 

  63. Chen SF, Zhou YQ, Chen YR, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):884–90.

    Article  Google Scholar 

  64. Li DH, Luo RB, Liu CM, Leung CM, Ting HF, Sadakane K, et al. MEGAHIT v1, 0: a fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods. 2016;102:3–11.

    Article  CAS  PubMed  Google Scholar 

  65. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34(18):3094–3100.

  66. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The Sequence Alignment/Map format and SAMtools, Bioinformatics. 2009;25(16):2078–9.

  67. Huson DH, Auch AF, Qi J, Schuster SC. MEGAN analysis of metagenomic data. Genime Res.. 2007;17(3)377–86..

  68. Noguchi H, Taniguchi T, Itoh T. Meta Gene Annotator: detecting species-specific patterns of ribosomal binding site for precise gene prediction in anonymous prokaryotic and phage genomes. DNA Res. 2008;15(6):387–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Li WZ, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22(13):1658–9.

    Article  CAS  PubMed  Google Scholar 

  70. Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40(D1):D290–301.

    Article  CAS  PubMed  Google Scholar 

  71. Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J Mol Biol. 2016;428(4):726–31.

    Article  CAS  PubMed  Google Scholar 

  72. Huerta-Cepas J, Szklarczyk D, Heller D, Hernandez-Plaza A, Forslund SK, Cook H, et al. eggNOG 5, 0: a hierarchical functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2019;47(D1):D309–14.

    Article  CAS  PubMed  Google Scholar 

  73. Menzel P, Ng KL, Krogh A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat Commun. 2016;7:9.

    Article  Google Scholar 

  74. Adriaenssens EM, Kramer R, Van Goethem MW, Makhalanyane TP, Hogg I, Cowan DA. Environmental drivers of viral community composition in Antarctic soils identified by viromics. Microbiome. 2017;5:14.

    Article  Google Scholar 

  75. Yu DT, He JZ, Zhang LM, Han LL. Viral metagenomics analysis and eight novel viral genomes identified from the Dushanzi mud volcanic soil in Xinjiang China. J Soil Sediment. 2019;19(1):81–90.

    Article  CAS  Google Scholar 

  76. Fuhrman JA. Marine viruses and their biogeochemical and ecological effects. Nature. 1999;399(6736):541–8.

    Article  CAS  PubMed  Google Scholar 

  77. Roux S, Hallam SJ, Woyke T, Sullivan MB. Viral dark matter and virus-host interactions resolved from publicly available microbial genomes. Elife. 2015;4:20.

    Article  Google Scholar 

  78. Roux S, Enault F, Hurwitz BL, Sullivan MB. VirSorter: mining viral signal from microbial genomic data. Peer J. 2015;3:20.

    Article  Google Scholar 

  79. De Sanctis M, Fanelli G, Gjeta E, Mullaj A, Attorre F. The forest communities of Shebenik-Jabllanice National Park (Central Albania). Phytocoenologia. 2018;48(1):51–76.

    Article  Google Scholar 

Download references


We thank Dr. Jinfeng Wang from Beijng Institute of life Science, China Academy of Sciences and Dr. Yu Shi from Institute of Soil Science, Chinese Academy of Sciences for the experimental design and valuable feedback.


This work was financially supported by National Key Research and Development Program of China (2018YFC1803100), NSF ERC on Nanotechnology-Enabled Water Treatment (EEC-1449500), the Youth Innovation Promotion Association of the Chinese Academy of Sciences (2018350), and the Jiangsu Provincial Natural Science Funds for Distinguished Young Scholar (BK20180110).

Author information

Authors and Affiliations



All authors contributed intellectually to and agreed to this submission. HD, YPF, and YM designed the experiments. HD conducted the experiments and collected the data. HD, YPF, and YM analyzed the data. HD, YPF, YM, and CS wrote the initial draft of the manuscript, while XJ and PA provided substantial feedback. The authors read and approved the final manuscript.

Corresponding authors

Correspondence to Pingfeng Yu or Mao Ye.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no conflicts of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Physical chemical parameters of the soil samples. Table S2. Integrase family genes (annotated using PFAM). Table S3. Membrane transporter genes (annotated using KEGG). Table S4. Reductase genes (annotated using KEGG). Table S5. Specific information on protein annotation of four lysogenic phages.

Additional file 2: Fig. S1.

Geographic information of contaminated sites and sampling locations. Fig. S2. Species richness of bacterial communities indicated by Chao1 index (A) and ACE index (B). Fig. S3. The phylogenetic tree of dominant bacterial species with relative abundance greater than 0.1%. Fig. S4. Profile of viral communities in Cr-contaminated soils. Fig. S5. This heat map indicates the distribution of top 20 viruses in different samples. Fig. S6. Redundant analysis (RDA) of virome. Fig. S7. The habitat composition of RefSeq viruses linked by viral contigs from Cr-contaminated soil. Fig. S8. Predicted virus-host linkages at the phylum level. Fig. S9. Predicted virus-host linkages at the genus level. Fig. S10. The non-metric multidimensional scaling (NMDS) of viral functional genes. Fig. S11. The relative abundance of virus-carried genes annotated by eggNOG database. Fig. S12. Epifluorescence-microscopy image of viruses.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Huang, D., Yu, P., Ye, M. et al. Enhanced mutualistic symbiosis between soil phages and bacteria with elevated chromium-induced environmental stress. Microbiome 9, 150 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: