ABSTRACT
A lack of knowledge of the microbial responses to environmental change at the species and functional levels hinders our ability to understand the intrinsic mechanisms underlying the maintenance of microbial ecosystems. Here, we present results from temporal microcosms that introduced inorganic and organic contaminants into agro-soils for 90 days, with three common legume plants. Temporal dynamics and assemblage of soil microbial communities and functions in response to contamination under the influence of growth of different plants were explored via sequencing of the 16S rRNA amplicon and by shotgun metagenomics. Soil microbial alpha diversity and structure at the taxonomic and functional levels exhibited resilience patterns. Functional profiles showed greater resilience than did taxonomic ones. Different legume plants imposed stronger selection on taxonomic profiles than on functional ones. Network and random forest analyses revealed that the functional potential of soil microbial communities was fostered by various taxonomic groups. Betaproteobacteria were important predictors of key functional traits such as amino acid metabolism, nucleic acid metabolism, and hydrocarbon degradation. Our study reveals the strong resilience of the soil microbiome to chemical contamination and sensitive responses of taxonomic rather than functional profiles to selection processes induced by different legume plants. This is pivotal to develop approaches and policies for the protection of soil microbial diversity and functions in agro-ecosystems with different response strategies from global environmental drivers, such as soil contamination and plant invasion.
IMPORTANCE Exploring the microbial responses to environmental disturbances is a central issue in microbial ecology. Understanding the dynamic responses of soil microbial communities to chemical contamination and the microbe-soil-plant interactions is essential for forecasting the long-term changes in soil ecosystems. Nevertheless, few studies have applied multi-omics approaches to assess the microbial responses to soil contamination and the microbe-soil-plant interactions at the taxonomic and functional levels simultaneously. Our study reveals clear succession and resilience patterns of soil microbial diversity and structure in response to chemical contamination. Different legume plants exerted stronger selection processes on taxonomic than on functional profiles in contaminated soils, which could benefit plant growth and fitness as well as foster the potential abilities of hydrocarbon degradation and metal tolerance. These results provide new insight into the resilience and assemblage of soil microbiome in response to environmental disturbances in agro-ecosystems at the species and functional levels.
INTRODUCTION
Agricultural ecosystems are experiencing increasing anthropogenic pressure and environmental perturbation, such as climate oscillation, plant invasion, and the accumulation of pollutants, pesticides, and antibiotics in soil (1–3). In particular, long-standing contamination is often associated with substantial changes in soil biodiversity (4, 5). In the terrestrial ecosystem, soil microorganisms play critical roles in driving the global biogeochemical cycles of carbon, nitrogen, and other inorganic elements (6), and they respond rapidly to environmental change caused by contamination (7, 8). Given their importance in ecosystem functioning and services, it is vital to determine the temporal dynamics of soil microbial communities and their functions in agro-ecosystems in response to such contamination. This realm of investigation could provide valuable insight into the restoration of polluted ecosystems and environmental management, and yet it remains understudied.
Plant invasion is an important biotic and environmental perturbation that can alter soil nutrient cycling at the global scale (9). Plant invasion may also influence the stability of soil microbial biodiversity by introducing colonized arbuscular mycorrhizal fungi, ectomycorrhizal fungi, and soilborne pathogens (10). Plant roots provide a considerable amount of nutrients to surrounding environments via their release of exudates and mucilage, which modify the local soil’s physiochemical properties, to subsequently shape the microbial communities residing in the soil (11–14). Moreover, specific microbes are able to assemble into plant-associated communities that influence terrestrial carbon and nutrient cycling, as well as host plants’ growth and health (15–17). Recently, host plant effects on the assembly of root-associated microbiomes were documented in Arabidopsis spp. (18, 19), rice (20), and some legumes (21). Microbes have evolved genes enabling them to adapt to plant environments; for example, more carbohydrate metabolism functions and fewer mobile elements characterized the plant-associated bacteria than with taxonomically related yet non-plant-associated ones (22). In addition, previous studies have demonstrated that soil contamination could alter the plant rhizosphere metatranscriptome, and genes related to hydrocarbon degradation were generally more expressed in contaminated soils; however, the exact complements of genes induced were different between bulk and rhizosphere soils (23, 24). Distinguishing them from other plants, legumes form nodules via symbiosis with N2-fixing rhizobia, generating specific associations between legumes and their root-associated microbiomes (21). A recent study showed that soybean roots imposed clear microbial selection at both the taxonomic and functional levels (25); the functional selection was related to the metabolism of nitrogen, iron, phosphorus, and potassium, all of which were enriched in the rhizosphere, thus potentially promoting plant growth and nutrition (25). Investigating the joint influence of soil and plant characteristics on the microbial community assemblage at the functional level is helpful for better understanding complex microbe-soil-plant interactions in agro-ecosystems, particularly in those with chemical contamination.
In the face of shifting soil conditions, microbes exhibit remarkable stability, largely due to their high degree of metabolic flexibility and physiological tolerance; their high abundances, widespread dispersal, and high growth rates; and their evolutionary adaptation via horizontal gene transfer (26, 27). Generally, the stability of microbial ecosystems is attributed to three mechanisms, resistance, resilience, and functional redundancy (27). In the first case, some microorganisms display a high degree of tolerance to disturbances; for example, those bacteria with resistance to heavy-metal or oil contamination (8). When the microbial community is changed by a disturbance yet rapidly recovers to its initial or alternative stable state, that is resilience (28, 29). When perturbed, the ecosystem processes retain similarity to those of the original state, even if the community is substantially altered without recovery, and this is then attributed to microbial functional redundancy (27). Soil functions show resilience to environmental change via acclimatization, as driven by three potential underlying mechanisms, as follows: (i) generalists with broad adaptations to new conditions are functionally dominant members of the community, (ii) specialists with distinct functional traits become active and predominate, and (iii) functionally dominant microbial members with rapid microbial adaptations are recruited under environmental disturbances (30). However, broad research interests in the microbial responses to environmental change largely focus on the whole-soil or community level, so we still lack much knowledge about species-level microbial responses and how they scale up to the functional level (30, 31).
The major aim of this study was to reveal the temporal dynamics of the soil microbial community in an agro-ecosystem in response to chemical contamination and to determine the microbial community assemblage occurring under the influence of different legume plants. We selected phenanthrene, n-octadecane, and cadmium as the pollutants, since they are prevalent in contaminated agro-soils. Three common legumes, Robinia pseudoacacia (woody), Medicago sativa (herbaceous), and Vicia villosa (herbaceous), were used in this study because of their wide distribution and suitability for a wide range of environments. We applied 16S rRNA gene amplicon sequencing to analyze the community structure and used the total DNA shotgun sequencing approach to assess their functional profiles. Our results demonstrated a strong resilience of the soil microbiome to chemical contamination and sensitive responses of taxonomic rather than functional profiles to selection processes induced by different legume plants, which provide valuable information for better understanding microbial resilience and assemblage patterns to environmental change at species and more deeply functional levels.
RESULTS
General characteristics of the sequencing data.In total, we obtained 21 soil samples contaminated with phenanthrene, n-octadecane, and CdCl2, including four time points (0, 10, 30, and 90 days), in three legume plant treatments with 90 days of growth (Table 1). Overall, the 16S rRNA gene amplicon sequencing yielded 1,045,206 high-quality sequences. The average number of sequences per sample (n = 21) was 49,772 (maximum, 59,291 sequences; minimum, 42,994 sequences; standard deviation [SD], 4,569 sequences). The total operational taxonomic unit (OTU) number was 4,577, defined at 97% sequence similarity. Proteobacteria (40.9%), Actinobacteria (27.7%), Acidobacteria (8.77%), and Chloroflexi (6.28%) accounted for the largest proportion of sequences. Meanwhile, the metagenomic shotgun sequencing yielded 256.79 Gb of data. After the quality trimming, a total of 164,030,284 sequences were obtained. The constructed metagenomic libraries were dominated by bacteria (80.98%), but there were also sequences matching to archaea (2.21%), viruses (1.22%), and eukaryotes (0.25%).
Microbial alpha-diversity characteristics for the 16S rRNA gene amplicon and the metagenomic shotgun sequencing data sets of soil samples
Temporal dynamics of soil microbial taxonomic and functional traits in response to contamination.During the incubation period, the alpha diversity for the taxonomic and functional traits showed clear successional patterns in response to the contamination, first decreasing and then increasing to basically the same as the day 0 samples, except for richness of functional traits (Fig. 1A and B). For beta diversity, a principal-coordinate analysis (PCoA) based on the Bray-Curtis dissimilarity revealed that the samples for different incubation time points formed distinct clusters in the ordination space (Fig. 1C and D), with significant differences found at both the taxonomic (RANOSIM = 1, P < 0.001; R2ADONIS = 0.8994, P < 0.001) and functional (RANOSIM = 0.9506, P < 0.001; R2ADONIS = 0.8133, P < 0.001) levels. Furthermore, to explore whether the soil microbial community recovered under disturbance from contamination, we estimated the dissimilarities in the taxonomic and functional profiles between the samples taken at the initial (day 0) and other time points (Fig. 1E and F). The dissimilarities significantly decreased from the day 10 to day 90 samples, which varied more rapidly for the functional than taxonomic profiles, as determined by the fitted quadratic ordinary least squares (OLS) models. In addition, the dissimilarities between the day 0 and day 90 samples were significantly lower for the functional than taxonomic profiles (P < 0.001), indicating a greater recovery of the functional profiles.
General patterns of microbial alpha and beta diversity in response to soil contamination and plant growth. The alpha diversity of richness and Shannon-Wiener index and the beta diversity based on Bray-Curtis distance between the samples were estimated. (A and B) Temporal changes in alpha diversity during incubation with organic and inorganic pollutants (phenanthrene + n-octadecane + cadmium) are shown for the taxonomic (A) and functional (B) traits. (C and D) PCoA of beta diversity among the samples of different incubation time points for the taxonomic (C) and functional (D) traits. (E and F) Dissimilarities of beta diversity between samples of day 0 and other time points for the taxonomic (E) and functional (F) traits, as estimated by the fitted quadratic OLS models. (G and H) PCoA of beta diversity among the samples of different plants grown in soil for the taxonomic (G) and functional (H) traits.
For the taxonomic profiles (Fig. 2), OTUs with an average relative abundance of >0.005% at all time points were classified into two clusters. Cluster 1 included 48 OTUs for which the abundance levels were first increased and then decreased; in contrast, cluster 2 included 884 OTUs that first decreased but then increased in abundance. Interestingly, the OTUs in both of these two clusters showed a trend of resiliency, with the lowest difference in abundance found between the day 0 and day 90 samples. The taxonomic distribution analysis showed that the OTUs in cluster 1 were mainly assigned to the genera Massilia, Lysobacter, Pseudoduganella, and Bacillus, while the genera Gaiella, Perlucidibaca, Sphingomonas, Nocardioides, and Aeromicrobium dominated in cluster 2.
Temporal dynamics of the microbial communities at the species level during incubation under soil contamination. This analysis was performed using the maSigPro method based on operational taxonomic units (OTUs) with an average relative abundance of >0.005%. (A) Temporal dynamic visualization of the significant OTUs was based on a cluster analysis that grouped OTUs with similar profiles and conveyed here in a heatmap. Each row in the heatmap has been standardized (to a mean of zero and a standard deviation of one), with its color intensity proportional to the standardized relative abundances of the taxa. The taxonomic distributions of the significant OTUs were estimated at the genus level for cluster 1 (B) and cluster 2 (C).
Similar to the taxonomic observations, two clusters were classified for the functional traits (see Fig. S1 and S2 in the supplemental material). In total, 855 Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology annotations (KOs) and 16 KEGG level 2 pathways were classified into cluster 1, whereas 322 KOs and 18 KEGG level 2 pathways were assigned to cluster 2. In cluster 1, the functional traits were primarily associated with signal transduction, energy metabolism, environmental adaptation, and the immune system. In cluster 2, the dominant functions were related to cell motility, some human diseases, amino acid metabolism, and the biodegradation and metabolism xenobiotics. Furthermore, we classified the KOs to their functional ontology as provided by Functional Ontology Assignments for Metagenomes (FOAM), which were relevant to environmental microorganisms. Based on FOAM level 2, totals of 55 and 57 traits were respectively grouped into clusters 1 and 2 (Fig. 3 and Table S1). The functions related to cellular response to stress, nucleic acid metabolism, saccharide and derivated synthesis, fatty acid oxidation, and the Embden-Meyerhof-Parnas (EMP) pathway were more prevalent in cluster 1 (Fig. 3A and B). Meanwhile, the functions related to transporters, hydrolysis of polymers, carbohydrate-active enzymes, and hydrocarbon degradation were most dominant in cluster 2 (Fig. 3A and C).
Temporal dynamics of the microbial communities at the functional level during incubation under soil contamination. This analysis was conducted using the maSigPro method based on the functional traits of FOAM level 2. (A) Temporal dynamic visualization of the significant OTUs was based on a cluster analysis that grouped OTUs with similar profiles and conveyed here in a heatmap. Each row in the heatmap has been standardized (to a mean of zero and a standard deviation of one), with its color intensity proportional to the standardized relative abundances of the taxa. (B and C) Functional distributions (Funs) of the significant traits were estimated at the FOAM level 1 for cluster 1 (B) and cluster 2 (C). Detailed information on the functional traits of FOAM level 2 is shown in Table S1.
Effect of plant selection processes on soil microbial taxonomic and functional traits.In the initial experimental design, we compared the day 10 and day 30 planted soils with the unplanted soils, finding no significant differences in microbial community composition (Fig. S3 and Table S2; by analysis of similarity [ANOSIM] and Adonis tests, P > 0.05). Therefore, we did not include day 10 and day 30 planted soils in the subsequent experiment and analysis. After 90 days of incubation, samples with and without plants formed distinct clusters in the ordination space (Fig. 1G and H), with significant differences in their taxonomic (RANOSIM = 0.6258, P = 0.003; R2ADONIS = 0.2571, P = 0.002) and functional (RANOSIM = 0.85, P = 0.005; R2ADONIS = 0.2284, P = 0.003) levels. Furthermore, the beta diversity for the taxonomic profiles was significantly different among the planted samples with different legumes (RANOSIM = 0.5391, P = 0.002; R2ADONIS = 0.3825, P = 0.002), while there were no significant differences detected among their functional profiles (RANOSIM = 0.1029, P = 0.197; R2ADONIS = 0.2859, P = 0.098). In addition, the alpha-diversity indices for taxonomic and functional traits were similar between the samples with and without plants after 90 days of incubation (P > 0.05).
To test whether the microbial populations were affected by the treatments, we identified the microbial taxa (OTUs) significantly enriched in each plant treatment relative to the unplanted soil using linear models. A bipartite association network visualized the associations between the significant OTUs and different plant species (Fig. 4). In clusters 1 to 3, a total of 62, 52, and 72 OTUs were associated with robinia, alfalfa, and vetch, respectively. In clusters 4 to 6, we found 46 OTUs simultaneously associated with two plant species. In cluster 7, there were 21 OTUs associated with three plant species. These results indicated a distinctly separate effect by different plants on the microbial communities, confirming the significant difference for beta diversity in the taxonomic profiles. Furthermore, Sphingopyxis, Phaselicystis, and Mesorhizobium were enriched in soil planted with robinia; Lysobacter, Flavisolibacter, Hydrogenophaga, and Bacillus dominated in soil planted with alfalfa; and Aeromicrobium, Marmoricola, Nitratireductor, and Rhizobium were more predominant in soil with vetch. The OTUs belonging to Pseudoxanthomonas, Devosia, and Chryseolinea were enriched by all the plant treatments.
The bipartite association network depicting the associations between the significantly enriched OTUs and the different plant species treatments (estimated via linear statistics in the “limma” package for R). Node sizes represent the relative abundance levels of the OTUs. Edges represent the association patterns of individual OTUs with different plant treatments. Circle-shaped nodes represent those OTUs only associated with one treatment. Diamond-shaped nodes represent the OTUs associated with two treatments. Triangle-shaped nodes represent the cross-combination OTUs associated with all three treatments. The number of OTUs and relative abundance levels are provided for each cluster (1 to 6), as are the taxonomic distributions of the dominant OTUs per cluster.
Some significant differences between the samples with and without plants were observed for the functional traits of FOAM level 1 (Fig. 5). The functions related to saccharide and derivative synthesis and the EMP pathway were more prevalent in the planted soils, while it was amino acid utilization biosynthesis metabolism, hydrocarbon degradation, and methanogenesis that were significantly overrepresented in the unplanted soils. A more detailed investigation was pursued for the functional traits based on FOAM level 2 (Fig. S4). Here, we found that functions associated with the response to nitrosative stress, the cellular response to oxidative stress, and pyrimidine and purine metabolisms were enriched in the planted soils, while phenylalanine, tyrosine and tryptophan biosynthesis, and glycosyltransferase functions predominated in the unplanted soils.
Effects of plant growth on the soil microbial community functional profile. (A) Relative abundance of the functional traits (FOAM database level 1) between samples with and without plants growing in the soil. The error bars show standard errors. Blue asterisks (*) indicate those categories significantly more abundant in samples with plant growth (P < 0.05, Wilcoxon rank-sum test) and orange asterisks (*) indicate the categories that were significantly more abundant in samples without plant growth. (B) PCoA plots of the functional genes classified into particular categories of FOAM database level 1. Similarity values between the samples with and without plant growth were examined via the ANOSIM test, which are shown in each plot. Only significant categories are displayed. The dashed ellipses in blue and orange represent the clustering of samples with and without plant growth, respectively.
Potential associations between soil microbial taxonomic and functional traits.To disentangle the potential associations between the microbial taxa and several important functional traits, we applied the network analysis based on the SpiecEasi method. At FOAM level 1, the network consisted of 335 nodes and 492 edges (Fig. S5). Most of these nodes belonged to Proteobacteria (n = 109), Actinobacteria (n = 94), and Acidobacteria (n = 39). The maximum number of edges (n = 42) occurred between Proteobacteria and Acidobacteria, and only four edges connected Actinobacteria to Acidobacteria. Nodes belonging to Proteobacteria were mostly connected with other phyla (138 edges), and the within-phylum associations were maximal for Actinobacteria, with 99 edges, followed by Proteobacteria, with 89 edges.
To explore the associations between microbial taxonomic and functional traits, we generated a subnetwork that only consisted of the connections involving functional traits (Fig. S6). Overall, 45 nodes (this included 24 OTUs and 21 functional traits) and 52 edges were selected for this subnetwork. There were 28 edges between the taxonomic and functional traits. For example, positive correlations were found for Limnobacter and Blastocatella and the tricarboxylic acid (TCA) cycle, Nocardioides and the EMP pathway, Solirubrobacterales and transporters; and Phenylobacterium and fatty acid oxidation. In contrast, Massilia was negatively correlated with the utilization of sugar and thiosulfate metabolism.
To further investigate the environmentally associated functions, we selected a functional subset (FOAM level 2) of cellular response to stress, hydrocarbon degradation, transporters, and the nitrogen cycle and examined their associations with the microbial taxa in soils (Fig. 6). Some interesting results were obtained. For example, Massilia was positively correlated with the major facilitator superfamily; Gaiellales was positively correlated with the ABC transporters; Gemmatimonadaceae was positively correlated with polycyclic aromatic hydrocarbon degradation; Gaiella and Acidobacteria were positively correlated with the regulation of translational initiation in response to osmotic stress; and Actinoplanes, Brevundimonas, Sphingomonadales, and Acidobacteria were all positively correlated with cellular response to cation stress.
Disentangling the potential associations between the soil microbial taxa and functional traits via a network analysis. This analysis was based on the SpiecEasi method. Functional traits were selected as a functional subset (FOAM level 2) of the cellular response to stress, hydrocarbon degradation, transporters, and nitrogen cycle. Microbial taxa were the OTUs with a relative abundance of >0.05%. Yellow nodes represent functional traits; other nodes were colored according to phylum. The size of each node is proportional to relative abundance. The taxonomic distributions of these OTUs are also given. BTEX, benzene, toluene, ethylbenzene, and xylene.
Furthermore, we applied the random forest (RF) analysis to identify the main microbial predictors of the functional traits in FOAM level 1 (Fig. S7 to S10). Different taxonomic phyla contributed to the important variable for distinct functional traits. The Betaproteobacteria were the most important variable for predicting key functional traits, including utilization of sugar (P < 0.05), amino acid utilization biosynthesis metabolism (P < 0.01), nucleic acid metabolism (P < 0.05), hydrocarbon degradation (P < 0.01), carbohydrate-active enzymes (P < 0.01), transporters (P < 0.01), and cellular response to stress (P < 0.05). Other important variables for predicting key functional traits were the Parcubacteria for fermentation (P < 0.01); the Bacteroidetes for the superpathway of thiosulfate metabolism, the EMP pathway, and gluconeogenesis (P < 0.01); the Chloroflexi for homoacetogenesis, fatty acid oxidation, the nitrogen cycle, and hydrogen metabolism (P < 0.01); the Actinobacteria for the TCA cycle (P < 0.01); the Thermomicrobia for sulfur compound metabolism and the hydrolysis of polymers (P < 0.05); and the Armatimonadetes for saccharide and derivative synthesis (P < 0.01).
DISCUSSION
Elucidating the mechanisms of microbial succession and their resilience in response to environmental disturbances, including contamination from chemical pollutants, is one of the central issues in microbial ecology (8, 32). It is challenging to robustly characterize how microbial temporal dynamics are affected by contamination at the species and more deeply functional levels (30, 31). Taxonomic and functional profiling in the present study was used to evaluate the differences in the structure and functional potential of the soil microbial community in response to disturbances. Herein, our results revealed distinct patterns of resilience in microbial diversity and specific traits at both the taxonomic and functional levels in response to inorganic and organic soil pollutants. Additionally, plants imposed stronger selection processes on the taxonomic profiles than on the functional profiles of the soil microbial community.
Microbes are able to respond rapidly to local environmental change, and pollutants often have a great influence on microbial community structure (33, 34). In temporal microcosms, microbial diversity significantly decreased throughout soil enrichment subcultures treated with pollutants (7, 8, 35). In the present study, soil microbial alpha diversity showed a clear succession pattern in response to contamination, which decreased at first and then recovered. This suggests the importance of resilience of soil microbiomes for the stability of microbial ecosystems. The rapid increase in functional richness in the early stage of contamination revealed the high influx of new functional traits, which may be attributed to the rapid evolutionary adaptation of soil microbes to environmental disturbances (27). The significant decrease in dissimilarities in the microbial community between samples of the initial and later stages also suggested the resilience of soil microbiomes while contaminated. In particular, the dissimilarities between the day 0 and day 90 samples were significantly lower for the functional than for taxonomic profiles, and the former varied more rapidly than the latter. This result suggested that functional traits exhibited more resilience than taxonomic ones. Microbes have high growth rates and undergo rapid evolutionary adaptation via horizontal gene transfer; therefore, their functioning could have shifted while they themselves did not change under environmental disturbances (26, 27). The resilience of soil function is attributable to acclimatization, and functional microbial members capable of rapid adaptation could quickly become abundant under environmental disturbances (30). Previous studies have demonstrated that the reshaping and constraining of soil functional responses to altered environments were driven by microbial community composition and their influence on the rates of ecosystem processes (36, 37).
Although the microbial alpha diversity apparently recovered in contaminated soils, substantial shifts in taxonomic and functional structure still occurred; this indicates that the recovery of soil microbial alpha diversity to a stable state was largely influenced by the pollutants we applied. The resilience of a community may also be defined and viewed as the recovery process following a disturbance leading to an alternative stable state (28, 29). Most of the traits were classified into two clusters, both of which exhibited a trend of resiliency, with the lowest difference in abundance levels found between the day 0 and day 90 samples. The OTUs belonging to Massilia, Lysobacter, Pseudoduganella, and Bacillus increased first but then decreased, pointing to their potential abilities of hydrocarbon degradation and metal tolerance. Massilia spp. are reportedly phenanthrene degraders (38), and Bacillus spp. were dominant in the n-octadecane degradation consortia (7, 39). In contrast, many more OTUs (n = 884) were classified into cluster 2, which decreased first and then increased. This result may be explained if most of the microbes had suffered immediate toxicity of heavy metal at the early stage of contamination, yet they could later adapt to this stress condition via high growth rates and evolutionary acclimatization during the incubation process. Our results are consistent with those from a prior study demonstrating that those microbial taxa which responded positively to disturbance eventually decrease in abundance to return to their original composition, allowing other negatively impacted taxa to then recover in abundance (27). Additionally, we also found that the functions related to cellular response to stress increased first and then decreased, thus indicating the microbes’ rapid functional responses to the selection pressure from inorganic pollutants (e.g., cadmium). In particular, the functions associated with hydrocarbon degradation showed opposing trends, which may be explained by members with the ability of hydrocarbon degradation being suppressed by the heavy-metal toxicity. Overall, these observations confirm the high resilience in the structure and functions of soil microbiomes in response to environmental contamination at not just the whole-community level but also at the species-specific level.
The selective effect of plant roots on the soil microbial community structure has been reported in many other studies and could be attributable to plant species identity (16, 25). Disturbance from pollutants reduces the number of active microbial taxa to only those that are pollutant tolerant; therefore, microorganisms that are generally associated with a particular plant species in uncontaminated soils may no longer be a relevant component of the community in an aged contaminated environment (40). In the present study, we observed no significant difference in microbial community compositions between the planted and unplanted samples after 10 or 30 days of incubation, while significant differences in the microbial taxonomic and functional profiles were detected after 90 days of incubation. Our results indicated that the plant-induced differences in the microbial community only occur after several months, whereas the effects of chemical contamination occur very quickly. This confirms that plants exerted a selection process on both taxonomic and functional traits of the soil microbial community in a contaminated environment. This also implies that plants growing in contaminated agro-soils have a significant impact on how the soil microbial community reassembles and recovers over time. Yet, our results showed that the different legume plants drove the distinct assembly of species composition rather than functional traits. There are two possible reasons for this result. First, different microbes in the distinct community may function similarly, thus resulting in the same soil functional profiling. Second, the microbial taxa selected by different plants may actually be functionally redundant, so that the soil microbial functions as a whole would not change when combined and considered at the community level. Prior studies have found that plant roots can exert significantly selective effects on microbial community structure, but to a lesser extent, or not at all, on the functional profiles (25, 41).
In contaminated soil environments, the influence of plant traits on the microbial community assemblage at the taxonomic and functional levels remains unclear. Here, we found that the enriched microbial taxa in soils with different plants were selectively assembled, and this might in turn benefit plant growth and health. These potentially beneficial microbes may support nutrient acquisition for promoting plant growth and are typically referred to as plant growth-promoting bacteria. For example, Mesorhizobium spp. enriched in soils planted with Robinia spp. have been mainly isolated from the black locust nodule (42). Rhizobium was the main rhizobial genus to form nodules with hairy vetch plants (43), and it was abundant in the vetch-grown soils. Ensifer spp. could perform nitrogen fixation for legumes, solubilize inorganic phosphate, produce indole acetic acid and siderophores, and induce systemic resistance to collectively promote plant growth (44, 45). Novosphingobium spp. can produce the phytohormones salicylic acid, gibberellins, indole-3-acetic acid, and abscisic acid (46). That said, in our study, the enrichment processes in the planted soils selected microbial functional genes specifically related to saccharide and derivated synthesis, the EMP pathway, and the response to nitrosative and oxidative stress. These particular functions appear to be relevant for interactions with the plants, some of which have been shown to be important in root-associated microbial communities (25, 41). For example, functions participating in the response to nitrosative and oxidative stress were crucial for the adaptation of plants to abiotic and biotic stresses (47). In particular, the presence of both organic pollutants and heavy metal (CdCl2) could function as environmental stresses for the selection of these bacteria with potential abilities of hydrocarbon degradation and metal tolerance. Indeed, Ensifer spp. have been shown to be capable of degrading phenanthrene (PHE) (48). In addition, a few genera enriched in soils planted with legumes, such as Novosphingobium, Nocardia, Sphingobium, Hydrogenophaga, and Sphingopyxis, were previously detected in microcosms where organic contaminants served as the sole sources of energy and carbon in the presence of CdCl2 (7, 8). Hence, our results regarding both species composition and functional traits suggest that plants may exert selection pressure on the soil microbial community based on its particular functional traits in the organic- and inorganic-contaminated environments. This selection process could benefit plant growth and fitness while fostering the potential abilities of hydrocarbon degradation and metal tolerance.
Network analysis based on SpiecEasi method provided an integrated understanding of the microbial community taxonomic and functional traits. Our analysis contained all possible positive and negative correlations between the taxonomic and functional traits, which could better reflect the many complex interaction relationships among species and potential functions of soil microbial community (25). In general, Proteobacteria had the most external connections, which indicated their adaptation across a wide range of ecological niches (49). The low interactions of Actinobacteria with other taxa may be due to their antibiotic activity (50). Furthermore, we estimated the associations of specific microbial taxa to certain environmental functional groups, including the cellular response to stress, hydrocarbon degradation, transporters, and the nitrogen cycle. Massilia spp. and Gaiellales were positively correlated with the major facilitator superfamily and ABC transporters, respectively, which encode important proteins with heavy-metal transportation ability (51, 52); this finding indicates the potential roles played by these two taxa in heavy-metal detoxification. Additionally, our results suggested that the bacteria belonging to Tumebacillus, Gemmatimonadaceae, and Acidobacteria may participate in hydrocarbon degradation, consistent with results from previous studies (53). Functions of cellular response to stress were associated with defending against threats from the environment (54). In our study, Gaiella spp., Actinoplanes spp., Brevundimonas spp., Sphingomonadales, and Acidobacteria all exhibited strong relationships with these functions, thus indicating their contributions to counteracting the stress from organic pollutants and heavy metal in agricultural soil.
Different taxonomic phyla likely explained the distinct functional traits of soil microbial communities. Betaproteobacteria, which are widely distributed across different ecological niches (49), are important predictors for many important functional traits, such as amino acid utilization biosynthesis metabolism, nucleic acid metabolism, and hydrocarbon degradation, highlighting their crucial roles in soil ecosystem processes. Parcubacteria (candidate phylum OD1) was found to predict fermentation functions; this is supported by other work inferring that Parcubacteria fermented various sugars to organic acids, with some species having the capacity to degrade complex carbon sources (55, 56). Bacteroidetes and Thermomicrobia contributed most to the functions related to sulfur cycling, and some of their members are known to be associated with sulfur metabolism (57). Overall, the combined network and RF analyses suggest that various taxonomic groups jointly contributed to the functional potential of the soil microbial community we studied, although specific phyla were the most important predictors for distinctive functional traits.
In this study, we observed clear succession and resilience patterns of soil microbial diversity and structure at the taxonomic and functional levels in response to soil contamination. In addition, different legume plants exerted stronger selection processes on taxonomic than functional profiles of soil microbiome in the organic- and inorganic-contaminated environments; these processes were potentially beneficial to the growth and health of plants with abilities of hydrocarbon degradation and metal tolerance. Importantly, the functional potential of the soil microbial community reflected contributions made by not one, but various cooccurring taxonomic groups, which did not preclude some phyla being the important predictors for specific functional traits. Overall, the results presented here provide valuable information for better understanding microbial resilience patterns to environmental contamination, not only on the community as a whole but also at the species and more deeply functional levels, as well as the selection processes acting on the taxonomic and functional profiles as induced by different plant hosts. Further work should investigate a range of land use types, environmental conditions, and soil properties to confirm the general and global applications of these findings. Ideally, such work could be extended to identifying the consequences of microbial temporal resilience on ecosystem functions and services.
MATERIALS AND METHODS
Experimental setup and sampling.In July 2014, 20 kg of soil samples (0 to 20 cm) was collected from a cornfield located in Yangling, Shaanxi Province, Northwest China (108°4′51″E, 34°17′31″N). This soil was of a sandy loam texture, with pH 8.16. Other soil properties were as follows: available nitrogen, 14.2 mg kg−1; total nitrogen, 1,013.9 mg kg−1; available potassium, 169.1 mg kg−1; available phosphorus, 18.2 mg kg−1; total phosphorus, 960.7 mg kg−1; and organic matter, 22.1 g kg−1. A subset of soils were stored at −80°C until the microbial analysis (referred to as day 0 of soil communities). The remaining soils were sieved (5-mm mesh size) to remove any plant debris and large clods. To prepare the contaminated soils, subsamples of the sieved soils were spiked with a mixture of phenanthrene and n-octadecane in dichloromethane at 1,000 mg/kg, plus cadmium chloride (CdCl2) in water at 50 mg/kg. The dichloromethane solution containing the organic pollutants was first mixed with 200 g of soil. After the complete evaporation of the dichloromethane under a fume hood, it was thoroughly mixed with a further 800 g of soil and CdCl2 solution.
In total, four time points for temporal microcosms without plants, 0, 10, 30, and 90 days, were selected to estimate the succession of soil microbiome in response to chemical contamination. Three legume plant treatments were sampled after 90 days of growth in the contaminated soils and used to determine the microbial community assemblage occurring under the influence of different legume plants. For all treatments, the contaminated soils (∼1 kg) were filled into pots (10-cm diameter) that had a depth of 10 cm. For the plant growth treatments, seeds of R. pseudoacacia (robinia), M. sativa (alfalfa), and V. villosa (vetch) were surface-sterilized and germinated at 28°C for 36 h under aseptic conditions. Five robinia plants, 20 alfalfa, and 20 vetch seedlings (each 1 cm in length) were sown per pot. Pots with and without plants were incubated in a greenhouse (16-h day [25°C]/8-h night [20°C]) for 90 days. Pots were given sterile water three times per week to maintain their soil moisture at ∼15%. Pots assigned to the different treatments were arranged randomly and rotated regularly throughout the 90-day period. Nine replicate pots were maintained for each of the four time points without plant (0, 10, 30, and 90 days) and for each of the three legume species at 90 days.
At the designated time points, the soils in pots without any plants were collected from a depth of 2 to 10 cm, which adopted destructive sampling to avoid disturbance. After the 90-day incubation, soils in the plant-grown pots were collected following the same procedure. Plant roots with soil attached were removed, and the rest of the soil without roots was mixed and collected. Although the rhizosphere consists of the soil most affected by plants, the sample volume collected here was too small. Therefore, we focused on the soil microbiome on a larger scale, that is, in the root zone impacted by plant growth. In total, 63 soil samples were obtained, and the three replications were pooled in groups, leaving 21 soil samples contaminated with phenanthrene, n-octadecane, and CdCl2, with three biological replicates for four time points (0, 10, 30, and 90 days), and three legume plant treatments for 90-day growth. All these soil samples were stored at −80°C for future use. At the end of the experiment, the remaining contaminated soils were collected and sent to an environmental protection company for processing via a chemical method.
DNA extraction and 16S rRNA gene amplicon and metagenomic sequencing.Genomic DNA was extracted from each of the 21 soil samples. The V4-V5 region of the 16S rRNA gene was amplified by using the primer pair 515F/907R. Sequencing was performed on a HiSeq 2500 (250-bp paired-end reads) platform (Illumina, Inc., San Diego, CA, USA). Metagenomic analysis of the 21 soil DNA samples was conducted with the sequencing libraries generated after implementing the NEBNext Ultra DNA library prep kit for Illumina (NEB, Ipswich, MA, USA). This sequencing was performed on a HiSeq 4000 platform (Illumina, Inc.). More detailed information is given in the supplemental material.
For the analysis of metabolism pathways, the KEGG database (58) was used, and amino acid alignment against the KEGG database was performed using BLASTP (E value ≤ 1 × 10−5) in DIAMOND (59). To annotate the set of genes related to environmental microorganisms, the identified gene families (specified by KEGG Orthology groups) were screened against the Functional Ontology Assignments for Metagenomes (FOAM) database and then grouped into different FOAM levels (60).
Statistical analyses.All statistical analyses were conducted using R version 3.2.2 (http://www.r-project.org), unless otherwise stated. For the taxonomic data sets, a subsample with a minimum of 42,994 sequences (according to the sample size) from each sample was used to eliminate all potential inaccuracies. The alpha diversity (richness and Shannon-Wiener index) and beta diversity (Bray-Curtis distance) were calculated for the taxonomic and functional traits of the soil microbial community based on their corresponding OTUs and open reading frames, respectively.
The dynamic patterns of the taxonomic and functional traits during the 90-day incubation period were identified using maSigPro (61) in the Bioconductor package, as described previously (8). A two-step regression approach was performed to select the traits with statistically significant stage changes (using false-discovery rate-corrected P values of <0.05). The visualization of significant traits was based on a cluster analysis for those group traits sharing a similar profile. The microbial taxa (i.e., OTUs) significantly enriched for each plant treatment relative to the soils lacking plants (day 90) were identified by linear statistics with the “limma” package. Specifically, a linear model was fitted for each OTU, and differentially abundant OTUs were estimated using moderated t tests. The moderated t test used a Bayesian model to shrink or expand the standard error of each OTU toward a common value, thereby borrowing information on the variance of other OTUs. This approach is especially suitable when the number of measurements per sample is large but the sample sizes are small. The resulting P values were adjusted for multiple-hypothesis testing by applying the Benjamini-Hochberg correction method. A bipartite association network was then used to visualize the associations between the significantly enriched OTUs and the different plant samples, which was generated in Cytoscape.
Network analyses were carried out to obtain a better understanding of the potential associations between the microbial taxa and functional traits of FOAM levels 1 and 2. Correlations among the OTUs with a relative abundance of more than 0.05% and the functional traits were calculated by sparse inverse covariance estimation for ecological association inference (SpiecEasi) using the Meinshausen and Bühlmann algorithm. Networks were visualized in the Cytoscape platform. The main microbial predictors for potential functions were identified by a classification RF analysis. In the RF models, the main microbial phyla served as predictors for all functional traits of FOAM level 1 (n = 21). The significance of the models and the cross-validated R2 values were assessed with 5,000 permutations of the response variable by using the A3 package. Similarly, the significance of the importance measures of each predictor on a given response variable (i.e., functional trait) was assessed by using the rfPermute package.
Data availability.The 16S rRNA gene amplicon sequencing and metagenomics data reported in this paper have been deposited in the Genome Sequence Archive in the Beijing Institute of Genomics (BIG) Data Center, Chinese Academy of Sciences, under accession number PRJCA001123 (http://bigd.big.ac.cn/bioproject/browse/PRJCA001123).
ACKNOWLEDGMENTS
This work was supported by the National Science Foundation of China (grants 31570493, 31270529, and 41807030), the Cheung Kong Scholars Program (grant T2014208), the National Postdoctoral Program for Innovative Talents of China (grant BX201700005), and the China Postdoctoral Science Foundation (grant 2018M630041).
S.J., W.C., and G.W. conceived and designed the experiments, and S.J. performed the experiments, analyzed the data, and wrote the paper.
We declare no conflicts of interest.
FOOTNOTES
- Received 15 October 2018.
- Accepted 3 January 2019.
- Accepted manuscript posted online 18 January 2019.
Supplemental material for this article may be found at https://doi.org/10.1128/AEM.02523-18.
- Copyright © 2019 American Society for Microbiology.
