Selection imposed by local environmental conditions drives differences in microbial community composition across geographically distinct groundwater aquifers.

Several studies have analyzed biogeographic distribution patterns of microbial communities across broad spatial scales. However, it is often unclear to what extent differences in community composition across different regions are caused by dispersal limitation or selection, and if selection is caused by local environmental conditions alone or additional broad-scale region-specific factors. This is especially true for groundwater environments, which have been understudied in this context relative to other non-subsurface habitats. Here, we analyzed microbial community composition based on exact 16S rRNA amplicon sequence variants (ASVs) from four geographically separated aquifers located in different regions along a latitudinal transect of ∼700 km across Germany. Using a combination of variation partitioning and ecological null models revealed that differences in microbial community composition were mainly the product of selection imposed by local environmental conditions and to a smaller but still significant extent dispersal limitation and drift across regions. Only ∼23% of the total variation in microbial community composition remained unexplained, possibly due to underestimated effects of dispersal limitation among local communities within regions and temporal drift. No evidence was found for selection due to region-specific factors independent of local environmental conditions. © FEMS 2019.

INTRODUCTION

The elucidation of processes that determine the biogeographic distribution of species is a central goal in community ecology and has further moved into the focus of microbial ecology over the past decade (Martiny et al. 2006; Nemergut et al. 2013; Meyer et al. 2018; Langenheder and Lindström 2019). Different theoretical concepts have been developed, which share overlapping perspectives on the processes that cause differences in community composition across space (Vellend 2010; Chase and Myers 2011; Leibold and Chase 2018). In essence, these processes can be broadly categorized as processes related to species sorting (selection), dispersal and random ecological and evolutionary drift due to stochastic migration, extinction and speciation events. These processes can act simultaneously and to varying degrees depending on the spatial scale at which communities are being investigated (Chase and Myers 2011). Within regions, sets of local communities that are linked by dispersal, and hence share a common regional species pool, have been defined as metacommunities (Leibold et al. 2004). Local communities may differ from each other as a result of species sorting caused by differences in local environmental conditions that select for distinct sets of species that are able to thrive under these conditions, provided that a certain degree of dispersal exists between communities to allow species to track environmental gradients. In such a case, differences in community composition would be strongly linked to differences in environmental conditions, while spatial distance between locations would have little effect (Leibold et al. 2004). However, community composition can be uncoupled from environmental conditions by processes that affect species dispersal. Apart from species sorting, differences between local communities can also arise due to dispersal limitation, in which case the impeded exchange of species in combination with random drift causes communities to diverge over time (Chase and Myers 2011). Additionally, dispersal limitation no longer allows species to track environmental differences and reach locations with their preferred environmental conditions, causing differences in community composition to be predominantly associated with spatial distance between locations rather than differences in environmental conditions. On the other hand, high similarities between communities can arise, and species sorting be over-ruled, under conditions with high-dispersal rates. Strong dispersal can homogenize local communities and allow species to occur even under unfavorable environmental conditions, if dispersal rates are sufficient to outpace species extinction caused by the unfavorable conditions. Also, in this case, differences in community composition would be less well predicted by differences in environmental conditions (Leibold et al. 2004).

Since local communities are subsets of a regional metacommunity, they are not independent of the regional species pool from which they are assembled. Understanding the processes that structure these metacommunities might, therefore, be an important step towards explaining local community diversity and biogeographic distribution patterns (Ricklefs 2008; Lindström and Langenheder 2012). Broad-scale region-specific factors like climate, geology, land use and historical events can leave a lasting imprint on community structure and thus may add to or even over-ride the effect of local contemporary environmental conditions (Martiny et al. 2006; Andersson et al. 2014; Fukami 2015; Stegen et al. 2016b; Vass and Langenheder 2017; Rummens, De Meester and Souffreau 2018; Svoboda et al. 2018). However, a common problem with studies that have investigated differences in community composition over large spatial scales is that distance decay relationships or the contribution of spatial distance in variation partitioning models can potentially be caused both by dispersal limitation over large spatial distances as well as by selection due to broad-scale regional factors (Leibold, Economo and Peres-Neto 2010; Hanson et al. 2012; Wang et al. 2013). Therefore, it is often unclear to what extent adaptation to these broad-scale factors, relative to dispersal limitation and drift, contribute to differences between communities across regions. Moreover, their effect on local community composition within regions, on top of dispersal, drift and selection imposed by the local environment, which may shape these communities as discussed above, is not well understood (Heino et al. 2017).

Several studies have investigated patterns of microbial biogeography over different spatial scales in various habitats, including soil, marine environments and surface freshwater systems like ponds, streams and lakes (for reviews see Hanson et al. 2012; Lindström and Langenheder 2012; Langenheder and Lindström 2019). However, the majority of studies so far have focused on local communities within a region, while only a few studies have explicitly investigated communities at larger spatial scales across regions (e.g. Martiny et al. 2011; Almasia et al. 2016; Comte et al. 2016; O'Brien et al. 2016; Ma et al. 2017; Hassell et al. 2018; Power et al. 2018; Shi et al. 2018). Moreover, a common conclusion from the reviews cited above is that the importance of selection relative to dispersal and drift not only changes depending on spatial scale, but also varies among habitat types (see also Wang et al. 2013). Strikingly, studies on groundwater ecosystems and related subsurface habitats are largely under-represented in this context compared to studies on surface habitats. Although recent studies have investigated the impacts of environmental conditions and selection processes on microbial community composition in subsurface environments (Shabarova et al. 2014; Stegen et al. 2015; Beaton et al. 2016; Graham et al. 2016a, 2017; Stegen et al. 2016a; Stegen et al. 2016b; Stegen et al. 2018; Savio et al. 2019), they mainly focused on local communities within a single aquifer, whereas studies that compared communities across aquifers from different regions are scarce (Ben Maamar et al. 2015; Danczak et al. 2018). Considering that groundwater-saturated zones of the terrestrial subsurface constitute the largest inland aquatic habitat for microorganisms on earth, and are estimated to harbor a significant fraction of the global microbial biomass (Griebler and Lueders 2009; McMahon and Parnell 2013; Magnabosco et al. 2018), the limited understanding of processes that shape microbial communities in these environments is a critical knowledge gap in microbial ecology.

In this study, we analyzed microbial community composition based on exact 16S rRNA amplicon sequence variants (ASVs; Callahan et al. 2016a) from four distinct shallow porous aquifers located in different catchment areas along a latitudinal transect of ∼700 km across Germany. We used a combined approach of variation partitioning and ecological null models to determine, first, the individual effects of local environmental conditions, spatial distance within regions and region identity on microbial community composition; second, the contributions of selection, dispersal and drift to community turnover within as well as across regions and third, the extent to which differences in community composition across regions were the result of dispersal limitation and selection caused by either broad-scale region-specific factors or local environmental conditions.

MATERIALS AND METHODS

Sample collection

A total of 45 samples were collected on single sampling campaigns between spring 2016 and summer 2018 from four distinct unconfined, shallow, porous aquifers, mainly consisting of unconsolidated gravel and sand, located in four different regions across Germany (Fig. 1). Region NOR (n = 12; September 2018; 53.72°N, 10.01°E) was located in Norderstedt near the city of Hamburg in the catchment of the Elbe River; region WUR (n = 13; May 2016; 49.77°N, 9.93°E) was located in Würzburg in the Main River catchment; region AUG (n = 12; June 2016; 48.25°N, 10.90°E) was located near Augsburg in the Lech River catchment; region MIT (n = 8; July 2018; 47.41°N, 11.26°E) was located near Mittenwald at the foothills of the German Alps in the Isar River catchment. Groundwater from the sampled areas of all aquifers was classified as non-contaminated and is used for drinking water production in their respective regions. Prevalent types of land use in each region were forests, grasslands and fallow agricultural land.

Figure 1.

Schematic map of Germany. Locations of the investigated regions are shown as black squares; rivers are shown in grey.

Groundwater samples were collected from fully screened monitoring wells using a submersible pump (MP1; Eijkelkamp Soil & Water, Giesbeek, the Netherlands) that was submerged at about half the depth of the water column in the well. Stagnant well water was purged prior to sample collection by pre-pumping to allow for approximately two volume exchanges and until physicochemical parameters (electrical conductivity, pH, temperature, dissolved oxygen concentration) had stabilized. All sample containers were sterilized prior to sampling and rinsed with sample water three times before sample collection. All samples were kept in the dark at 4°C for transport to the lab and until processing. Samples for dissolved organic carbon (DOC) measurements were collected in glass bottles that were baked at 450°C for 4 h or soaked in 10% w/v sodium persulfate overnight. DOC samples were passed through a 0.45 µm filter (Millex-HV; Merck-Millipore, Carrigtwohill, Ireland) rinsed once with sample water prior to sample collection and were acidified on-site to a final pH ≤ 2 with HCl. Samples for total prokaryotic cell counts were collected in Falcon tubes and fixed with 2.5% v/v glutardialdehyde (final concentration) immediately after sampling. Samples for DNA extraction (5 L) were collected in autoclaved glass bottles or plastic containers rinsed three times with 1 M HCl followed by three washing steps with 80% v/v ethanol (residual ethanol was allowed to evaporate overnight). Cells were collected on a 0.2 µm polycarbonate filter membrane (Merck-Millipore) within 48 h after sample collection and stored at −20°C until DNA extraction using the protocol by Pilloni et al. (2012).

Measurements of environmental variables

To estimate local environmental conditions, we measured 13 physicochemical parameters for each groundwater sample in addition to total prokaryotic cell counts. Electrical conductivity, pH, temperature and dissolved oxygen concentrations were measured on-site using field sensors (WTW, Weilheim, Germany). DOC concentrations were measured as non-purgeable organic carbon using high-temperature combustion (680°C) and infrared detection of CO2 on a TOC-V CPH Analyzer coupled to an ASI-V autosampler (Shimadzu, Kyoto, Japan). Inductively coupled plasma atomic emission spectrometry (ARCOS; Ametec-Spectro, Kleve, Germany) was used for the determination of cations (calcium (measured spectral element line: 183.801 nm), magnesium (279.079 nm), potassium (766.491 nm) and sodium (589.592 nm)) with radio frequency power set to 1400 W and argon as plasma gas at a flow rate of 15 L min−1. Samples were introduced by a peristaltic pump connected to a micromist nebulizer with a cyclon spray chamber. Anion concentrations (chloride, nitrate, orthophosphate, sulfate) were determined by ion chromatography (Dionex ICS-1500; pre-column: Dionex AG4; analytical column: Dionex AS4; Thermo Scientific, Idstein, Germany) with Na2CO3 (1.8 mM) + NaHCO3 (1.7 mM) as eluent at a flow rate of 1 mL min−1. Total prokaryotic cell counts were determined by flow cytometry (FC500 CYTOMICS; Beckman Coulter, Brea, CA, USA) with instrument settings as in Bayer et al. (2016). Cells were fluorescently stained in a 500 µL sample aliquot with SYBR Green I (Invitrogen, Darmstadt, Germany) at a ratio of 1:10 000 and incubated in the dark at 37°C for 13 min. 100 µL suspension of fluorescent beads (Trucount Tubes; BD Biosciences, San Jose, CA, USA) was added to each sample as internal standard for quantification. Measurements were done in biological and technical duplicates.

16S rRNA amplicon sequencing and data processing

DNA concentrations in raw extracts were determined using the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Paisley, UK). DNA extracts were diluted to 1 ng µL−1 with EB buffer (Qiagen, Hilden, Germany) and used as template (1 µL) for amplification of the V4 region of 16S rRNA genes using the primer pair 515FB (5′-GTGYCAGCMGCCGCGGTAA) (Parada, Needham and Fuhrman 2016) and 806RB (5′-GGACTACNVGGGTWTCTAAT) (Apprill et al. 2015) extended with Illumina adapters. Each reaction (25 µL) contained 12.5 µL NEBNext high-fidelity 2X PCR master mix (New England Biolabs, Ipswich, MA, USA), 3.75 µL 2% w/v BSA (Roche Diagnostics, Mannheim, Germany), 0.5 µL of each primer (10 µM) and 6.75 µL nuclease-free water. Initial denaturation was achieved at 98°C for 30 s followed by 25 amplification cycles (98°C, 10 s; 50°C, 30 s; 72°C, 30 s) and final elongation at 72°C for 5 min. Each sample was amplified in independent triplicate reactions; triplicates were pooled after amplification. Pooled amplicons were purified using magnetic beads (AMPure-XP; Beckman Coulter) at a bead:sample ratio of 0.8 and an incubation time of 5 min at room temperature. After washing (twice; 200 µL 80% v/v ethanol) and air-drying (10 min, room temperature), amplicons were eluted from the beads with 30 µL EB buffer. Amplicon size and concentration were determined by capillary gel electrophoresis (Fragment Analyzer; Agilent Technologies, Santa Clara, CA, USA) using the DNF-473 standard sensitivity NGS fragment analysis kit (Agilent Technologies). 10 ng amplicons were used as template for index PCR using Illumina Nextera XT Index Kit v2 primers (Illumina, San Diego, CA, USA) according to the manufacturer's specifications and with the same polymerase as above. After purification and electrophoresis as above, barcoded amplicons were pooled in equimolar concentrations (4 nM) and used for paired-end sequencing (2 × 300 bp) on an Illumina MiSeq platform.

Sequence data were processed in R (version 3.5.0) (R Core Team 2018) using DADA2 (version 1.10.1) (Callahan et al. 2016a) for quality filtering, merging of paired reads, inference of ASVs and chimera removal, according to the workflow by Callahan et al. (2016b) with slight modifications as described in the following. Truncation length during quality trimming was set to 280 bp and 200 bp for forward and reverse reads, respectively, after primer trimming. Negative controls were excluded from the model building step to infer error rates from the sequence data. ASVs were inferred across all samples using pseudo-pooling, which prevents discarding ASVs as singletons based on the occurrence in a single sample, if it is represented by at least two error-free reads in at least two samples in the full data set. ASVs found in negative controls or with sequence lengths < 261 bp were discarded. In addition, to facilitate downstream processing and reduce sparsity of the data, ASVs with an abundance < 0.001% across all samples were removed. Taxonomic assignment was done using the online implementation of IDTAXA (Murali, Bhargava and Wright 2018) by mapping ASV sequences against the SILVA SSU reference database (release 132) (Quast et al. 2013) with a 50% confidence threshold. ASVs that were classified as mitochondria or chloroplasts were discarded, as well as ASVs that could neither be classified as bacteria nor archaea. To infer phylogenetic relationships, ASV sequence alignments were obtained using the ‘DECIPHER’ package (version 2.10.1) (Wright 2015) and passed on to FastTree (Price, Dehal and Arkin 2009) for building a midpoint-rooted phylogenetic tree. The final ASV table contained 9153 ASVs; abundances were rarefied to 6281 reads per sample, which was the lowest number observed in a single sample. For a number of samples rarefaction curves did not reach saturation at this depth, which led to an underestimation of ASV richness by ∼30% on average compared to the total expected richness in the communities based on estimated asymptotes obtained from extrapolating rarefaction curves according to Chao et al. (2014) using the ‘iNEXT’ package (version 2.0.19) (Hsieh, Ma and Chao 2016). However, comparing the ASV richness after rarefaction to the estimated asymptotic richness by linear regression showed that this underestimation was uniform across samples, and furthermore that Shannon diversity was almost unaffected by the rarefaction. Additionally, differences in community composition based on all three beta diversity metrics relevant to this study (Bray–Curtis dissimilarity, β-mean nearest taxon distance and β-mean pairwise distance; see below) were well maintained in the rarefied data set compared to the original unrarefied data. Therefore, we can assume that rarefaction did not distort the overall structure of the data and still allowed meaningful comparisons between samples within our data set (Fig. S1, Supporting Information). Sequence data are publicly available at the NCBI sequence read archive (accession no. SRP191753).

Data analysis

All analyses were done in R. ASV richness and Faith's phylogenetic diversity (PD) were calculated using the ‘picante’ package (version 1.7) (Kembel et al. 2010). Differences in microbial community composition were analyzed based on β-mean nearest taxon distance (β-MNTD), which is the mean phylogenetic distance of species in one community to their closest relatives in another community and thus focuses on short phylogenetic distances, i.e. the tips of a phylogenetic tree, indicating turnover of lineages that have diverged relatively recently in evolutionary history. We additionally used β-mean pairwise distance (β-MPD), which is the overall mean phylogenetic distance between species in two communities and thus also captures deeper phylogenetic distances, indicating turnover of deeper branching phylogenetic lineages (Fine and Kembel 2011; Liu et al. 2017). β-MNTD and β-MPD were calculated with abundance weighting using the functions ‘comdistnt’ and ‘comdist’, respectively, of the ‘picante’ package. Differences in microbial community composition were illustrated by non-metric multidimensional scaling (NMDS) using the ‘metaMDS’ function of the ‘vegan’ package (version 2.5.3) (Oksanen et al. 2018). Environmental variables were standardized to z-scores for all analyses. Variables containing censored data, i.e. values below the detection limit (nitrate: < 0.131 mg L−1; orthophosphate: < 22.2 µg L−1), were handled according to Helsel (2011) using rank-transformation with tied ranks for values below the detection limit. Overall environmental differences were calculated as standardized Euclidean distances between samples considering all measured variables. Permutational analysis of multivariate dispersion (PERMDISP; Anderson 2006) was used for pairwise tests of differences between regions in microbial community beta diversity, environmental differences and spatial distance between sampling locations, respectively, using the ‘betadisper’ function (‘vegan’) with 10 000 permutations. Differences in ASV richness were assessed using Kruskal–Wallis non-parametric analysis of variance with Dunn's mean rank sum tests for pairwise comparisons and Holm correction for multiple testing.

We applied the null model approach by Stegen et al. (2012, 2013) to study the contributions of selection and dispersal processes on community turnover within as well as across regions. This approach is based on a two-step procedure: first, under the assumption that phylogenetic similarity between closely related taxa approximates ecological similarity, and that dispersal between communities exists at least to a minimal degree over evolutionary time scales to allow species sorting to act and outpace the evolution of distinct communities in situ (Stegen et al. 2013), the strength of species sorting is evaluated in the first step based on the β-nearest taxon index (β-NTI). β-NTI is the standardized effect size of β-MNTD, which indicates how much the observed difference between a pair of communities differs from a null distribution of β-MNTD calculated with randomized phylogenetic relationships for which species labels and abundances are repeatedly shuffled across the tips of the phylogenetic tree. β-NTI < −2 (β-NTI > +2) indicates that species in two communities are phylogenetically significantly more (less) closely related than expected by chance, suggesting selection of similar (different) species in both communities (referred to as homogeneous and variable selection, respectively, sensu Dini-Andreote et al. (2015) and Stegen et al. (2015). β-NTI < 2 indicates no significant deviation from the null distribution, suggesting that processes other than selection are responsible for the observed differences in community composition, i.e. dispersal, dispersal limitation and drift. In this case, the RCBray index is used in the second step to identify these alternative processes. Because phylogeny is assumed to be irrelevant to the chance of species being subject to dispersal, dispersal limitation, or random drift, RCBray does not consider phylogenetic relationships to calculate differences between communities, but only uses information on species occurrence and abundance. RCBray is a measure for the departure of the observed Bray-Curtis dissimilarity between two communities from a null distribution of dissimilarities between probabilistically assembled communities, which include species proportional to their respective abundances in the two compared communities and their occurrence frequencies in the data set, while maintaining local species richness and the number of individuals. RCBray < −0.95 (RCBray > +0.95) indicates that two communities share significantly more (less) species than expected by chance, which is interpreted as homogenizing dispersal (dispersal limitation and drift) being responsible for the observed differences between a pair of communities. |RCBray| < 0.95 indicates that two communities share as many species as expected by chance, indicating drift acting alone. β-NTI and RCBray were calculated as in Stegen et al. (2012, 2013) with 999 randomizations. For the analyses within regions, β-NTI and RCBray were calculated for each region separately based on null distributions that only considered ASVs found within a given region. For the analysis across regions, β-NTI and RCBray were calculated across all samples with ASVs found in the full data set.

As mentioned above, the ecological inference drawn from β-NTI regarding the influence of selection on differences between communities is based on the assumption that phylogenetic similarity between species across short phylogenetic distances approximates ecological similarity. This requires that phylogenetic distance between species correlates positively with differences in environmental optima (i.e. environmental optima have a phylogenetic signal). We tested this assumption for our data set using Mantel correlograms as done by others (Wang et al. 2013; Dini-Andreote et al. 2015). Differences in environmental optima between ASVs were estimated as standardized Euclidean distances between relative abundance-weighted means for environmental variables that were shown to have a significant effect on microbial community composition by distance-based redundancy analysis (db-RDA; see below) (Stegen et al. 2012; Dini-Andreote et al. 2015). The phylogenetic signal was evaluated at phylogenetic distance steps of 2% of the maximum phylogenetic distance with Mantel correlograms using Pearson correlation and 999 permutations for significance testing; P-values were adjusted for multiple testing using progressive Holm correction (‘mantel.correlog’; ‘vegan’). The analysis was done for each region separately only considering ASVs found within a given region, as well as with ASVs found across regions in the full data set. For the latter, we randomly selected 4500 ASVs similar to Dini-Andreote et al. (2015) since an analysis comprising all 9153 ASVs was computationally unfeasible. In all cases significant positive correlations were found mainly over short phylogenetic distances (12%–18% of the maximum phylogenetic distance) confirming that the assumption of a phylogenetic signal was met for our data set (Fig. S2, Supporting Information).

We used variation partitioning based on db-RDA to examine the relative importance of local environmental conditions, spatial distance with regions, and region identity on differences in microbial community composition and the degree to which these components were responsible for species sorting (Legendre and Anderson 1999; Legendre 2007). To study the effect on differences in community composition, abundance weighted β-MNTD and β-MPD were used as response matrices in the db-RDA models, respectively. Local environmental conditions were represented by standardized environmental variables. To reduce variance inflation, collinear environmental variables (electrical conductivity, sodium, calcium, magnesium, chloride and sulfate concentrations) were replaced by the first principal component resulting from a PCA of these six variables (referred to as ionPC1). IonPC1 was significantly positively correlated with all six variables and explained 74% of the variance. Environmental variables were selected by forward selection using the adjusted R2 of a full db-RDA model containing all environmental variables as stopping criterion (Blanchet, Legendre and Borcard 2008) (function ‘ordiR2step’ with 10 000 permutations; ‘vegan’). Calculation of variance inflation factors (VIF) (function ‘vif.cca’; ‘vegan’) confirmed low degrees of redundancy among the selected variables in all models (all VIF < 2). The marginal significance of each selected environmental variable was assessed using permutation tests (function ‘anova.cca’ with 10 000 permutations; ‘vegan’). Spatial distance within regions and region identity were included as independent components in the db-RDA models following the approach used by Declerck et al. (2011). Region identity representing spatial distance across regions was included as a dummy-coded variable matrix. Spatial distance between sampling locations within regions were represented by a staggered matrix of Moran's eigenvector maps (MEMs) as described by Declerck et al. (2011), where MEM vectors were arranged in blocks such that each block represents the spatial variation between locations within a given region, while locations from different regions are assigned a value of 0. MEMs are sets of orthogonal vectors derived from principal coordinate analysis on Euclidean geographic distances between connected sites, where individual vectors represent distances between sites at different spatial scales (Dray, Legendre and Peres-Neto 2006). The MEM matrix was constructed using the ‘create.dbMEM.model’ function in the ‘adespatial’ package (version 0.3–2) (Dray et al. 2018). Permutation tests were used as above to assess the overall significance of each component (i.e. region identity, MEM matrix, and the set of selected environmental variables) in individual db-RDA models and only significant components were subsequently used for variation partitioning (function ‘varpart’; ‘vegan’).

To study the effect on selection processes, we repeated the db-RDA including the same three explanatory components as above, but this time using the β-NTI matrix calculated across all samples as response matrix. The rationale behind this approach is that changes in β-NTI should only result from selection since the effects of dispersal, dispersal limitation and drift are accounted for in the null distribution by maintaining species abundances within samples during the randomization of phylogenetic relationships (Stegen et al. 2013; Wang et al. 2013). Accordingly, the fraction of variation in β-NTI explained by variables used to estimate local environmental conditions indicates that these variables impose selection, whereas a significant effect of spatial distance or region identity would indicate selection by spatially structured unmeasured environmental variables or broad-scale region-specific factors, respectively, rather than dispersal limitation. Since db-RDA requires only positive distance values, β-NTI was scaled to range between 0 and 1 as in Stegen et al. (2013).

RESULTS

Differences in microbial community composition and environmental conditions within and across regions

Analyses of differences in environmental conditions and microbial community composition by PCA and NMDS based on β-MNTD, respectively, revealed distinct clustering of samples by region with little overlap of samples from different regions (Fig. 2). In terms of environmental conditions, regions along the North–South transect were mainly separated along the second PCA axis, mostly influenced by differences in pH and concentrations of oxygen, potassium and DOC. Samples from the WUR region additionally separated from the other three regions along the first PCA axis, mostly influenced by temperature, electrical conductivity and concentrations of various ions (summary statistics of individual environmental variables are listed in Table S1, Supporting Information). In contrast to the turnover of closely related ASVs across regions indicated by β-MNTD, clustering of communities was weaker when differences in community composition were assessed across broader phylogenetic scales based on β-MPD (Fig. S3, Supporting Information). Congruently, while > 65% of the ASVs were exclusively detected within a single region, the majority of higher taxonomic groups from phylum to genus level (∼69%–73%) were observed across more than one region, further corroborating that differences between regions were mainly caused by turnover of related ASVs within broader clades such as genus or family. Regardless of the taxonomic level, taxa that occurred in more than one region also showed higher average relative abundances suggesting that local communities were dominated by more widespread taxa (Fig. S4, Supporting Information). The most dominant taxonomic groups in all four regions were Alpha, Delta and Gammaproteobacteria, in addition to Bacteroidia, Actinobacteria and taxonomically unclassified bacteria (Fig. S5, Supporting Information). Despite the dominance of these classes, community evenness calculated at the ASV level was high in all regions with Pielou's index values ranging between 0.8 and 0.9 on average (Fig. S6, Supporting Information). Accordingly, average relative abundances of the most dominant individual ASVs within a single region were relatively low ranging between 0.6% and 6%. These dominant ASVs were predominantly found within the families Burkholderiaceae, Caulobacteraceae, Pseudomonadaceae and Rhodocyclaceae in the WUR, AUG and MIT regions, and Thiovulaceae, Gallionellaceae as well as members of the Thaumarchaeota in the NOR region (Table S2, Supporting Information).

Figure 2.

(A), PCA showing differences in local environmental conditions (O2: dissolved oxygen; DOC: dissolved organic carbon; K+: potassium; TCC: total prokaryotic cell counts; PO43−: orthophosphate; Na+: sodium; Temp: temperature; Cl−: chloride; SO42−: sulfate; EC: electrical conductivity; Ca2+: calcium; NO3−: nitrate; Mg2+: magnesium). (B), NMDS showing differences in community composition based on abundance-weighted β-MNTD.

Comparing the dispersion of samples in the PCA and NMDS analyses shown in Fig. 2 suggested that regions with larger differences in environmental conditions also displayed larger differences in microbial community composition. This was confirmed by significant differences of within-region environmental heterogeneity and differences in community composition, respectively, revealed by pairwise PERMDISP tests (Fig. 3). Also in this case, patterns of observed differences in community composition based on β-MNTD matched the patterns of differences in environmental conditions better than β-MPD (Fig. S7, Supporting Information).

Figure 3.

(A), Spatial distance between sites within regions. (B), ASV richness within regions (total number of ASVs observed within each region is given in parentheses). (C), Differences in local environmental conditions (standardized Euclidean distance considering all environmental variables) within regions. (D), Differences in microbial community composition (abundance-weighted β-MNTD) within regions. Asterisks indicate significant differences inferred from PERMDISP tests (10 000 permutations) (A, C, D) and Dunn's rank sum tests (B) (*P < 0.05; **P < 0.01; ***P < 0.001). Note: we chose to display distances on their original scales as distances to group centroids obtained from PERMDISP revealed the same patterns as shown in A, C and D.

Since we had to rely on access to pre-installed monitoring wells during the sampling campaigns, it was unfortunately not possible to obtain samples from each region with the same spatial coverage. However, these differences in spatial coverage did not seem to have biased the estimates of ASV richness (Faith's PD showed the same pattern as richness, Fig. S8, Supporting Information), differences in microbial community composition, or environmental differences. For example, even though region WUR had the smallest spatial coverage, it displayed the second highest alpha and beta diversity estimates as well as the second largest environmental differences (Fig. 3).

Effect of selection and dispersal processes on community turnover inferred from null models

When evaluated within the individual regions as well as for pairwise comparisons of communities across regions, median β-NTI values were not significantly different from the null expectation, except for the NOR region, thus indicating no significant effect of selection on community assembly on average. However, the distributions of β-NTI in all regions as well as for comparisons across regions were strongly positively skewed (Fig. 4). Calculations of the fractions of pairwise community comparisons indicative of the different turnover processes showed that the contribution of selection to the observed differences between communities varied for each region between 32% (WUR) and 75% (NOR) (Fig. 4). In most cases, variable selection was the dominating process, indicating that communities were more different than expected by chance, except for the AUG region, where homogenous selection was the dominating selection process, suggesting that communities were more similar than expected. The fractions not accounted for by selection processes were largely dominated by dispersal limitation and drift in most cases, or drift acting alone in the MIT region. Across regions, variable selection was the dominating process, accounting for 69% of the observed differences between communities while the remaining fraction was indicated to result from dispersal limitation and drift.

Figure 4.

(A), Distribution of β-NTI values for pairwise community comparisons. Dashed lines indicate the range of β-NTI under the null expectation of no significant effects of selection (< |2|). Solid lines within violins represent quartiles (1st, median, 3rd). (B), Contribution of individual turnover processes to observed differences in microbial community composition derived from null models as classified by Dini-Andreote et al. (2015) and Stegen et al. (2015). Null models were run for each region separately; for the analysis across regions, null models were run on the full data set, and only results for pairs of communities from different regions are shown.

Variation partitioning of differences in microbial community composition and changes in selection

We applied db-RDA and variation partitioning to identify environmental variables that shaped microbial community composition (β-MNTD), and to dissect the individual contributions of these variables relative to spatial distance within regions and region identity (Table 1). Contradictory to the null model results for the individual regions that hinted at dispersal limitation, spatial distance between sites within regions represented by MEMs did not have a significant effect on differences in microbial community composition in an individual db-RDA model (adjusted R2 = −0.03; P = 0.916) and were therefore not considered for variation partitioning. In contrast, region identity and variables representing local environmental conditions (pH, ionPC1, and concentrations of dissolved oxygen, orthophosphate and DOC) together explained ∼77% of the variation in community composition, of which the majority (i.e. ∼41%) was shared between both components. The effect of environmental variables alone was still significant after controlling for region identity (pH, dissolved oxygen) and explained ∼27% of the variation, whereas region identity alone explained only ∼9% after controlling for the effect of environmental variables. In contrast to the results obtained for β-MNTD, < 7% of the total variation could be explained for β-MPD, and the individual fractions explained by local environmental conditions and region identity were almost equally low (∼2%) (Table S3, Supporting Information). Hence, together with the results described above, differences in community composition both in response to local environmental and regional differences were best reflected by turnover across short phylogenetic distances represented by β-MNTD compared to turnover across broader phylogenetic scales captured by β-MPD.

Table 1.

Partition of variation in microbial community composition (abundance-weighted β-MNTD) and selection (β-NTI) between local environmental conditions (Env; significant environmental variables are listed in the rightmost column) and region identity (Reg). Env+Reg represents the total variation explained by both components; Env|Reg (Reg|Env) represents the marginal fraction of variation explained by each component after controlling for the other; Env∩Reg represents the fraction of explained variation shared between both components. The explained variation is given as adjusted R2. Significance of each component and individual variables was tested using 10 000 permutations (note: significance of Env∩Reg cannot be tested). Spatial distance between sites within regions represented by MEMs was not significant in either case (adj. R2 = 0, P > 0.9) and was therefore not included in the analyses.

Response matrix	Component	df	Adj. R2	P	Significant variables (P <0.05)
β-MNTD	Env	5	0.6772	0.0001	pH, O2, ionPC1*, PO43−, DOC
	Reg	3	0.4972	0.0001	Dummy-coded region identity
	Env+Reg	8	0.7691	0.0001
	Env∩Reg	0	0.4053
	Env|Reg	5	0.2719	0.0001	pH, O2
	Reg|Env	3	0.0919	0.0001
	Residuals	36	0.2309
β-NTI	Env	5	0.6618	0.0001	pH, O2, PO43−, NO3−, DOC
	Reg	3	0.3747	0.0001	Dummy-coded region identity
	Env+Reg	8	0.6238	0.0001
	Env∩Reg	0	0.4127
	Env|Reg	5	0.2492	0.0001	pH, O2, DOC
	Reg|Env	3	−0.0380	0.0022
	Residuals	36	0.3762

Principal component representing 74% of the variance in electrical conductivity and concentrations of sodium, calcium, magnesium, chloride and sulfate (all positively correlated with ionPC1).

Given that region identity alone explained a significant amount of the variation in community composition, we further explored to which extent this variation was due to dispersal limitation or caused by selection either by local environmental conditions or region-specific factors. To this end, we used variation partitioning as above, only this time using the β-NTI matrix calculated for the full data set as response matrix in the db-RDA (Table 1). As in the analysis above, spatial distance between sites within regions, which would reflect the contribution of spatially structured unmeasured environmental variables, did not have a significant effect (adjusted R2 = −0.04; P = 0.996), whereas region identity together with local environmental conditions (pH and concentrations of dissolved oxygen, orthophosphate, nitrate and DOC) explained ∼62% of the variation in β-NTI. However, the effect of region identity was strongly tied to the effect of local environmental conditions, such that the variation explained by region identity alone dropped to zero after controlling for the effect of environmental variables (note that negative adjusted R2 although significant is interpreted as zero (Legendre 2007)). In contrast, environmental variables alone were still significant (pH, dissolved oxygen, DOC) and explained almost 25% of the variation in β-NTI after controlling for region identity. About 38% of the variation was unexplained, representing regionally and spatially unstructured, unmeasured environmental variables that imposed selection.

DISCUSSION

The aim of our study was to establish the relative contributions of processes that cause variation in microbial community composition in groundwater environments across distinct aquifers located in different regions. We hypothesized that variation in community composition can be due to species sorting imposed by local environmental conditions measured at the time of sampling, and potential broad-scale region-specific factors like climate, geology or historical events, in addition to processes related to dispersal and drift within as well as across regions. Our analyses showed that differences in local environmental conditions were well reflected by differences in microbial community composition within regions. This observation points towards the influence of species sorting, where stronger environmental gradients within a region are predicted to increase niche diversity in a metacommunity, and hence cause different species to sort into local communities along these environmental gradients (Langenheder and Lindström 2019).

The results obtained from the null models only partially agreed with this observation. On the one hand, the different degrees to which selection was indicated to be responsible for the differences in community composition in the NOR region compared to the AUG and MIT region did match the observed differences in environmental heterogeneity for these regions. This would support the hypothesis outlined above that stronger environmental gradients increase the influence of species sorting. On the other hand, contradictory results were found for the WUR region, which showed the second largest environmental differences, but exhibited the lowest contribution of selection. However, we have to mention that parts of the aquifer in the WUR region are artificially recharged with treated river water during the summer months (i.e. May to October) but not during the rest of the year. The samples for this study were collected at the early stage about 2 weeks after the start of the annual infiltration period, which may have constituted a perturbation to the microbial communities. It has been shown that random colonization through dispersal and drift can gain importance on community assembly in disturbed environments (Ferrenberg et al. 2013; Zhou et al. 2014; Fukami 2015; Langenheder and Lindström 2019), which could explain the relatively low contribution of species sorting in the WUR region. Furthermore, the null models indicated relatively strong contributions of dispersal limitation acting alongside drift, especially in the WUR and AUG region. Although comparable results have been obtained in previous studies on microbial community assembly in groundwater environments (Stegen et al. 2013; Beaton et al. 2016; Graham et al. 2017), this is at odds with our observation that spatial distance within regions did not have a significant effect on differences in community composition in the db-RDA.

Such apparently conflicting results between distance-based regression approaches and ecological null models have previously been reported by Langenheder et al. (2017) in a study on community assembly in lake biofilms. There are two possible explanations for these observations. One is that the inferences drawn from the null models might be an oversimplification of the actual ecological processes that shape microbial communities. The null model approach assumes that selection should mainly manifest itself in stronger or weaker phylogenetic community turnover than expected by chance (i.e. significant values for β-NTI). The basic underlying assumption that phylogenetic relatedness tends to approximate ecological similarity between microbial taxa has been confirmed by previous studies (Stegen et al. 2012; Wang et al. 2013; Dini-Andreote et al. 2015; Martiny et al. 2015; Liu et al. 2017), and was further suggested by a significant phylogenetic signal of environmental differences between ASVs in our data set (Fig. S2, Supporting Information). In the light of these findings, inferring the effect of selection from phylogenetic community turnover seems valid. However, it is also known that certain microbial traits are phylogenetically not well-conserved (Martiny et al. 2015), as we will further discuss below, and therefore selection involving such traits would not be reflected by phylogenetic turnover metrics like β-NTI, but could still result in higher than expected taxonomic community turnover reflected by RCBray, which does not consider phylogenetic relationships. Thus, a significant deviation from the null expectation of RCBray could still be the result of selection processes even if phylogenetic community turnover does not deviate from the null expectation of β-NTI (Langenheder et al. 2017). Alternatively, it is possible that mere spatial distance does not appropriately reflect actual groundwater flow paths via which microorganisms may disperse in porous aquifers (Freimann et al. 2015; Schmidt, Cuthbert and Schwientek. 2017; Smith et al. 2018). In this case, differences in community composition would not necessarily correlate with spatial distance even if dispersal was limited between local communities. For our study, this seems to be the more likely explanation, as we did not find significant correlations between changes in RCBray and differences in environmental conditions within the individual regions (based on Mantel correlation tests with 10 000 permutations; all P > 0.05; data not shown). This suggests that we may have underestimated the effect of dispersal limitation between local communities in the db-RDA, although we cannot fully rule out that selection involving phylogenetically non-conserved traits may have played a role as well.

In this context, we further have to point out that the ability to detect effects of species sorting and dispersal in our study was limited to differences in community composition that could be resolved based on 16S rRNA sequences. It has been shown that even closely related strains of the same species with near identical 16S rRNA genes can differ significantly in their ecological preferences and show distinct biogeographic distributions (Hahn et al. 2016; Larkin and Martiny 2017; Chase et al. 2018; Choudoir and Buckley 2018). Therefore, our results need to be interpreted with the necessary caution, bearing in mind that the high degree of conservation of 16S rRNA genes did not allow complete differentiation between microbial ecotypes beyond the level of ASVs in our study.

Even though we might have missed the variation in community composition caused by dispersal limitation between local communities, the majority of the total variation in community composition evaluated across all regions (> 75%) could still be explained by local environmental conditions and region identity. Variation partitioning of β-MNTD revealed a larger marginal effect of local environmental conditions (27%) compared to the effect of region identity (9%). This strongly indicates that microbial communities were shaped by these local environmental conditions, whereas region-specific factors and dispersal limitation between regions only played a secondary albeit still significant role. This was furthermore supported by the large contribution of selection to differences in community composition across regions inferred from the null models, similar to findings reported by Danczak et al. (2018). It is worth noting that the four different regions in our study were sampled at different time points due to logistic constraints, although all sampling campaigns were conducted roughly in the same season, i.e. late spring and summer. Still, the variation in community composition possibly caused by temporal drift independent of environmental conditions, in addition to possibly undetected effects of dispersal limitation within regions, which both would be represented by the residual fraction of unexplained variation, was comparatively small (∼23%).

Interestingly, the effect of environmental conditions and region identity were mainly reflected by turnover of closest relatives between communities, i.e. turnover over short phylogenetic distances measured as β-MNTD, but not by turnover of deeper branching phylogenetic lineages measured as β-MPD. Responses of microbial communities to environmental conditions have previously been shown to affect turnover across short phylogenetic distances, which indicate evolutionary relatively recent adaptations (Wang et al. 2013; Liu et al. 2017). Deep phylogenetic distances on the other hand capture more distant evolutionary events (Fine and Kembel 2011), which we hypothesized may include region-specific adaptations or evolutionary origins of phylogenetic clades within regions (Ricklefs 2006). This, however, was not the case as > 93% of the variation in β-MPD could not be explained by region identity and local environmental conditions, which both had equally miniscule individual effects.

The dominance of species sorting by local environmental conditions is in line with the meta-analysis by Hanson et al. (2012), who compared studies on microbial communities across various habitats and spatial scales and found that environmental conditions explained most of the variation in microbial community composition in the majority of analyzed studies. Similar conclusions were drawn from a literature review by Lindström and Langenheder (2012). Additional evidence for the importance of local environmental conditions on microbial community structure in groundwater environments in particular was provided by Ben Maamar et al. (2015) who reported similarities in microbial community composition in relation to similar environmental conditions across three unconnected fractured groundwater aquifers, as well as by other studies on single aquifers within a region (Stegen et al. 2013; Beaton et al. 2016; Graham et al. 2017). Nevertheless, region identity still explained a significant fraction of the variation in community composition after controlling for local environmental conditions in our study, comparable to previous studies that compared microbial community composition over broad spatial scales in various aquatic and terrestrial habitats (Souffreau et al. 2015; Power et al. 2018; Plassart et al. 2019), or similar examples from studies on larger organisms (Declerck et al. 2011; Viana et al. 2016; Heino et al. 2017). However, in these studies it largely remained unclear whether such large-scale distance decay relationships were the result of dispersal limitation across regions or selection by regionally structured factors. Using the standardized effect size of differences in community composition obtained from null models like β-NTI in addition to raw metrics like β-MNTD allows making such a distinction, because β-NTI quantifies the degree to which the phylogenetic turnover between two communities is stronger (or weaker) than expected given the observed differences in species richness, occupancy and abundance caused by dispersal and drift (Stegen et al. 2013; Wang et al. 2013). By partitioning the variation in β-NTI between region identity and local environmental conditions, we could show that local environmental conditions, both explained by measured variables and by unmeasured, spatially unstructured variables represented by the residual fraction, explained most of the variation in selection, whereas region identity alone did not have a significant effect. Combined with the results obtained for β-MNTD, this leads to the conclusion that the variation in β-MNTD explained by region identity was mainly due to dispersal limitation and drift across regions rather than species sorting imposed by broad-scale regional factors.

CONCLUSION

Our study showed that differences in microbial community composition across distinct aquifers from different geographic regions were mainly the product of species sorting imposed by local environmental conditions, with a relatively smaller but still significant contribution of dispersal limitation and drift across regions. However, we did not find evidence for significant selection effects caused by region-specific factors independent of local environmental conditions (represented by both measured and unmeasured variables). Although species sorting also played a determining role in structuring local microbial communities within the individual regions, we found partially inconsistent results between distance-based analyses and ecological null models regarding the contribution of dispersal limitation and drift within regions. Hence, combining microbial community analyses with hydrological models to map groundwater flow paths and identify possible dispersal routes for microorganisms will be important for future research to allow for more accurate estimates of the contribution of dispersal to microbial community assembly in groundwater environments. This in turn would be an important step towards a better understanding of the link between microbial community composition and biogeochemical functions in these ecosystems (Graham et al. 2016b; Graham and Stegen 2017).

Click here for additional data file.