Environmental connectivity controls diversity in soil microbial communities

Soil resident microbes

We chose sand as a realistic source of a mixed microbial community (which is referred to as the sand community or SC). Because the sand community cannot be preserved as a whole by freezing, we collected fresh material for better consistency for each experiment from the same spot in St. Sulpice near Lake Geneva (GPS coordinates: 46.508032N, 6.544050 E) as described in Moreno et al.⁵¹. Sampled sand at different seasons thus likely carried slightly different starting communities and cell densities. The sand was sieved through 2 mm² pores to remove large particles. The sieved sand was stored at room temperature and used within 7 days for extraction of resident microbial cells.

Microbial cells were extracted from four aliquots of 200 g of sand. Each 200 g aliquot was transferred into a 1-l conical flask and submerged in 400 ml of 21 C minimal media salts (MMS) (containing, per litre: 1 g NH₄Cl, 3.49 g Na₂HPO₄·2H₂O, 2.77 g KH₂PO₄, pH 6.8)²¹. Flasks were incubated at 25 °C under rotary shaking at 120 rpm for 1 h. The sand was allowed to settle and the supernatant was decanted into a set of 50 ml Falcon tubes, which were centrifuged at 800 rpm with an A-4-81 rotor and a 5810R centrifuge (Eppendorf AG.) for 10 min to precipitate heavy soil particles. Supernatants were decanted into clean 50 ml Falcon tubes and centrifuged at 4000 rpm for 30 min to pellet cells. The supernatants were carefully discarded and the cell pellets were resuspended and pooled from the four aliquots (i.e., from the initial 800 g of sand) in one tube using 5 ml of MMS. The pooled liquid suspension was further sieved through a 40 µm Falcon cell strainer (Corning Inc.) in order to remove any particles and large eukaryotic cells that may obstruct flow cytometry analysis (see below). A small proportion of the sieved liquid suspension was used to quantify the numbers of recovered cells (see below); the remainder was used within 12 h for bead encapsulation or for mixed liquid suspended growth (see below). With this gentle method, we extracted approximately 3 × 10⁵ cells g⁻¹ of sand.

Flow cytometry cell counting

Cell numbers in extracted soil communities and in the mixed liquid suspended growth experiments were counted by flow cytometry. SC-suspensions were diluted 100 times in MMS and stained in 200 µl aliquots with 2 µl of diluted SYBR Green I solution (1:100 in DMSO; Molecular Probes) in the dark for 30 min at room temperature. In some experiments, cells were additionally stained with 2 µl propidium iodide solution (10 µg ml^–1, Molecular Probes). Aliquots of 20 µl were aspired at 14 µl min^–1 on a Novocyte flow cytometer with absolute volumetric cell counting (ACEA Biosciences, USA). Cells were thresholded above a forward scatter signal (FSC-H) of 20 and further gated for propidium iodide-staining (excited at 535 nm and its fluorescence was collected at 617 ± 30 nm) and for SYBR Green I (excitation 488 nm, 530 ± 30 nm band-pass filter; channel voltage at 441 V) above values of 1000 (Supplementary Fig. 5).

Cell samples from the mixed liquid suspension growth experiments were diluted to approximately 10⁶ ml⁻¹ and subsampled to aliquots of 100 µl. The subsamples were then mixed with 100 µl of 8 g l^–1 sodium azide in phosphate buffered saline and incubated for 1 h at 4 °C to arrest cell respiration and growth. Samples were then stained with SYBR Green I as above and quantified by flow cytometry using the same thresholds and gates as describe above.

Bacterial strains and pre-culturing procedures

P. veronii 1YdBTEX2 is a toluene, benzene, m-xylene and p-xylene degrading bacterium isolated from contaminated soil²⁰. The strain was tagged with a single-copy chromosomally inserted mini-Tn7 transposon carrying a P_tac–mCherry cassette (Pve, strain 3433) as described in the ref. ⁵². A single P. veronii colony from a selective plate with toluene as the sole carbon substrate after 48 h incubation at 30 °C was inoculated into 10 ml of liquid MMS containing 5 mM sodium succinate as the sole carbon source and grown for 24 h at 30 °C with rotary shaking at 180 rpm²¹. After 24 h, the cells were harvested and washed for bead encapsulation or for comparative liquid mixed suspension growth, as described below.

Agarose bead encapsulation

SC cell suspensions containing between 2 × 10⁷ to 10⁸ cells ml^–1 were encapsulated in agarose using rapid mixing with pluronic acid in dimethylpolysiloxane and subsequent cooling, followed by sieving to achieve beads with a diameter range of 40–70 µm⁵³. The entire procedure was carried at room temperature and near a gas flame to maintain antiseptic conditions. 1% (w/v) low melting agarose (GEPAGA04-62, Eurobio ingen, France) was prepared in PBS solution (PBS contains per L H₂O: 8 g NaCl, 0.2 g KCl, 1.44 g Na₂HPO₄, 0.24 g KH₂PO₄, pH 7.4) and dissolved by heating in a microwave. The molten agarose solution was cooled down and equilibrated in a 37 °C water bath. Separately, 15 ml of dimethylpolysiloxane (Sigma-Aldrich, DMPS5X-500G) was poured in a 30 ml glass test tube. 1 ml of the 37 °C-agarose solution was mixed with 30 µl of pluronic acid (10% Pluronic® F-68, Gibco, Life Technologies) by vortexing at the highest speed (Vortex-Genie 2, Scientific Industries, Inc.) for a minute. Into this mixture of agarose and pluronic acid, 200 µl of prepared SC cell suspension at 0.2–1.0 × 10⁸ cells ml^–1 was pipetted and vortexed again at the highest speed for another minute. Five hundred microliter of this mixture was added drop-wise into the glass tube with dimethylpolysiloxane that was being vortexed at maximum speed. Vortexing was continued for 2 min. The tube was then immediately plunged into crushed ice and allowed to stand for a minimum of 10 min. After this, the total content of the tube was transferred into a 50 ml Falcon tube. The tube was centrifuged for 10 min at 2000 rpm using an A-4-81 swinging-bucket rotor (Eppendorf). The oil was carefully decanted while retaining the beads pellet. Fifteen milliliter of sterile PBS was added to the pellet and the beads were resuspended by vortexing at a speed set to 5. The tubes were again centrifuged at 2000 rpm for 10 min and any visible oil phase on the top was removed using a pipette. The process was repeated once more to remove any visible oil phase. Beads of diameter between 40 and 70 µm were then recovered by passing the PBS-resuspended bead content of the tube first over a 70-µm cell strainer (Corning Inc.). A further 5 ml of PBS was added to the cell strainer to flush remaining beads (<70-µm) into the filtrate. The collected bead filtrate was subsequently passed over a 40-µm cell strainer (Corning Inc.) to remove beads smaller than 40 µm. Recovered beads on the sieve were washed with an additional 5 ml of PBS, and any smaller beads in the filtrate that stuck to the bottom side of the cell strainer were gently removed by absorption with a Whatman 3M filter paper. After this, the sieve was inverted and placed on top of a clean 50 ml Falcon tube. 1.5 ml of incubation medium (MMS with the respective carbon substrates, see below) was used to collect the beads from the sieve into the tube. Two tubes were processed in parallel, which were pooled in the same final Falcon tube to yield a total volume of 3 ml. This volume was then split in three aliquots of 1 ml to make triplicate incubations. The encapsulation procedure produced ~1.2 × 10⁶ beads per ml, with an effective volume of 10% of the total volume of the liquid phase in the incubations. For comparative growth experiments, the procedure was repeated with the P. veronii pure culture. A visual guide to the procedure is presented in Supplementary Methods, Section 4.

Cell-in-bead growth incubations

Each of the incubation vials (with 1 ml bead solution) was further complemented with 4 ml MMS containing either mixed carbon substrates (“Mixed-C”, 0.1 mM C) or sand extract (see below) as growth substrates. Mixed-C solution was prepared by dissolving 16 individual compounds (Supplementary Table 2) in milliQ-water (Siemens Labostar) in equimolar concentration such that the total carbon concentration of the solution reached 10 mM C. These compounds are also listed in EcoPlates^TM (Biolog Hayward, CA, USA) and have been previously used as a soil representative substrates⁵⁴.

Sand extract was prepared by extraction with pre-warmed (70 °C) sterile milliQ-water. A quantity of 100 g sand was mixed with 200 ml milliQ-water in a 250 ml Erlenmeyer flask and swirled on a rotatory platform for 15 min, after which it was subjected to 10 min sonication in an ultrasonic bath (Telesonic AG, Switzerland). Sand particles were sedimented and the supernatant was decanted, and passed through a 0.22-µm vacuum filter unit (Corning Inc.). This formed the “sand extract”, of which 4 ml was added directly to the 1 ml bead suspension in the vials.

At the amended substrate concentrations, on average between 8 and 18% of SYTO-9-stained cells were detected on images outside (automatically detected) beads. As these proportions did not increase over time, we conclude that these were not cells growing outside beads but cells escaped from beads squeezed and broken between the microscope slide and coverslip during imaging. We therefore did not further take them into consideration when judging intrabead cell growth.

Mixed liquid suspension growth

For growth of SC-cells or of P. veronii in regular mixed liquid suspension, 2 ml of MMS with 0.1 mM mixed-C substrates were inoculated in quadruplicate with 5 × 10⁵ cells ml^–1 (prepared as described above for bead encapsulation and growth). As a negative control for background growth, MMS without added carbon substrate was inoculated. Assays were incubated at 25 °C with rotary shaking at 120 rpm, sampled daily (100 µl) for cell fixation, staining and flow cytometry counting of SC-cell numbers.

Bead sampling and microscopy

For sampling, the vials were removed from the incubator and beads were collected at 1200 rpm for 1 min using a swinging-bucket A-4-81 rotor (Eppendorf). An aliquot of 10 µl of bead suspension was carefully sampled from the bottom of the vials, mixed with 0.6 µl of 50 µM SYTO-9 solution to stain cells, and incubated for 20 min at room temperature in the dark. Vials were vortexed and placed back into the incubator. Five microliter of sterile milli-Q water was added to the stained beads and the complete aliquot (15 µl) was spread on a regular microscope glass slide to minimise aggregation of beads. A coverslip (24 × 50 mm) was gently placed to avoid air bubbles and excessive squeezing of the beads. Ten random positions on the slide were imaged with the 20× objective (NA 0.35) using an inverted AF6000 LX epifluorescence microscope system (Leica AG, Germany) equipped with a DFC350FXR2 camera. Every position was imaged in four sequential channels (phase contrast, 25 ms; mCherry, Y3-cube, 750 ms; SYTO-9, GFP-cube, 50 and 340 ms). The 50 ms-SYTO-9 channel exposure was used for analysis while the 340-ms exposure was used for verification of weak signals, if necessary. Images were recorded as 16-bit TIF-files and further processed using a custom MATLAB routine (described below).

Microscopy image analysis

A custom MATLAB image processing and analysis routine was developed to segment beads and microcolonies inside beads from the image-series⁵⁵. For each time-point and experimental replicate, the phase contrast, mCherry, and SYTO-9 images were read using the imread function built in MATLAB (version 2016b, MathWorks inc., USA). To identify the beads on each image, sharp changes in intensity were detected in the phase-contrast images using the edge function. Individual beads within a specific radius range were then identified using the imfindcircles function. In the next step, the microcolonies inside each bead were identified by thresholding and segmenting the mCherry and SYTO-9 images exclusively within the identified bead areas. mCherry and SYTO-9 images were further aligned to identify microcolonies in SYTO-9 having mCherry signal, which corresponds to Pve. Overlapping signals were considered to originate from a Pve colony if the area overlap between two channels was at least 30% or larger. If this were not the case, the areas were considered to consist of both Pve and SC microcolonies. Object areas smaller than 2 pixels were discarded. All microcolonies were thus differentiated as corresponding to Pve (mCherry plus SYTO-9 signal) or SC (SYTO-9 only), after which their area, fluorescence intensity and geometric distance (within the bead) were calculated.

Results were summarized for each incubation and time point to comprise the following information: (i) total number of beads containing microcolonies; (ii) the mean per bead productivity (PBP) (the product of the identified areas times their mean SYTO-9 fluorescence intensity, normalized over all analyzed beads); (iii) the normalized PBP distribution and its bin sum; and (iv) the PBPs for beads with single SC-cell occupancy or with SC-cell pairs.

Community diversity analysis

The diversity of the starting and growing SC communities was analyzed by high-throughput sequencing of amplified V3–V4 regions of the 16S rRNA gene. DNA was isolated from cell pellets from 1.5 ml liquid SC suspensions or from beads (2 ml) kept at –80 °C until analysis using a FastDNA Spin kit for soil according to the manufacturer’s instructions (MPBio). We targeted the V3–V4 hypervariable region of the 16S rRNA gene by amplification with the 341f/785r primer set and appropriate Illumina adapters and barcodes, following recommendations of the reagent supplier (https://support.illumina.com/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf). The 16S Amplicon PCR Forward Primer = 5′

TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG. The 16S Amplicon PCR Reverse Primer = 5′

GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC. Equal amounts of amplified DNA from each sample were pooled and sequenced bidirectionally on the Illumina MiSeq platform at the University of Lausanne Genome Technologies Facilities. Raw 16S rRNA gene amplicon sequences were quality filtered, concatenated, verified for absence of potential chimeras, dereplicated and mapped to known bacterial taxonomy (OTUs) using QIIME2 at 99% similarity to the SILVA taxonomic reference gene database on a UNIX platform⁵⁶. OTU-assigned reads were normalized to the sum of reads per sample and plotted by multi-dimensional scaling based on calculated Bray–Curtis distances as implemented in the R phyloseq package (v1.16.2).

Community growth modeling in high and low environmental connectivity

We developed a custom model (MATLAB) to test the effects of environmental connectivity and population growth on community diversity in presence or absence of interspecific interactions. Simulations were used to derive OTU diversity at stationary phase (OTU distributions, paired-growth, and diversity measures). Detailed description of the model is provided in Supplementary Methods, Sections 1–3. Code examples are available⁵⁵. Briefly, we modeled community growth from 200,000 in silico cells at the beginning (corresponding to the inoculated 2 × 10⁵ cells ml⁻¹ in the experiments), either in a high connectivity (as in our liquid experiments) or low connectivity (as in the agarose bead-encapsulated cells) evironment. In the latter case, we simulated both single (75% of all in silico beads) or pairs of founder cells (25%). Single cells were simulated as a one-dimensional vector with length of 200,000 and pairs as two one-dimensional vectors (2 × 200,000).

The initial OTU composition of the vector was derived from the experimentally determined OTU distribution from 16S rRNA amplicon sequencing in the washed soil suspensions before inoculation. The probability of occurrence of each OTU is computed as:

and was used in sampling the 200,000 initial cells at t = 0 h (Supplementary Methods, Section 1). To simulate the dynamics of growth, we employed Monod kinetics, described as:

where µ_sp,i and µ_max,sp,i denote the (maximum) specific growth rate of species i, (S) is the prevailing substrate concentration, and ({Ks}) is the concentration of the substrate at the point, where the specific growth rate is half of the maximum growth rate. All growth kinetic parameters except the growth rate were kept the same for all species. For simplicity and in absence of empirical kinetic data on SC cells, we assumed a general ({K}_{{rm{S}}}) = 0.3 × 10^–6 g ml^–1, constant ({{rm{yield}}},=,0.3) carbon to biomass conversion in g g^–1 and constant ({{rm{mass}}},=,120) fg cell^–1 for all species. The yield factor accounted for CO₂ losses from carbon metabolism. In contrast, we varied ({mu }_{{max }}) between 0.01 and 0.6 h^–1, and conditionally included growth penalty factors on single founder cells in the low connectivity environment (as explained below). We further examined the effect of having genotype-independent (randomly attributed) death at the beginning or death biased to fast growing genotypes (µ_max > 0.25). Finally, we tested further interaction terms that influenced attributed growth rates. The initial carbon concentration was set to 50 mg ml^–1, which allowed similar community development in terms of size (i.e., cell numbers) as in the experiments.

Growth was simulated for 120 time steps, corresponding to 60 h in the experiments, at which point the substrate is depleted and cells stop dividing (stationary phase). Based on the attributed growth rates to every OTU (i.e., every cell and genotype of the vector or double vector), the model calculates per time step how much substrate is converted into biomass (we allow continuous biomass formation) and lost in form of CO₂, which is subtracted to calculate the remaining substrate concentration for the next time step. When the overall substrate concentration is lower than ({S}_{{min }}) = 3 × 10^–6 g ml^–1, growth stops. The production of cell biomass is converted to cell numbers, which is then subsampled at the last time step per OTU (to an equivalent of 50,000 sequence reads), per single or pair (to an equivalent of 5000 beads) to calculate developed microcolony sizes, diversity measures, and interaction effects (as in Figs. 5 and 6).

The following interaction scenarios were simulated. Although we allowed growth penalties and interspecific interactions to influence attributed OTU growth rates, a threshold of ~0.6 h^–1 was imposed as the maximum individual OTU growth rate in all simulations.

High connectivity random vs. OTU-abundance growth rates. OTU-specific growth rates (µ_max,sp) were drawn randomly between 0.01 and 0.6 h^–1. Alternatively, growth rates were attributed to the vector of OTUs according to the probability distribution function reflecting their measured log₁₀ empiric abundance at t = 0 h (Supplementary methods, Section 2).

Low connectivity single founder cell growth penalty. We contrasted simulations with single founder cells growing according to their OTU-proportional attributed growth rate and those in which that growth rate was multiplied by a penalty, composed by a factor equal to the inverse proportion of the initially attributed µ_max,sp per OTU. The assumptiton is that the slower the inherent growth rate, the more likely that OTU is penalized when it is alone (Supplementary methods, Sections 2 and 3). This was combined with testing the effect of random or biased death on the starting community.

Low connectivity paired interspecific interaction effects. We further tested different assumptions on the nature of interspecific interactions and simulated how these affected community growth rates and diversity outcomes. These effects directly influenced the OTU attributed growth rates in the doubles (Supplementary methods, Sections 3.1–3.3). In the bimodal scenario (Supplementary methods, Section 3.4.1), we assumed that the community is composed of two underlying distributions; rare and abundant members (the threshold being placed at log₁₀ measured OTU relative abundance = 2.8), with abundant members having a higher probability to be positively influenced in pairs. The probability is drawn from a bimodal interaction curve that attributes an interaction factor (between 0.01 and 2.2), which is multiplied with the assigned OTU growth rate at start.

In the biased positive model (Supplementary methods, Section 3.4.2) we allowed a 40% chance for an interaction term imposed independently on each founder cell in a pair to lower the attributed OTU growth rate (factor range 0.4–0.6), and 60% chance for a factor in between 0.6 and 1.4 to modulate or increase the growth rate.

In the positive on slow model (Supplementary methods, Section 3.4.3), the attributed OTU growth rates on each founder cell in a pair had a chance of 40% to become improved inversely proportionally to its initial growth rate, thereby favoring slow growers

In the biased negative model (Supplemtary methods, Section 3.4.4), we attributed OTU-abundance proportional growth rates to each partner of the founder pair, but penalized faster growers (µ > 0.15) at 20% chance and the others at 40% chance that their growth rate would be multiplied by a negative interaction factor (range 0.01–0.1).

Finally, in a random model (Supplementary methods, Section 3.4.5), we allowed OTU-abundance proportional growth rates in pairs to be multiplied with a factor randomly drawn in the range of 0.01–1.25, independently for each partner in a pair. The models were contrasted to those without any assumed interspecific interactions, and without or with assumed random or fast-growing genotype biased cell death at start (Supplementary methods, Section 2.3.1).

All simulations were run five times from the beginning, independently producing five derived parameter values for alpha-diversity, OTU- and microcolony size distributions in stationary phase and partner interactions.

Statistics and reproducibility

Liquid suspension growth experiments were carried out in biological quadruplates and all bead experiments were carried out in biological triplicates. Total numbers of analyzed bead and those of beads with single or double occupancy are reported. Derived community growth rates and P. veronii-normalized yields were compared using t-tests (n as reported, two-sided test, unequal variance). Normalized PBP bin-size distributions were globally compared using Fisher’s exact test implemented in R (2000 replicates). Median and 75th percentile aggregate PBPs across different experiments were compared using the non-parametric Wilcoxon signed-rank test. Correlations between simulated species abundance distributions and empirical OTU relative abundances were calculated by bootstrapping (n = 1000) in MATLAB. Correlation coefficients from five independent simulations were compared using t-tests. The proportion correctly predicted OTU abundances by simulation was calculated as the ratio to observed values within a two-fold or four-fold range, and compared by two-sided t-tests on five independent simulations. Simulated and observed microcolony size distributions for single or paired founder cells among different models were compared by principal component analysis in MATLAB (pca), and by Spearman rank correlation (spear) from five independent simulations. Single and paired productivity was then compared between each other using two-sided t-tests of the 75th percentiles of microcolony size distribution (n = 5). Simulated and observed paired growth (excluding pairs with dead cells) was categorized and counted in a grid of 12 × 12 (each bin covering 0.5 log₁₀-distance) using MATLAB’s hist3d function, and then compared by pca from five independent simulations. Confidence intervals on ratios of paired simulated microcolony sizes (excluding those pairs with a non-growing or dead partner) were determined by subsampling (n = 1000) from mean ratio distributions, which were then used to calculate the fractions deviating from experimentally observed paired size ratios.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com