General characterization of soil profiles
In this study, contaminated soil profiles were sampled from two pollution conditions (i.e., CP1 and CP2 for more contaminated sites and UCP for an uncontaminated site). The experimental design and sample collection details are summarized in Figure S1. A statistical comparison of three profiles was conducted by a one-way ANOVA (Figure 1A). Overall, the contaminated sites (CP1 and CP2) were characterized by low pH, high Eh, and high concentrations of sulfate, ferrous (Fe(II)) iron, ferric (Fe(III)) iron, Sbtot, and Astot (p < 0.01). Additionally, CP1 and CP2 formed a similar cluster pattern of geochemical parameters and distinctly separated from UCP (Figure 1B). In the following analyses, CP1 and CP2 were grouped as “CP” due to their similar soil properties. At each level, CP demonstrated a significantly lower pH, but higher Eh and sulfate concentrations in all profile depths (p < 0.0001) (Figure 2). TOC contents and Fetot concentrations were significantly different in deep soils (1.5–2 m, p < 0.001) and surface (0–0.5 m, p < 0.01), respectively, while nitrate and Fe(II) concentrations were not significantly different at any depth.
The depth-resolved differences of Sb and As concentrations in each profile were also compared (Figure 3). In surface soil, CP was severely contaminated by smelting actives and the Sbtot and Astot concentrations were much higher than in UCP. The average values of Sbtot concentration in CP was 10–20 fold higher than in UCP. The highest value of Sbtot concentrations in CP was 815 mg/kg (0–0.5 m), which was much higher than the secondary standard of Environmental Quality Standards for Soils in China (10 mg/kg; GB15168-2008). The bioaccessible fractions (-exe and -srp) of Sb in surface soil of CP were also higher than UCP. Likewise, the highest Astot concentration in CP was 70.8 mg/kg, which also exceeded the natural background values of agricultural land (1–40 mg/kg; GB15168-2018). The bioavailable fractions (-exe and -srp) of As in surface soil were not detectably different between UCP and CP.
As soil depths increased, Sbtot and Astot dramatically decreased to 141.9 and 38.1 mg/kg in deep soil (1.5–2 m), respectively. Similar descending trends from surface to deep soil were also detected in the -exe and -srp fractions of Sb and As. The concentrations of most Sb and As fractions significantly (p < 0.01) negatively correlated with the depths of the vertical profiles (except Assrp, which was not significant). Results suggested that Sb and As contamination mainly originated from smelting deposition in surface layer (< 0.5 m), with a downward migration in the vertical profiles (e.g., atmospheric deposition) . Consequently, the effects of contamination gradually declined as soil depth increased. As expected, different contamination conditions at each soil depth had different effects on soil microorganisms.
Composition of microbial communities at different profile depths
The compositions of microbial communities were characterized by 16S rRNA amplicon sequencing. The Bray–Curtis-based principal coordinate analysis (PCoA) revealed that the distribution of microbial communities at each profile depth was strongly affected by Sb and As contamination in CP (Figure S2). Smelting contamination reduced the alpha diversity of microbial communities at each profile (Figure 4A and B). Random forest model suggested the low alpha diversity could be attributed to the contents of Sbtot and Astot (Figure 4C). The relative abundance of the top 10 phyla detected in the soil profiles are presented in Figure S3. The taxonomy analyses indicated that the soil microbiota in CP at different depths were dominated by Acidobacteria (18.8%–44.7%), Chloroflexi (8.7%–42.4%), and Proteobacteria (11.4%–27.1%). The distribution of these predominant phyla was not similar to those in UCP. Several phyla were more abundant in UCP that were not detected in CP throughout the soil profiles, including Planctomycetes, Bacteroidetes, Gemmatimonade, and Nitrospirae.
This study also characterized the depth-resolved microbial communities and found that many phyla exhibited different dynamics throughout soil depths (Figure 5). A previous study demonstrated that surface soil microbial communities are naturally distinct from deep soil communities in their composition and structure . Consistently, the relative abundance of the predominant phyla in UCP, Acidobacteria, Actinobacteria, and Bacteroidetes, gradually decreased from surface soil (0–0.5 m) to deep soil (1.5–2 m), while the relative abundance of Gemmatimonadetes and Nitrospirae dramatically increased as soil depth increased. However, the depth-resolved microbial community compositions in CP were distinct compared to UCP. As soil depth increased, Chloroflexi and Thaumarchaeota were the most variable phyla. Their relative abundances gradually increased from surface to deep soil, ranging from 8.69% to 42.43% and 0.49% to 20.17%, respectively, but both phyla were found in low abundance in UCP at each depth. In contrast, the relative abundance of Acidobacteria and Proteobacteria in CP gradually decreased as soil depth increased. The distinct cluster pattern of the microbial communities (Figure S2, S3), as well as the highly variable abundance of microorganisms (Figure 5), indicated that the variation of Sb and As contamination along the vertical soil profiles shaped the predominant microbial communities and had a highly significant effect on the depth-resolved microbial communities in soils. The co-occurrence analysis supported this hypothesis (Figure S4). It was also observed that Sbtot was one of the most important environmental drivers that affected soil microbiota due to its highly connected properties.
Microbial interaction networks at different profile depths
The above results provide important insight into how single members of soil microbial communities respond to soil properties. Therefore, the multitude of direct and indirect interactions that occur at the whole microbial community level were further investigated. Here, an applied co-occurrence network analysis uncovered a 16S rRNA gene data set containing > 4,800,000 sequences from all soil profile samples (n = 60) and calculated robust associations between microbial taxa. A total of eight co-occurrence networks were constructed (four CP networks and four UCP networks) to compare microbial interactions at different profile depths (Figure 6A). In UCP, the pattern of microbial interactions remained relatively stable from surface soil (0–0.5 m) to deep soil (1.5–2 m), as indicated by the connections between the nodes and edges. In contrast, the microbial interactions in CP showed an obvious changing pattern as soil depth increased. The potential interactions in deeper networks (less contaminated) were stronger than surface networks (heavily contaminated).
Results also revealed that there were remarkable differences in the topological parameters of the networks (Figure 6B and 6C; Table S1). The node number (representing keystone microbiomes) and cluster number (representing function modules of microbial communities) in CP were always lower than UCP, indicating weak interactions among soil microbiota in contaminated profiles. In particular, the node number in CP gradually increased from surface to deep soil, indicating that individual species became more connected with other patterns to perform syntrophic functions in deep soils than those in surface soil. Therefore, the cluster number of microbial communities gradually decreased and grouped into less modules. Additionally, the ratio of positive/negative edge numbers of each network was calculated (Figure 6D and 6E; Table S1). Both positive and negative edge number ratios in UCP were ~50%. As a comparison, Sb and As contamination reduced the proportion of negative correlations in microbial co-occurrence networks, as the positive edge ratio was much higher than the negative edge in CP (70% vs. 30%). A previous modeling study suggested that communities with a large proportion of members connected through positive links were unstable. In such communities, its members may respond in tandem to environmental fluctuations, resulting in positive feedback and co-oscillation. In contrast, negative links may stabilize co-oscillation in these communities and promote network stability . This theory is supported by the observations of this study, where the networks were loosely connected in heavily contaminated surface soil but gradually recovered and were well connected in the less contaminated deep soil.
A network analysis of the significant taxon co-occurrence patterns, beyond the basic inventory descriptions of the composition and diversity of the microbial communities, may help explore the structure of complex microbial communities across spatial and temporal gradients [23, 24]. We hypothesized that Sb and As concentrations, rather than soil resource availability, is the main factor responsible for the observed changes in microbial interactions throughout the vertical soil profiles. Previous studies demonstrated that Sb and As contamination had considerable effects on soil microbial community dynamics [1, 25, 26]. However, understanding how these complex microbial communities recover from disturbances in situ remains a major challenge. In this study, different components of the microbial communities responded differently to Sb and As contamination at different soil depths. As soil depth increased, Sb and As contamination significantly declined from surface to deep soil (p < 0.0001, Figure 3). Soil microbiota in CP were generally more resilient, but less resistant, than UCP. Thus, the potential interactions among soil microbiota gradually shifted from weak (in heavily contaminated surface soil) to strong (in less contaminated deep soil). This distinction between surface and deep soil can be associated with the direct or indirect effects of Sb and As contamination.
Sb and As contamination may influence C and N cycling
A total of 12 metagenomes from CP and 2 metagenomes from UCP were constructed to investigate microbial metabolic functions (Figure S1). Based on the shotgun metagenomic sequencing, further in-depth investigations of the fundamental biogeochemical processes uncovered how soil microbiota adapted to different levels of Sb and As contamination throughout soil profiles. Results revealed that microbial metabolic potential changed with soil depth. In the following sections, the genes involved in the C and N metabolism pathways were analyzed against the KEGG database.
C metabolism pathway: Generally, the relative abundances of the genes encoding the Arnon–Buchanan reductive citrate cycle (rTCA) pathway in both UCP and CP were much higher than the five carbon fixation pathways (Figure 7A). This observation may be attributed to the high energy efficiency of rTCA, as 2 ATP equivalents are required to form pyruvate, while > 5 ATP equivalents are required in other C fixation pathways, except the Wood–Ljungdahl pathway, which requires < 1 ATP equivalent [27, 28]. Additionally, most genes involved in C metabolism from CP samples were much lower than UCP throughout all soil profiles. Since C bioavailability is one of the most important determinants for microbial activity , the low abundances of these genes indicated microbial growth at each depth was strongly restricted by smelting contamination.
Besides, several pathways in CP gradually increased from surface to deep top soil, including the rTCA, dicarboxylate-hydroxybutyrate cycle (DC/4-HB), and hydroxypropionate-hydroxybutylate cycle (3-HP/4-HB) pathways. However, a previous study reported that C metabolic potential and diversity decreased from surface to deep soil (negatively correlated with soil depth), as surface soil harbored a higher proportion of pre-adapted inhabitants for substrate metabolism . Additionally, surface soils are rich in available C substrates due to the input of root exudates, surface litter, and root detritus. Because the quantity and quality of C substrates decline in deeper soil layers, the rates of C input at lower depths are generally low and C tends to be in limited availability .
Such discrimination may be explained by C cycling that was strongly inhibited by Sb and As contamination at each soil depth in CP, especially surface soil. As soil depth increased, Sb and As concentrations declined and the microbes involved in carbon cycling gradually recovered. These findings indicate that the effects of contamination were much greater than C availability in surface soil. This could be supported by the co-occurrence network analysis, where 4 major C metabolism pathways, rTCA, reductive pentose phosphate cycle (Calvin-Benson-Bassham cycle; CBB), 3-HP/4-HB, and DC/4-HB, were all negatively correlated with Sbtot and Astot concentrations (Figure 7B).
N metabolism pathway: Microbial-mediated N transformations are important for N cycling in the soil. Genes encoding the assimilatory nitrate reduction (nitrate → ammonia) and denitrification (nitrate → nitrogen) pathways were more abundant in CP samples (Figure S5). Although the relative abundances of most genes did not directly correlate with soil depths, nitric oxide reductase (norB) and nitrate/nitrite transporter (narK) positively correlated with Astot or Sbtot, suggesting that smelting contamination may promote the denitrification pathway. The enrichment of the assimilatory nitrate reduction and denitrification pathways reduce nitrate concentrations. During this process, nitrate serves as an electron acceptor and couples with As(III) for anaerobic oxidization . These findings provide insight into microbial adaptation in As contaminated soils, as toxicity was largely reduced from As(III) to As(V). In contrast, genes encoding the nitrification (ammonia → nitrite) pathway were restricted by smelting contamination, as ammonia monooxygenase (pmoABC) was much lower in CP samples and gradually increased from surface (0–0.5 m) to deep soil (1.5–2 m).
Microbial resistance to Sb and As contamination
The depth-resolved microbial activity characterization indicated that innate microbiota was greatly affected by the contamination of Sb and As. Because As and Sb have similar chemical structures, microbes may use similar metabolic pathways to transform both As and Sb .
As resistance genes in UCP versus CP: Results revealed that genes related to the arsenic resistance transcriptional regulator (arsR) were more abundant than other genes (Figure 8A). A set of ars genes [encoding As(V) reduction] were enriched by contamination and widely detected in CP, including arsR, arsB, arsC, ACR3, and arsH. These genes were previously found to be prevalent among sediment, activated sludge, and unpolluted soil [32, 33]. However, the genes encoding As(III) oxidation (i.e., aoxAB) were quite low in both UCP and CP. Genes encoding As(V) respiration reduction (arr genes) and As(III) methylation (arsM) were undetectable. The prevalence of ars genes could be ascribed to the versatility of As(V) reducing bacteria . For example, one study found a broad phylogenetic distribution of arsC genes in soil bacteria, including Rhizobiales, Burkholerials, and Geobacillus . In addition to ars genes, the relative abundance of arsenite-transporting ATPase (ASNA1) and arsenite transporter, the ACR3 family (ACR3) also increased in CP. Four gene types were solely detected in CP samples, including arsenical resistance protein (arsH), arsenical-resistance protein 2 (ACR2), arsenite methyltransferase (AS3MT), and arsenate-mycothiol transferase (arsC) (Figure 8B). The relatively high abundances of these genes in CP suggest that bacterial resistance in response to Sb and As contamination was enhanced , including the ability to confer resistance to arsenate and arsenite via operons (e.g., arsR, arsB, arsC) (Ji and Silver, 1992).
As resistance genes at different profile depths: The distributions of these As metabolism genes were highly variable throughout soil profiles in CP. The abundance of these genes decreased with soil depth (e.g., arsB, ACR3, and arsC, Figure 8A and B). The enrichment of these genes in surface soil may be attributed to the higher concentrations of As and Sb in top soil layers. Elevated As and Sb concentrations exerted strong natural selection on the present microbiota and, in turn, these microbiotas expressed more functional genes that encode the As detoxification process to relieve such perturbations. In contrast, several other genes increased with soil depth (e.g., arsR). The detection of highly abundant As resistance genes in deep soil layers implied that Sb and As contamination affected the whole soil profile, even at a depth of 2 m. Thus, these findings indicate that Sb and As pollution treatments should not be limited to surface soil remediation and that the distribution, migration, and toxicity of Sb and As pollution in deep layers should also be considered.
Correlation of As resistance genes and contaminations: In addition to the highly variable abundances of As-related genes in CP, strong correlations were detected among these genes and concentrations of As fractions (Figure 8C), such as ACR3 and arsB. A previous study detected a positive correlation between the total abundances of As-related genes and Astot concentrations . A separate study found that the abundances of genes involved in different As biotransformation processes also positively correlated with different As species in soil samples . Furthermore, in this study, strong correlations of these As-related genes and Sb fractions were observed. For example, Sbtot positively correlated with arsB and ACR3, suggesting the potential roles of these genes in Sb metabolism and contamination may promote the proliferation of such genes. Because As and Sb share structural similarities, microorganisms may use similar metabolic pathways to transform both As and Sb. A previous study reported that As and Sb induced arsenic-resistant operons, which contain an arsenical pump membrane protein (arsB) . Similarly, YL Meng, Z Liu and BP Rosen  reported that the ArsB family can catalyze the transport of Sb fractions. Genomic sequencing of Comamonas sp. S44 (an Sb-oxidizing bacterium) revealed the presence of genes encoding ArsB . Therefore, it is proposed that Sb redox may be catalyzed by enzymes that are encoded by the same As-related functional genes. However, due to the limited database of genes encoding Sb cycling and resistance in KEGG functional categories, this hypothesis requires further investigation based on cultured isolations of Sb-metabolizing bacteria.