Global patterns in gene content of soil microbiomes emerge from microbial interactions

Abstract

Microbial metabolism sustains life on Earth. Sequencing surveys of communities in hosts, oceans, and soils have revealed ubiquitous patterns linking the microbes present, the genes they possess, and local environmental conditions. One prominent explanation for these patterns is environmental filtering: local conditions select strains with particular traits. However, filtering assumes ecological interactions do not influence patterns, despite the fact that interactions can and do play an important role in structuring communities. Here, we demonstrate the insufficiency of the environmental filtering hypothesis for explaining global patterns in topsoil microbiomes. Using denitrification as a model system, we find that the abundances of two characteristic genotypes trade-off with pH; nar gene abundances increase while nap abundances decrease with declining pH. Contradicting the filtering hypothesis, we show that strains possessing the Nar genotype are enriched in low pH conditions but fail to grow alone. Instead, the dominance of Nar genotypes at low pH arises from an ecological interaction with Nap genotypes that alleviates nitrite toxicity. Our study provides a roadmap for dissecting how global associations between environmental variables and gene abundances arise from environmentally modulated community interactions.

Figure 1:. pH is associated with covariation in denitrification pathway composition in the global topsoil microbiome.

(A) Topsoils sampled at n = 189 globally-distributed sites were chemically characterized (pH, Ca, Mg, …), sequenced via shotgun metagenomics, and functionally annotated by Bahram et al. (3). (B) The relative abundance of denitrification reductases in each soil sample (relative to total gene content) are plotted in order of increasing total relative abundance. Reductase color legend indicated in panel A. (C-D) Unit-invariant singular value decomposition (uiSVD) (30) was used to decompose the data in panel B into contributions due to pathway magnitude (d_i) and pathway composition $\left(\overrightarrow{c_i}\right)$. The results of this decomposition are plotted in panel D. (E-F) Pathway magnitude values and principal component (PC) scores for pathway composition obtained via uiSVD (Methods) are compared with 17 environmental variables, and squared Pearson correlation coefficients are shown. Pathway magnitude is most correlated with C/N ratio (ρ = 0.59, p < 10⁻⁶ via one-tailed randomization test; Fig. S1), while scores of PC2 are most correlated with pH (ρ = 0.64, p < 10⁻⁶ via one-tailed randomization test; G). (H) Loadings of PC2 are shown, where positive values indicate reductase content that increases with pH, and vice versa. See also Figs. S1, S2.

Figure 2:. Enrichment cultures recapitulate global patterns in denitrification gene content.

(A) Denitrifying communities were enriched from soil samples in two pH conditions (6.0 and 7.3), and then sequenced via shotgun metagenomics to measure taxonomic composition and genotypes. Six soil samples were mechanically homogenized and used to inoculate a serial dilution experiment in denitrifying (anaerobic) conditions. After 72 h growth cycles, cultures were repeatedly passaged (12×) into a defined medium containing 2mM nitrate via a 1/8 dilution factor (Methods). (B) Endpoint community compositions are shown at the level of taxonomic order for enrichments in pH 6.0 and 7.3. Compositions were determined by extracting 16S rRNA fragments (miTAGs) from shotgun metagenome data (39), and then taxonomically annotating miTAGs via RDP (40). Taxa present at a relative abundance less than 1% are omitted. (C) Median denitrification reductase genotypes are shown for metagenome-assembled genomes (MAGs) corresponding to the four most abundant taxa present in the composition data in panel B. MAGs were functionally annotated using RAST (41), and the median genotype was computed over MAGs obtained in different samples. (D) Strains were isolated from cryopreserved samples from the enrichment endpoint (Methods). The strain Pseudomonas sp. PD Nar⁺ represents the dominant taxa present across samples at pH 6.0, and Rhizobiales sp. RH Nap⁺ represents the dominant taxa across samples at pH 7.3. See also Figs. S3, S4.

Figure 3:. Environmental filtering does not explain outcome of acidic enrichments.

(A) Schematic of the experimental design. PD Nar⁺ (blue) and RH Nap⁺ (orange) were passaged repeatedly in monoculture under denitrifying (anaerobic) conditions. After each 72 h growth cycle, cultures were passaged (3×) into a defined medium containing 2mM nitrate using a 1/8 dilution factor. Biomass is measured via 600 nm absorbance at the end of each cycle, and nitrate and nitrite concentrations are measured throughout each cycle using a Griess assay (Methods). (B and C) Endpoint biomass is shown for each strain at both pH 6.0 (panel B) and 7.3 (panel C). Lines connect the average across replicates. RH Nap⁺ (orange) produces more biomass at both pH levels, while PD Nar⁺ (blue) appears to decay in abundance at pH 6.0 (despite PD Nar⁺ community enrichment in this condition, Fig. 2B). (D-G) Nitrate and nitrite concentration dynamics are shown at six timepoints throughout each cycle in each condition. PD Nar⁺ nitrate and nitrite reduction rates slow at pH 6.0 as the growth-dilution cycles progress, consistent with its reduction in biomass (Panels D & F, blue). PD Nar⁺ (blue) accumulates nitrite at each pH condition, whereas RH Nap⁺ (orange) does not (Panels F & G). Technical replicates (n = 4) are shown for each strain in each experimental condition, with lines connecting the averages of these replicates

Figure 4:. Nitrite toxicity impacts denitrification activity of isolates at low pH.

(A) PD Nar⁺ and RH Nap⁺ were inoculated in anaerobic monoculture at pH 6.0 with 1.75mM nitrate and varying nitrite levels (colorbar). Nitrate and nitrite concentration dynamics, shown for each strain, were measured via Griess assay (Methods). As the initial nitrite concentrations ([NO₂⁻]₀) rose, the growth of both strains was inhibited. PD Nar⁺ was unable to fully reduce nitrate when [NO₂⁻]₀ > 0.35mM (blue and green curves, top left panel) and was unable to fully reduce nitrite when [NO₂⁻]₀ > 0mM (all except dark purple curve, bottom left panel). Similarly, RH Nap⁺ was unable to fully reduce either nitrate or nitrite when [NO₂⁻]₀ > 0.35mM (blue and green curves, right panels). Mean and standard deviation of technical replicates (n = 3) are shown. (B) PD Nar⁺ and RH Nap⁺ were again grown in anaerobic monoculture at pH 6.0 with varying nitrite levels, but nitrate was not supplied to prevent growth (Methods). Density of colony forming units was measured via plating with replicates (n = 3) for each condition, with points indicating means across replicates. Error bars were calculated by weighting across dilution levels, as described in Methods. Lines are log-linear fits, and the (negative) slope of the line indicates mortality rate (Methods). For PD Nar⁺, inferred death rates were 0.008 ± 0.004, 0.023 ± 0.004, and 0.039 ± 0.004 h⁻¹ for [NO₂⁻]₀ = 0, 0.875, and 1.75mM, respectively. For RH Nap⁺, inferred death rates were 0.013 ± 0.004, −0.001 ± 0.001, and 0.007 ± 0.001 h ¹ for [NO₂⁻]₀ = 0, 0.875, and 1.75mM, respectively. Uncertainties indicate standard deviations of log-linear fit parameters. (C) Nitrate (top) and nitrite (bottom) metabolite dynamics for monocultures of PD Nar⁺, supplied with initial nitrate concentrations of 0.875mM (light blue) and 1.75mM (dark blue). A consumer resource model was fit to the low initial nitrate data (solid light blue lines) and used to predict high initial nitrate data (dashed dark blue lines). The prediction failed at intermediate time points (bottom panel), indicating that the accumulation of nitrite slowed the metabolic rate. Inset in bottom panel shows final optical density (OD). Endpoint biomass levels were greater in the 0.875mM nitrate condition than in the 1.75mM condition, suggestive of mortality induced by nitrite accumulation. Mean and standard deviation of technical replicates (n = 3) of the data are shown.

Figure 5:. Community alleviates nitrite toxicity under acidic conditions.

(A-B) Monoculture PD Nar⁺ growth metabolite dynamics (shown previously in Fig. 4C) compared with the those of a 1:1 PD Nar⁺ and RH Nap⁺ co-culture at pH 6.0. Panel A shows nitrate concentration and B shows nitrite concentration. Initial biomass (OD₆₀₀ = 0.01) and nitrate (1.75mM) was the same in each condition, but the co-culture (purple) exhibited significantly less nitrite accumulation despite only slightly decreased nitrate reduction rate, as well as significantly higher biomass production (purple bar, inset, panel B). Points and error bars indicate the means and standard deviations, respectively, across technical replicates (n = 3). (C) Schematic of proposed interaction between PD Nar⁺ and RH Nap⁺. PD Nar⁺ (blue cell) reduces nitrate quickly and accumulates nitrite (thick blue arrow). In monoculture, this nitrite accumulation is self-inhibitory (dashed red line). In co-culture, however, RH Nap⁺ (orange cell) quickly reduces nitrite (thick orange arrow), alleviating the inhibitory effects of nitrite on PD Nar⁺ (purple line). (D) Schematic illustrating the multi-cycle co-culture experiment shown in panels E and F. Mixtures of PD Nar⁺ and RH Nap⁺ were prepared across a range of ratios spanning 0.03:0.97 to 0.97:0.03 (distinguished by color), with total biomass held constant. These mixtures were then transferred to fresh media buffered at pH 6.0 or pH 7.3, with 2mM nitrate supplied. Cultures were grown under anaerobic conditions for 72 h and passaged 1:8 into fresh media for a total of four growth cycles. 16S amplicon sequencing was used to infer PD Nar⁺ relative abundance at the end of each cycle in pH 6.0 conditions, and at the end of four cycles in pH 7.3. (E) PD Nar⁺ relative abundance dynamics at pH 6.0 are shown. f_Nar+_,0 = 0.03,0.5,and 0.97 are highlighted by darker lines. After a transient decrease in the f_Nar+_,0 > 0.1 conditions, co-cultures approach a value of f_Nar+ ≈ 60–0.75. (F) PD Nar⁺ relative abundance dynamics at pH 7.3 are shown. The relative abundance of PD Nar approached zero after four cycles in all conditions, indicating that it was competitively excluded by RH Nap⁺ in this pH condition. All relative abundance values and uncertainties are calculated as described in Methods.

Figure 6:. Nar⁺ and Nap⁺ phenotypes are conserved across diverse taxa.

Additional Nar⁺ and Nap⁺ soil isolates were grown in anaerobic conditions at pH 6.0, and their denitrification phenotypes were measured. (A) Denitrification phenotypes were summarized by area under the nitrate (top) and nitrite (bottom) curves (AUC). Metabolite dynamics for PDM04 Nar⁺ (blue) and XNM01 Nap⁺ (orange) monocultures are shown, with area under the curve (AUC) for each metabolite illustrated by the colored shaded regions. Smaller AUC indicates faster metabolite reduction. (B) AUC for nitrate (top) and nitrite (bottom) for three Nar⁺ and four Nap⁺ strains isolated in Ref. . Nar⁺ strains exhibit faster nitrate reduction (blue bars, top panel) and slower nitrite reduction (blue bars, bottom panel), and Nap⁺ strains are the opposite. The mean and standard deviation of the AUC across technical replicates (n = 2) are shown.

Figure 7:. Proteobacteria dominate low-richness denitrifying populations in topsoils.

Soil incubation experiments were performed to assess the response of endogenous communities to nitrate spike-ins. (A) Ten soil samples spanning a range of pH values (5.0–7.1) were collected from sites within the Cook Agronomy Farm (Pullman, WA, USA). Samples were processed to obtain soil slurries, which were then incubated in triplicate under anaerobic conditions for 96 h with spike-ins of nitrate which yielded 2mM final concentrations (nitrate⁺) and without (control). DNA was extracted from frozen endpoint samples for 16S V3-V4 amplicon sequencing, and differential abundance analysis of amplicon sequence variants (ASVs) was performed to statistically determine which taxa were significantly enriched under addition of nitrate (red points, scatter plot). Details describing these procedures are given in Methods. (B) The number of ASVs (grouped by phylum) significantly enriched by nitrate spike-ins are shown for each soil (ordered by increasing pH). Pie charts show the fraction of ASVs classified by PICRUSt2 (50) as possessing the enzymes to perform denitrification (blue), DNRA (red), or neither (grey). Median ASV counts and genotype fractions measured across technical replicate incubations (N = 3) are shown. Across the range of pH values, Proteobacterial ASVs are most frequently enriched, with a substantial fraction of these ASVs likely possessing the genes to perform denitrification. (C) Bars show the percentage of enriched denitrifying ASVs out of total detected ASVs across the range of soils, and fractions are shown above the bars. Median values across technical replicate incubations (N = 3) are given. The small fraction of significantly enriched denitrifying ASVs suggests that denitrifying communities in topsoils are low-richness.