Ecological mechanisms underlying soil bacterial responses to rainfall along a steep natural precipitation gradient

Ecological mechanisms underlying soil bacterial responses to rainfall along a steep natural... Abstract Changes in the structure and function of soil microbial communities can drive substantial ecosystem feedbacks to altered precipitation. However, the ecological mechanisms underlying community responses to environmental change are not well understood. We used an 18-month soil reciprocal transplant experiment along a steep precipitation gradient to quantify how changes in rainfall affected bacterial community structure. We also conducted an enhanced dispersal treatment to ask whether higher immigration rates of taxa from the surrounding environment would accelerate community responses to climate change. Finally, we addressed how the composition of soil bacteria communities was related to the functional response of soil respiration to moisture in these treatments. Bacterial community structure (OTU abundance) and function (respiration rates) changed little in response to manipulation of either rainfall environment or dispersal rates. Although most bacteria were ecological generalists, a subset of specialist taxa, over 40% of which were Actinobacteria, tended to be more abundant in the rainfall environment that matched their original conditions. Bacteria community composition was an important predictor of the respiration response to moisture. Thus, the high compositional resistance of microbial communities dictated respiration responses to altered rainfall in this system. dispersal, drought, generalist, nestedness, reciprocal transplant, spatial turnover, specialist INTRODUCTION Soil microbes mediate enormous fluxes of carbon (C) among plant detritus, soil and the atmosphere via their roles in organic matter decomposition and soil respiration. Therefore, microbial responses to altered temperature, precipitation or nutrient inputs may drive substantial feedbacks to global climate change. We previously found that microbial respiration was a persistent function of historical climate rather than current conditions: soils from historically wetter sites exhibited higher respiration rates regardless of contemporary water availability (Hawkes et al.2017). Such functional legacies are likely created by microbial community resistance to environmental change (Hawkes and Keitt 2015), which has been observed in other long-term studies (e.g. Bond-Lamberty et al.2016). Nevertheless, high beta diversity in microbial communities across sites or through time is commonly reported (Freedman and Zak 2015; Kivlin and Hawkes 2016), consistent with microbial communities rapidly tracking environmental conditions. Here, we suggest that a better understanding of microbial community responses to the environment can be achieved by considering community structure beyond diversity or composition to include species nestedness, turnover and specialization. Most studies examining microbial responses to climate change focus on patterns in alpha and beta diversity. Soil bacterial communities exhibit extremely high species richness (alpha diversity) at the local scale and dramatic compositional turnover (beta diversity) among sites (Hughes et al.2001; van der Gucht et al. 2007). Few studies, however, have separately examined two independent components of microbial community beta diversity: spatial turnover and nestedness. Species turnover among sites is often equated with beta diversity, but nestedness is a component of beta diversity where the species pool at one local site is a subset of the species pool at a more diverse site (Baselga 2009). Partitioning the turnover and nestedness components of beta diversity is important because these factors respond to different ecological drivers. For example, community nestedness is highly sensitive to variation in local extinction rates and is often associated with environmental harshness (Ulrich, Almeida-Neto and Gotelli 2009). High species turnover, in contrast, is often driven by species sorting over environmental gradients (van der Gucht et al.2007). Whether spatial turnover or nestedness dominates community responses to climate change may depend on the interaction between two factors: habitat specialization within the microbial community and rates of dispersal. A microbial species can be considered a ‘habitat specialist’ if it is encountered more often in one type of habitat than any other (Mariadassou, Pichon and Ebert 2015). If the regional species pool is primarily composed of habitat specialists, then spatial turnover across environmental gradients should be high (Baselga 2009). By contrast, if habitat generalists predominate in the regional species pool, then nestedness may account for more of the variation in community composition among local sites, due to stochastic extinction of particular taxa at individual sites (Ulrich, Almeida-Neto and Gotelli 2009). In either scenario, high rates of dispersal can homogenize community structure across environmental gradients (Cadotte 2006). Dispersal may speed community acclimatization to altered environments if immigrants are better suited to contemporary conditions than residents (van der Gucht et al.2007; deVries and Shade 2013; Gianuca et al.2017). In this scenario, dispersal limitation would underlie microbial community resistance to change (Cline and Zak 2013; Lawrence, Bell and Barraclough 2016). Few studies have quantified the relative abundance of habitat specialists and generalists in soil microbial communities, although model simulations suggest that ecological specialization may drive community responses to environmental change (Hawkes and Keitt 2015). Evolutionary models demonstrate that large differences in the ‘optimal’ phenotype across different habitats tend to favor the evolution of ecological specialists, whereas generalists are expected in conditions that promote high phenotypic plasticity (van Tienderen 1997). High beta diversity of microbial assemblages suggests that most soil bacteria are specialists, yielding a strong positive correlation between taxon abundance and ecological specificity at broad spatial scales (Mariadassou, Pichon and Ebert 2015). Similarly, the degree of dispersal limitation within soil microbial communities is largely unquantified. Laboratory studies suggest that variation in dispersal rates can affect bacterial community structure and function (Lindström and Östman 2011), but there is little evidence for dispersal limitation in intact bacterial communities (van der Gucht et al.2007). Our goal in this study was to experimentally test microbial community responses to altered precipitation, focusing on how dispersal and habitat specialization shape these responses. To do so, we conducted a soil reciprocal transplant experiment across a 45-cm rain gradient in central Texas to expose soils from three dry and three wet sites to either to their normal ‘home’ or a changed ‘away’ precipitation regime (Fig. 1). In addition, half of the soil cores were subjected to elevated dispersal rates (i.e. inoculation) of microbes from the surrounding precipitation region. The experiment was designed to test four hypotheses: First, we expected to find differences in microbial community composition based on soil origin, driven by differences in historical climate and edaphic conditions at each unique site (Hypothesis 1). Next, we anticipated that the degree of microbial habitat specialization would structure communities both along the precipitation gradient and in response to transplantation treatments. If specialists are more abundant than generalists, spatial turnover should be high across the rainfall gradient; if generalists are more abundant, nestedness should play a larger role in generating beta diversity (Hypothesis 2). We also hypothesized that the enhanced dispersal treatment would reduce beta diversity (both turnover and nestedness) among soils of different origins (Hypothesis 3). Finally, we expected that changes in microbial community structure in response to transplantation and dispersal would drive the patterns of respiration in each community (Hypothesis 4) that were previously observed (Hawkes et al.2017). Figure 1. View largeDownload slide Schematic of the field reciprocal transplant experiment. Figure 1. View largeDownload slide Schematic of the field reciprocal transplant experiment. METHODS Study site and experimental design We worked in sites located across the Edwards Plateau in central Texas, where mean annual precipitation (MAP) decreases steeply from east to west (∼90 to 45 cm) while soil parent material and savanna grassland plant communities (dominated by Juniperus ashei, Bouteloua spp., Aristida spp. and Schizachyrium scoparium) remain relatively invariant. At these sites, mean rainfall peaks in early spring (April/May) and in the fall (September/October), but year-to-year variance is high. Mean annual temperature is approximately 20°C. Sample locations within sites were selected to minimize differences in topography (<2% slope) and the plant community (native grassland with <50% woody plant cover), and all soils are rocky, calcareous Mollisols. We conducted a reciprocal transplant experiment of soil cores from the drier west end to the wetter east end of the gradient and vice versa. To do so, we established two 6 m × 6 m common garden plots, one located at an eastern site where MAP = 89 cm (Ladybird Johnson Wildflower Center, or WFC, Lat: 30.18, Lon: –97.88) and one located 300 km to the west where MAP = 63 cm (MOR Ecolab, Lat: 29.44, Lon: –100.06) (Table S1, Supporting Information). Soils from three wetter eastern sites and three drier western sites were placed at each common garden in PVC cores (n = 6), creating a full factorial combination of home region and transplanted treatments. Moreover, six cores containing soils sampled on-site were also incubated in each garden as a control treatment (Fig. 1). Soils from all eight sites (three eastern sites, three western sites and the two common garden sites) were collected in July 2011, sieved to 2 mm to remove rocks and roots, packed into PVC cores at approximately native bulk density (150 g at ∼0.95 g cm−3) and installed in each common garden site within a week. Cores were arranged in rows, 50 cm apart, in random order with respect to treatment. The PVC cores measured 5 cm in diameter and 15 cm in depth, with 38-μm mesh bottoms to hold soil in while allowing for drainage. All cores were open at the top, allowing for passive dispersal of microbial propagules from the surrounding environment. The 38-μm pore size of the mesh core bottom meant that bacteria and fungi could also enter from surrounding soil. To examine the effects of dispersal rate on microbial community responses to the transplant treatment, half the cores at each garden (i.e. n = 3 in each treatment combination) received an ‘active dispersal’ treatment. We intended to simulate enhanced immigration rates of taxa from the surrounding precipitation region. Therefore, cores receiving the active dispersal treatment in the eastern common garden were inoculated with ‘eastern’ microbial assemblages, and cores in the western common garden were inoculated with ‘western’ assemblages. To create the inoculants, soils were freshly collected from the three sites in the surrounding precipitation region and the common garden site itself, thoroughly mixed in sterile water in equal proportions and filtered through 38-μm mesh to remove solid particles; 500 μL of the resulting microbial slurry was applied to each core in the enhanced dispersal treatment. A corresponding volume of sterile water was applied to the remaining experimental units. Enhanced dispersal cores were inoculated at establishment in June 2011, and again in November 2011, April 2012 and October 2012 to mimic naturally occurring variation in arrival rates (Jones et al.2017). With six experimental replicates for each factorial combination of common garden, soil origin and dispersal treatments, there were a total of 72 cores at each common garden site in addition to the six control cores (n = 78 per garden, 156 for the entire experiment). We measured a suite of physicochemical properties for each of the eight soil types included in our reciprocal transplant experiment. Soil C and N content were measured via combustion (Apollo TOC/TN, Teledyne Tekmar, Mason, OH, USA). Soil texture (% sand, silt and clay) was determined via the hydrometer method. Soil pH was measured in by mixing samples with distilled water in a 1:2 ratio. A summary of soil chemical data is given in Table S1, Supporting Information. Characterization of microbial community structure Cores were harvested in December 2012, after 18 months in the field, at which time they were removed from the ground intact and transported to the laboratory on ice. Soils from each core were sieved to 2-mm and homogenized prior to DNA extraction with a MoBio PowerSoil extraction kit (MoBio Laboratories, Carlsbad, CA). DNA extracts were amplified with bacteria-specific SSU primers Bact-0341 and Bact-0785 (Klindworth et al.2012), and the reverse primer was barcoded with an 11-bp treatment-specific barcode sequence. Amplicons were cleaned using the MoBio UltraClean PCR Clean-Up kit (MoBio Laboratories, Carlsbad, CA), quantified via Qubit dsDNA High Sensitivity assays (Life Technologies, Grand Island, NY) and pooled in equimolar mixtures. Pooled libraries were purified of fragments <150 bp in size using Agencourt AmPure XP magnetic beads (Beckman Coulter, Brea, CA). Libraries were sequenced on a 454 GS FLX sequencer with titanium chemistry (Roche, Branford, CT) at the University of Texas Genome Sequencing and Analysis Facility. We analyzed the resulting sequences using the QIIME pipeline v. 1.5.0 (Caporaso et al.2010). We eliminated sequences with quality scores below 25, fewer than 150 bases, homopolymer runs of more than 10 bp, with any ambiguous bases, or with an anomalous barcode. We also removed all singletons (unique to the entire dataset) and eliminated chimeras with USEARCH 6.1 (Edgar 2010). We used the UCLUST algorithm (Edgar 2010) to define operational taxonomic units (OTUs) at 97% sequence similarity. Representative sequences from each OTU were searched against the SILVA (Quast et al.2012) and RDP (Wang et al.2007) databases to identify and eliminate any non-bacterial taxa. To avoid bias in sampling effort, we rarefied to 1000 sequences, which was the smallest number of sequences in a sample. Sequences reads are archived in the NCBI Sequence Read Archive, Accession PRJNA368746. Quantifying microbial functional responses to soil moisture To determine whether the transplant or dispersal treatments affected microbial functional responses to water availability, soil respiration rate in the field common gardens was measured under ambient moisture conditions in November 2011, April 2012, October 2012 and December 2012. We capped cores with septa-embedded lids and sampled the headspace after 30 min. Headspace samples were stored in gas-tight vials until CO2 was quantified on a SRI 8610C gas chromatograph with FID and methanizer (SRI Instruments, Santa Monica, CA). We also examined soil respiration in response to controlled moisture conditions in the lab using data from Hawkes et al. (2017). Briefly, soils from the common garden cores were subsampled at harvest and assigned to high (25%) and low (10%) gravimetric soil moisture treatments in a laboratory incubation experiment, with microbial function (CO2 production) measured over 8 weeks. Headspace samples were analyzed for CO2 content as described above for field samples. Composition statistical analysis To examine microbial community structure and diversity across all experimental cores (Hypothesis 1), we used permutational ANOVAs to examine variation in bacterial community composition among treatments. Henceforth, we will use the term ‘region of origin’ to indicate the climatic province from which soils were originally sampled (wetter east vs. drier west), and ‘site of origin’ to indicate the specific site from which they were sampled (Table S1, Supporting Information). Therefore, in our models, site is nested within region of origin. Permutational ANOVAs were conducted using a combination of adonis in the vegan package (Oksanen et al.2016), with each term entered last in the model to obtain Type III sums of squares, and nested.npmanova in the BiodiversityR package (Kindt and Coe 2005) to obtain correct F-ratios and P values; 999 permutations were used per run. We also conducted planned pairwise comparisons among four treatment groups: soils that originated from the wetter eastern end of the gradient, incubating in home and away rainfall environments; and soils that originated from the drier western end of the gradient, incubating at home and away. A Bonferroni correction was applied to control for multiple comparisons. To visualize differences in community composition, we calculated non-metric multidimensional scaling (NMDS) coordinates based on Bray-Curtis dissimilarity matrices in PC-ORD (v6.2; McCune and Mefford 2011). Although statistics reported here are based on relative abundance data, all effect sizes and P values are nearly identical for presence-absence data. To examine variation in bacterial community composition in relation to soil properties, historical precipitation regime and recent drought stress, we used partial regularized discriminant analysis (RDA; Borcard, Legendre and Drapeau 1992; Legendre, Borcard and Peres-Neto 2005). Collinear variables were removed based on variance inflation factors >2.5, resulting in the following predictors in the final model: latitude, elevation, MAP, MAT, soil moisture in each core at harvest, the Palmer Drought Severity Index (PDSI) for 18 months in the home site prior to collection, PDSI for the experimental duration in the common garden sites, soil total organic carbon content (TOC), soil pH and soil % clay. Ecological specificity analysis We calculated the ecological specificity index (SI) following Mariadassou, Pichon and Ebert (2015) for the 1000 most abundant OTUs across the entire experiment. The index ranges between 0 and 1, with 0 indicating a perfect generalist (equal abundance across all sites) and 1 indicating a perfect specialist (found only in one of the eight sites). Specificity was calculated relative to the site of soil origin in which that OTU was most abundant; that is, we calculated local abundance in site X relative to average local abundance across all sites of soil origin, combined across gardens and dispersal treatments. Following these calculations, we performed a one-way ANOVA to determine whether the distribution of SIs varied among sites of soil origin. We also calculated a community-weighted mean SI for each site of soil origin, based on the presence or absence of each OTU within cores belonging to each site of soil origin. To determine whether the relative abundance of specialists or generalists shaped community responses to the transplant treatment (Hypothesis 2), we identified ecological specialists as those in the top decile of the SI distribution and ecological generalists as those in the bottom decile. Then, for each taxon, we calculated log response ratios (LRRs) for abundance in home vs. away transplant common gardens. We used a t test to ask whether the LRRs were greater for specialists vs. generalists, indicating higher abundance of specialists in their home precipitation region. Beta diversity partitioning We quantified community dissimilarity due to the turnover vs. nestedness components of beta diversity using the betapart package in R (Baselga and Orme 2012). We used the beta.pair function to generate dissimilarity matrices for the turnover and nestedness components of beta diversity across all experimental cores. Multi-response Permutation Procedure (MRPP) tests were conducted to determine whether dissimilarity associated with turnover and nestedness varied among soil sites of origin, transplant treatments and dispersal treatments (Hypothesis 3). To understand how turnover and nestedness were related to any shift in composition associated with transplantation, we also quantified these components of beta diversity across cores in home vs. away treatments within each site of soil origin. Finally, turnover and nestedness within each soil type were regressed against community mean SI indices to determine whether a higher relative abundance of ecological specialists is correlated with greater species turnover in response to the transplant treatment. Linking microbial community composition with functional responses to moisture To determine whether microbial community structure was a significant predictor of patterns of respiration in response to current moisture (Hypothesis 4), we used a regression approach to identify drivers of the respiration moisture response in the lab incubation experiment. Laboratory data were chosen because moisture conditions were more carefully controlled than in the field, but differences in CO2 fluxes between eastern- and western-origin soils were identical to those observed in situ (see ‘Results’). First we calculated ΔResp (microbial functional response to soil moisture) as the difference in mean CO2 flux in high vs. low soil moisture regimes in the laboratory incubation experiment. This calculation was performed for each of the 26 unique treatment groups (unique combination of soil site of origin × common garden location × dispersal treatment; and two common garden controls). Then we modeled ΔResp as a function of climatic, edaphic and microbial community variables. For our climatic predictors at each site of origin, we included the total rainfall experienced the year before transplant (June 2010 – 2011), henceforth referred to as ‘antecedent’ precipitation, and mean temperature at each site of origin during that time period (‘antecedent’ temperature). We excluded 30-year MAP and MAT because they were highly correlated with multiple predictor variables; see Table S3, Supporting Information. Because there was no significant main effect of common garden identity in analyses of bacteria diversity and composition, we did not include climate variables associated with the transplant period in the gardens. To explore the effect of soil physicochemical traits on respiratory responses, we used total organic C, pH and % clay in soils from each site of origin (total N was excluded because was highly correlated with soil C: R2 = 0.95). Finally, we included four predictor variables related to microbial community composition and structure: mean bacterial community richness and phylogenetic diversity in each treatment group, and mean NMS axis 1 and 2 scores in each treatment group. The latter were derived from the Bray-Curtis distance matrix of bacterial community composition. We used the regsubsets command in the R package leaps (Lumley 2017) to iterate through all possible combinations of predictors. The best models were selected by comparing BIC scores; pseudo R2 for each generalized linear model were calculated by comparing residual deviance vs. null deviance. RESULTS Diversity and community composition During the 18 months that soil cores were incubated in the two common gardens, cumulative rainfall was 54% higher in the eastern vs. western garden location (Fig. S1, Supporting Information), although both gardens experienced a 5-month severe drought at the start of the experiment. At the time of harvest, soil bacterial communities included 734 OTUs per sample on average (after rarefaction to 1000 sequences). However, alpha diversity was significantly affected by the interaction of site of soil origin and dispersal treatment (Table S2, Supporting Information): in soils originating from four sites, alpha diversity increased modestly by 2%–4% in the active dispersal treatment. All communities were generally dominated by Acidobacteria, Actinobacteria and Proteobacteria (Fig. S2a, Supporting Information). Bacterial community composition varied primarily by site of soil origin, explaining about 13% of the compositional variance despite 18 months of incubation in a new location (Table 1, Fig. 2). Site origin effects differed slightly between the two common gardens, but this explained only 2.6% of compositional variation (Table 1). The dispersal treatment had no effect on community composition, nor did region (drier west vs. wetter east) of soil origin. In planned pairwise comparison tests, eastern-origin and western-origin bacteria communities were always significantly different from each other (P < 0.01). However, there was no difference in community structure between eastern-origin soils in their home vs. away rainfall environment (F1,73 = 1.14, P = 0.294), and the same was true for western soils incubating at home vs. away (F1,73 = 1.14, P = 0.354). This indicates that bacterial community structure changed little in response to transplantation. Figure 2. View largeDownload slide NMS ordination of bacterial community composition in soils from each site, garden and dispersal treatment. (Stress = 22.44 after 56 iterations.) Only site and the interaction of site with common garden had significant effects. The dashed line separates soils originating from eastern and western regions. Average NMS scores for each site are plotted ±1SE (n = 4–6). Figure 2. View largeDownload slide NMS ordination of bacterial community composition in soils from each site, garden and dispersal treatment. (Stress = 22.44 after 56 iterations.) Only site and the interaction of site with common garden had significant effects. The dashed line separates soils originating from eastern and western regions. Average NMS scores for each site are plotted ±1SE (n = 4–6). Table 1. Results of permutational ANOVA examining bacteria community composition as a function of region of soil origin (Region), common garden location (Garden), dispersal treatment (Dispersal), soil moisture treatment (Moisture), site nested in region (Site|Region) and their interactions. Factors significant at P < 0.05 are in bold.   df  MS  F  R2  P  Region  1  1.516  1.446  0.031  0.163  Garden  1  0.457  1.426  0.009  0.918  Dispersal  1  0.322  1.029  0.007  0.993  Region*Garden  1  0.326  1.016  0.007  0.450  Region*Dispersal  1  0.300  0.959  0.006  0.557  Garden*Dispersal  1  0.307  0.993  0.006  0.998  Region*Garden*Dispersal  1  0.300  0.970  0.006  0.740  Site(Region)  6  1.048  3.446  0.128  0.001  Garden*Site(Region)  4  0.320  1.053  0.026  0.037  Dispersal*Site(Region)  4  0.313  1.030  0.026  0.183  Garden*Dispersal*Site(Region)  4  0.309  1.017  0.025  0.359  Error  117  0.305  0.724        df  MS  F  R2  P  Region  1  1.516  1.446  0.031  0.163  Garden  1  0.457  1.426  0.009  0.918  Dispersal  1  0.322  1.029  0.007  0.993  Region*Garden  1  0.326  1.016  0.007  0.450  Region*Dispersal  1  0.300  0.959  0.006  0.557  Garden*Dispersal  1  0.307  0.993  0.006  0.998  Region*Garden*Dispersal  1  0.300  0.970  0.006  0.740  Site(Region)  6  1.048  3.446  0.128  0.001  Garden*Site(Region)  4  0.320  1.053  0.026  0.037  Dispersal*Site(Region)  4  0.313  1.030  0.026  0.183  Garden*Dispersal*Site(Region)  4  0.309  1.017  0.025  0.359  Error  117  0.305  0.724      View Large In the RDA analysis, all environmental and spatial variables together explained only 10.7% of the variance in community composition (Table 2), despite wide variation in soil chemical properties (soil % C ranging from 2.4% to 10.9%; clay ranging from 18.5% to 39.7%; pH ranging from 7.4 to 8.2; Table S1, Supporting Information). Precipitation-related variables (MAP, gravimetric soil moisture and PDSI) accounted for 2.7%, soil characteristics (TOC, pH and clay) explained 2.1% and spatial effects explained only 0.1% of the variation in bacteria community composition among sites. Table 2. Results of partial RDA parsing how bacteria community composition relates to climatic (Clim), edaphic (Soil) and spatial (Space) factors. Most of the variation in the bacteria community was unexplained.   Inertia  R2  Adj R2  P  All vars  23.45  0.157  0.107  <0.001  Clim  20.19  0.054  0.027  <0.001  Soil  11.69  0.042  0.021  <0.001  Space  4.49  0.016  0.009  <0.001  Clim|Space + Soil  8.65    0.012  <0.001  Soil|Space + Clim  1.28    0.005  <0.001  Space|Clim + Soil  1.84    0.001  <0.001  Unexplained  125.73          Inertia  R2  Adj R2  P  All vars  23.45  0.157  0.107  <0.001  Clim  20.19  0.054  0.027  <0.001  Soil  11.69  0.042  0.021  <0.001  Space  4.49  0.016  0.009  <0.001  Clim|Space + Soil  8.65    0.012  <0.001  Soil|Space + Clim  1.28    0.005  <0.001  Space|Clim + Soil  1.84    0.001  <0.001  Unexplained  125.73        View Large Beta diversity partitioning Across all cores in the experiment, differences in bacteria community structure attributable to site of soil origin and transplantation treatments were generated by species turnover, not nestedness (Table 3). Neither turnover nor nestedness responded to the dispersal treatment because this treatment had no effect on community composition. Similarly, within sites of soil origin, changes in community structure in response to transplantation were dominated by species turnover (Table 4). The degree of nestedness (i.e. local extinction) across transplant treatments was generally quite small. Table 3. Dissimilarity among communities in different treatments due to species turnover and nestedness. P-values are based on MRPP tests on community dissimilarity matrices. Bold type indicates significance at P < 0.05.   Observed δ, (significance of δ)    Dissimilarity due to species turnover Expected δ 0.644  Dissimilarity due to nestedness Expected δ 0.021  Region of soil origin  0.628 (0.001)  0.021 (1.000)  Site of soil origin  0.570 (0.001)  0.022 (0.990)  Common garden site  0.642 (0.008)  0.021 (0.708)  Dispersal treatment  0.645 (0.562)  0.021 (0.634)    Observed δ, (significance of δ)    Dissimilarity due to species turnover Expected δ 0.644  Dissimilarity due to nestedness Expected δ 0.021  Region of soil origin  0.628 (0.001)  0.021 (1.000)  Site of soil origin  0.570 (0.001)  0.022 (0.990)  Common garden site  0.642 (0.008)  0.021 (0.708)  Dispersal treatment  0.645 (0.562)  0.021 (0.634)  View Large Table 4. Community-weighted mean specialization indices (SIs) for bacterial communities in each soil type; and the turnover and nestedness components of beta diversity between home and away precipitation treatments within each soil type.   Community-weighted mean SI  Turnover  Nestedness  Eastern sites  COL  0.208  0.123  0.010  ING  0.229  0.142  0.005  SCO  0.201  0.123  0.035  WFC (common garden site)  0.246  NA  NA  Western sites  DRS  0.218  0.136  0.040  SCS  0.218  0.141  0.005  SIE  0.213  0.150  0.005  MOR (common garden site)  0.247  NA  NA    Community-weighted mean SI  Turnover  Nestedness  Eastern sites  COL  0.208  0.123  0.010  ING  0.229  0.142  0.005  SCO  0.201  0.123  0.035  WFC (common garden site)  0.246  NA  NA  Western sites  DRS  0.218  0.136  0.040  SCS  0.218  0.141  0.005  SIE  0.213  0.150  0.005  MOR (common garden site)  0.247  NA  NA  View Large Ecological specialization Ecological specialization indices (SI) varied widely among taxa, ranging from 0.030 to 0.688, but the mean value of 0.184 suggests most taxa were not strong habitat specialists. However, a significant positive relationship between the specificity index and taxon rank abundances (R2 = 0.06, P < 0.001; Fig. 3) indicates numerical dominance by the most specialized taxa. There was significant variation in the degree of ecological specificity relative to site of soil origin (one-way ANOVA, F7,985 = 19.70, P < 0.001). In post-hoc tests, the highest SIs were found for taxa specialized to both home common garden sites (MOR, WFC; P < 0.001; Fig. 4a). Specificity indices were also significantly stronger for taxa specialized to soils originating from the drier western regions vs. the wetter eastern region (t = 2.051, P = 0.041, Fig. 4b). Community-weighted mean specialization indices in each site of soil origin followed a similar pattern, and were 18% higher in home common garden soils than in eastern-origin soils sampled from the wetter end of the rainfall gradient (Table 4). There was no relationship among community-weighted mean SIs and either the turnover or nestedness components of beta diversity. Figure 3. View largeDownload slide Relationship between OTU abundance and ecological specificity index (SI). Figure 3. View largeDownload slide Relationship between OTU abundance and ecological specificity index (SI). Figure 4. View largeDownload slide Indicators of bacteria specialization. Mean ecological specialization indices (+1 SE) for the most abundant bacterial taxa for (a) soils from each site of origin (arranged by MAP from left to right) and (b) in each precipitation region (wetter east vs. drier west). SI indices range from 0 (no specificity) to 1 (total specificity, i.e. the taxon is found only at one site). (c) Log response ratios (LRRs) for abundance in the home precipitation region vs. away region for specialist vs. generalist taxa. Positive LRRs indicate the taxon in question is more abundant in its home precipitation region. Figure 4. View largeDownload slide Indicators of bacteria specialization. Mean ecological specialization indices (+1 SE) for the most abundant bacterial taxa for (a) soils from each site of origin (arranged by MAP from left to right) and (b) in each precipitation region (wetter east vs. drier west). SI indices range from 0 (no specificity) to 1 (total specificity, i.e. the taxon is found only at one site). (c) Log response ratios (LRRs) for abundance in the home precipitation region vs. away region for specialist vs. generalist taxa. Positive LRRs indicate the taxon in question is more abundant in its home precipitation region. Specialist bacterial taxa (mean SI of 0.502) tended to be more abundant in home vs. away environments in comparison with generalist taxa (mean SI of 0.062; t = –1.795, P = 0.074; Fig. 4c). The taxonomic affiliation of specialists and generalists also varied significantly (χ2 = 25.729, P = 0.003): the largest fraction of generalist OTUs (28%) belonged to Proteobacteria, whereas 42% of specialist OTUs belonged to Actinobacteria (Fig. S2b, Supporting Information). Links between composition and function As previously reported (Hawkes et al.2017), rates of respiration were higher in soils that originated from sites located in the wetter eastern vs. drier western end of the gradient, in addition to current soil moisture conditions. This pattern was observed for in situ respiration rates measured in the field (Table S4a and b, Supporting Information) as well as in the laboratory microcosm experiment, where soil moisture was carefully controlled (Hawkes et al.2017). Our model selection procedures indicated that climatic, edaphic and microbial community variables could predict the respiratory response to soil moisture (ΔResp) moderately well: pseudo-R2 ranged from 0.18 to 0.33 (Table 5). As expected, the explanatory power of each model generally increased with the number of predictor variables. However, the two best (lowest BIC) models included only antecedent rainfall and NMS axis 1 scores as predictors. The magnitude of the respiratory response to soil moisture increased with the amount of rainfall soils received in the past year (Fig. 5a) and increased with NMS axis 1 scores (Fig. 5b), which were linked with the abundance of taxa in the Acidobacteria and Actinobacteria. NMS axis 1 scores, antecedent rainfall, soil pH and soil clay content were the variables most frequently included in the top-ranked models (Table 5); ΔResp decreased with soil clay content (Fig. 5c) and had a very weak negative relationship with soil pH (data not shown). Figure 5. View largeDownload slide The relationship between functional responses to soil moisture in the laboratory incubation (ΔResp) and climatic, edaphic and microbial community variables: (a) annual rainfall soils were exposed to prior to the transplant treatment; (b) bacteria community composition (as summarized by NMS 1 scores) and (c) soil clay content in each site of origin. Figure 5. View largeDownload slide The relationship between functional responses to soil moisture in the laboratory incubation (ΔResp) and climatic, edaphic and microbial community variables: (a) annual rainfall soils were exposed to prior to the transplant treatment; (b) bacteria community composition (as summarized by NMS 1 scores) and (c) soil clay content in each site of origin. Table 5. The five best-ranked models (as determined with BIC scores) predicting soil respiration responses to moisture in the laboratory incubation experiment. ‘Rain’ and ‘temp’ refer to antecedent rainfall and temperature at each site of origin in the 12 months prior to the transplant treatment; clay refers to soil % clay; NMS1 and NMS2 refer to mean scores on the first and second axis of an NMDS of bacteria community composition.   Predictors    Pseudo  Model rank  Biotic  Climatic  Edaphic  BIC  R2  1    Rain    – 109.46  0.182  2  NMS1  Rain    – 107.61  0.226  3  NMS1    Clay; pH  – 106.82  0.296  4  NMS1  Temp  Clay; pH  – 104.62  0.324  5  NMS1  Rain; temp  Clay; pH  – 101.58  0.329    Predictors    Pseudo  Model rank  Biotic  Climatic  Edaphic  BIC  R2  1    Rain    – 109.46  0.182  2  NMS1  Rain    – 107.61  0.226  3  NMS1    Clay; pH  – 106.82  0.296  4  NMS1  Temp  Clay; pH  – 104.62  0.324  5  NMS1  Rain; temp  Clay; pH  – 101.58  0.329  View Large DISCUSSION We found evidence for consistent differences in bacterial community composition among sites from across a natural precipitation gradient; this compositional variation was weakly structured by soil properties as well as historical and recent rainfall (Hypothesis 1). Although most taxa were habitat generalists, ecological specialists were locally abundant and tended to be more responsive to the transplantation treatments. The distribution of these specialist taxa likely contributed to the high spatial turnover observed across experimental cores (Hypothesis 2). Contrary to Hypothesis 3, however, the enhanced dispersal treatment had no effect on bacterial community structure or patterns of beta diversity. Finally, we found strong evidence that bacterial community composition mediated functional responses to soil moisture variation (Hypothesis 4). The characteristic respiratory patterns of soils from the wetter eastern and drier western regions of the rain gradient were maintained over the 18-month experiment, paralleling the compositional resistance of bacterial communities to environmental change. The most striking pattern we observed in our study was the high degree of resistance of bacteria community composition to novel rainfall regimes. This compositional resistance affected functional responses to water availability: eastern-origin (wetter) and western-origin (drier) soils exhibited characteristic patterns in soil respiration, which were unchanged following transplant and dispersal treatments. These functional patterns were correlated with variation in bacterial community structure among sites of soil origin. However, the relationship between climate, bacteria composition and respiration was not straightforward in our system, with MAP at each site explaining only ∼3% of the community compositional variance, but recent rainfall driving ∼18% of the functional response. There may be physiological adaptation of microbial taxa independent of composition, or indirect effects of rainfall on substrate availability for microbial growth (Moyano, Manzoni and Chenu 2013). Additionally, soil fungi, which we did not measure, likely also contributed to the observed patterns of respiration. Finally, the strong relationships between annual precipitation and microbial community composition that were observed along the gradient in 2008 (Hawkes et al.2017) may be lost over time, reduced by the disturbance required to construct cores or diminished under conditions of drought (we do not have initial microbial community data from the date of transplantation to address this issue definitively). Based on our results, the degree of environmental change relative to long-term climate may play a role in how resistant microbial communities are to perturbation, and therefore how community-level functions (like respiration) track changing environments. In our system, microbial communities at both ends of the gradient experience drought on a regular basis despite variation in MAP, potentially explaining why most bacterial taxa were not specialized to a specific rainfall regime. Similarly, in other soil reciprocal transplant experiments, fairly modest shifts in bacterial community composition were observed over 1 to 17 year study intervals when climate envelopes in home and away treatments were relatively similar (Lazzaro, Gauer and Zeyer 2011; Bond-Lamberty et al.2016). For example, reciprocal transplants of soil cores between forest stand types or soil types had no effect on community composition (Hannam, Quideau and Kishchuk 2007; Bond-Lamberty et al.2016). Bottomley et al. (2006) found that bacterial communities in meadow soils were largely unchanged following transplantation to a forest, although the reverse treatment did modify community structure. By contrast, moving soils between oak and grassland rapidly altered oak microbial communities that experienced hotter, drier conditions outside their historical climate (Waldrop and Firestone 2006). Microbial community composition also responded rapidly to transplantation along a dramatic (4.5°C) temperature gradient (Vanhala et al.2011). In summary, we suspect that more extreme climatic shifts (relative to natural climate variability) may be required to alter the relative abundance of generalist bacterial taxa. The compositional resistance of bacterial communities in both common gardens helps explain why patterns of respiration across soil types were unaltered by the transplant and dispersal treatments. Differences in community structure among common gardens, sites of soil origin and transplant treatments were driven by species turnover, not nestedness (local extinction), suggesting that changes in the identity or relative abundance of taxa drove the subtle compositional responses to the experimental treatments. Local specialization is a potential ecological mechanism for these compositional changes following transplantation. In support of this idea, OTUs with the highest ecological specificity tended to be more abundant in the garden matching their home precipitation region (drier west vs. wetter east), and were highest for the garden soils that were located at their specific site of origin. Moreover, decreases in the abundance of specialist OTUs following transplantation to a novel precipitation region (without changes in the chemical composition of their surrounding soil matrix) suggest that many of these taxa were responding directly to water availability. These results are consistent with laboratory studies on cultured microbes that have demonstrated ecological specialization related to moisture niche (Lennon et al.2012). Additionally, field studies have found strong niche conservatism in microbial responses to rapid change in soil moisture (Placella, Brodie and Firestone 2012). In our study, over 40% of ecological specialists were Actinobacteria, which are often found to be especially responsive to drought (Barnard, Osborne and Firestone 2013; Evans, Wallenstein and Burke 2014). This may explain why ecological specialization indices were higher in soils from drier western sites and in common garden (WFC and MOR) soils, where Actinobacteria comprised a relatively greater proportion of observed taxa. Despite the fact that most taxa were habitat generalists, specialists tended to be the most locally abundant OTUs. Specialization indices were calculated with respect to site of soil origin, not MAP alone, and so may also capture microbial responses to multiple non-climatic site factors. Soil characteristics are often reported as microbial drivers; however, soil organic C content, texture and pH explained only 2% of the variation in community structure. Although it is likely that unmeasured soil factors also shape bacterial communities, soil properties do not appear to be major controllers in the current system. Our results directly contrast with other studies that report major effects of pH and carbon availability on microbial community structure (Lauber et al.2008; Fierer et al.2009), perhaps because the relatively uniform geology of the Edwards Plateau constrains the amount of variation in these properties. We found no evidence that enhanced dispersal of microbes from the surrounding environment substantially affected bacterial community structure, similar to previous studies demonstrating a lack of dispersal limitation in soil bacteria (van der Gucht et al.2007). Cores subject to the enhanced dispersal treatment tended to have slightly more diverse communities based on number of taxa, but this response was not observed in all cases. Although these results are consistent with high microbial resistance, it is also possible that the microbial inoculation treatments failed due to environmental conditions: both common gardens experienced a historic drought at the start of the experiment, and bacteria in the soil inoculum may not have survived long under exceedingly hot, dry conditions. Alternatively, successful establishment of immigrating taxa—rather than immigration rate itself—may be the limiting factor for community-level shifts in response to environmental change (Fukami et al.2010). Establishment of new taxa may be impeded by biotic interactions with resident microbial species (Jones et al.2017), or by edaphic conditions unrelated to rainfall. It is possible that these issues were further exacerbated by low microbial numbers in our inoculum treatments; however, given typical microbial biomass values for these soils, each of the four repeated inoculum treatments likely delivered many thousands of viable cells. CONCLUSIONS Historical rainfall regime shapes bacterial community structure and function along a steep precipitation gradient (Hawkes et al.2017). We found that differences in community composition among sites are driven by species turnover rather than nestedness, consistent with some degree of ecological specialization based on precipitation regime. Although specialist taxa tended to exhibit more pronounced changes in abundance when transplanted to regions of higher or lower rainfall, compositional and functional responses to the transplant treatment were minimal. This suggests that the habitat generalists which dominated these microbial communities confer a high degree of ecological resistance to altered environmental conditions, with potentially important consequences for ecosystem-scale carbon cycling. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements We are grateful for assistance with fieldwork and labwork provided by Jim Dula, Nathan Christianson and Sarah Chen. Site access was kindly permitted by the Texas Parks and Wildlife Department, the Texas Ecolab Program and the Texas Historical Commission. Previous drafts of the manuscript were improved by comments from Jennifer Rocca and Stephanie Kivlin. FUNDING This material is based upon work supported by Texas Ecolab and the University of Texas funds to CVH. BGW was supported by a NSF Graduate Research Fellowship (DGE-1110007) and a NSF Doctoral Dissertation Improvement Grant (DEB-1210361). Conflict of interest. None declared. REFERENCES Barnard R, Osborne C, Firestone M. Responses of soil bacterial and fungal communities to extreme desiccation and rewetting. ISME J  2013; 7: 2229– 41. Google Scholar CrossRef Search ADS PubMed  Baselga A. Partitioning the turnover and nestedness components of beta diversity. Global Ecol Biogeogr  2009; 19: 134– 43. Google Scholar CrossRef Search ADS   Baselga A, Orme CDL. betapart: an R package for the study of beta diversity. Methods Ecol Evol  2012; 3: 808– 12. Google Scholar CrossRef Search ADS   Bond-Lamberty B, Bolton H, Fansler S et al.   Soil respiration and bacterial structure and function after 17 years of a reciprocal soil transplant experiment. PLoS One  2016; 11: e0150599. Google Scholar CrossRef Search ADS PubMed  Borcard D, Legendre P, Drapeau P. Partialling out the spatial component of ecological variation. Ecology  1992; 73: 1045– 55. Google Scholar CrossRef Search ADS   Bottomley PJ, Yarwood RR, Kageyama SA et al.   Responses of soil bacterial and fungal communities to reciprocal transfers of soil between adjacent coniferous forest and meadow vegetation in the Cascade Mountains of Oregon. Plant Soil  2006; 289: 35– 45. Google Scholar CrossRef Search ADS   Cadotte MW. Metacommunity influences on community richness at multiple spatial scales: a microcosm experiment. Ecology  2006; 87: 1008– 16. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Cline LC, Zak DR. Dispersal limitation structures fungal community assembly in a long-term glacial chronosequence. Env Microbiol Rep  2013; 16: 1538– 48. Google Scholar CrossRef Search ADS   deVries FT, Shade A. Controls on soil microbial community stability under climate change. Front Microbiol  2013; 4: 265. Google Scholar PubMed  Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics  2010 26: 2460– 1. Google Scholar CrossRef Search ADS PubMed  Evans S, Wallenstein M, Burke I. Is bacterial moisture niche a good predictor of shifts in community composition under long-term drought? Ecology  2014; 95: 110– 22. Google Scholar CrossRef Search ADS PubMed  Fierer N, Strickland MS, Liptzin D et al.   Global patterns in belowground communities. Ecol Lett  2009; 12: 1238– 49. Google Scholar CrossRef Search ADS PubMed  Freedman Z, Zak DR. Soil bacterial communities are shaped by temporal and environmental filtering: evidence from a long-term chronosequence. Env Microbiol Rep  2015; 17: 3208– 18. Google Scholar CrossRef Search ADS   Fukami T, Dickie IA, Paula Wilkie J et al.   Assembly history dictates ecosystem functioning: evidence from wood decomposer communities. Ecol Lett  2010; 13: 675– 84. Google Scholar CrossRef Search ADS PubMed  Gianuca AT, Declerck SAJ, Lemmens P et al.   Effects of dispersal and environmental heterogeneity on the replacement and nestedness components of β-diversity. Ecology  2017; 98: 525– 33. Doi: /101002/ecy1666. Google Scholar CrossRef Search ADS PubMed  Hannam KD, Quideau SA, Kishchuk BE. The microbial communities of aspen and spruce forest floors are resistant to changes in litter inputs and microclimate. Appl Soil Ecol  2007; 35: 635– 47. Google Scholar CrossRef Search ADS   Hawkes CV, Waring BG, Rocca JD et al.   Historical climate controls soil respiration responses to current soil moisture. P Natl Acad Sci USA  2017; 114: 6322– 7. Google Scholar CrossRef Search ADS   Hawkes CV, Keitt TH. Resilience vs historical contingency in microbial responses to environmental change. Ecol Lett  2015; 18: 612– 25. Google Scholar CrossRef Search ADS PubMed  Hughes J, Hellmann J, Ricketts T et al.   Counting the uncountable: statistical approaches to estimating microbial diversity. Appl Environ Microb  2001; 67: 4399. Google Scholar CrossRef Search ADS   Jones ML, Ramoneda J, Rivett DW et al.   Biotic resistance shapes the influence of propagule pressure on invasion success in bacterial communities. Ecology  2017; 98: 1743– 9. Google Scholar CrossRef Search ADS PubMed  Kindt R, Coe R. Tree Diversity Analysis: A Manual and Software for Common Statistical Methods for Ecological and Biodiversity Studies . Nairobi, Kenya: World Agroforestry Centre, 2005. Kivlin SN, Hawkes CV. Temporal and spatial variation of soil bacteria richness, composition, and function in a Neotropical rainforest. PLoS One  2016; 11: e0159131. Available at: https://doi.org/10.1371/journal.pone.0159131. Google Scholar CrossRef Search ADS PubMed  Klindworth A, Pruesse E, Schweer T et al.   Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res  2013; 41. doi: 101093/nar/gks808. Lauber CL, Strickland MS, Bradford MA et al.   The influence of soil properties on the structure of bacterial and fungal communities across land-use types. Soil Biol Biochem  2008; 40: 2407– 15. Google Scholar CrossRef Search ADS   Lawrence D, Bell T, Barraclough TG. The effect of immigration on the adaptation of microbial communities to warming. Am Nat  2016; 187: 236– 48. Google Scholar CrossRef Search ADS PubMed  Lazzaro A, Gauer A, Zeyer J. Field-scale transplantation experiment to investigate structures of soil bacterial communities at pioneering sites. Appl Environ Microb  2011; 77: 8241– 8. Google Scholar CrossRef Search ADS   Legendre P, Borcard D, Peres-Neto PR. Analyzing beta diversity: partitioning the spatial variation of community composition data. Ecol Monogr  2005; 75: 435– 50. Google Scholar CrossRef Search ADS   Lennon J, Aanderud Z, Lehmkuhl B et al.   Mapping the niche space of soil microorganisms using taxonomy and traits. Ecology  2012; 93: 1867– 79. Google Scholar CrossRef Search ADS PubMed  Lindström ES, Östman Ö. The importance of dispersal for bacterial community composition and functioning. PLoS One  2011; 6: e25883. Google Scholar CrossRef Search ADS PubMed  Lumley T. leaps: regression subset selection, including exhaustive search. 2017. https://cranr-projectorg. Version: 3.0. Mariadassou M, Pichon S, Ebert D. Microbial ecosystems are dominated by specialist taxa. Ecol Lett  2015; 18: 974– 82. Google Scholar CrossRef Search ADS PubMed  McCune B, Mefford MJ. PC-ORD: Multivariate Analysis of Ecological Data . Gleneden Beach, Oregon, USA: MjM Software, 2011. Moyano FE, Manzoni S, Chenu C. Responses of soil heterotrophic respiration to moisture availability: an exploration of processes and models. Soil Biol Biochem  2013; 59: 72– 85. Google Scholar CrossRef Search ADS   Oksanen J, Blanchet FG, Friendly M et al.   Package ‘vegan.’ 2016. https://cranr-projectorg. Version: 2.4-1. Placella SA, Brodie EL, Firestone MK. Rainfall-induced carbon dioxide pulses result from sequential resuscitation of phylogenetically clustered microbial groups. P Nat Acad Sci USA  2012; 109: 10931– 6. Google Scholar CrossRef Search ADS   Quast C, Pruesse E, Yilmaz P et al.   The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res  2012; 41: D590– 6. Google Scholar CrossRef Search ADS PubMed  Ulrich W, Almeida-Neto M, Gotelli NJ. A consumer's guide to nestedness analysis. Oikos  2009; 118: 3– 17. Google Scholar CrossRef Search ADS   van der Gucht K, Cottenie K, Muylaert K et al.   The power of species sorting: local factors drive bacterial community composition over a wide range of spatial scales. P Natl Acad Sci USA  2007; 104: 20404– 9. Google Scholar CrossRef Search ADS   van Tienderen PH. Generalists, specialists, and the evolution of phenotypic plasticity in sympatric populations of distinct species. Evolution  1997; 51: 1272– 380. Google Scholar CrossRef Search ADS   Vanhala P, Karhu K, Tuomi M et al.   Transplantation of organic surface horizons of boreal soils into warmer regions alters microbiology but not the temperature sensitivity of decomposition. Glob Change Biol  2011; 17: 538– 50. Google Scholar CrossRef Search ADS   Waldrop MP, Firestone MK. Response of microbial community composition and function to soil climate change. Microb Ecol  2006; 52: 716– 24. Google Scholar CrossRef Search ADS PubMed  Wang Q, Garrity GM, Tiedje JM et al.   Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microb  2007; 73: 5261– 6267. Google Scholar CrossRef Search ADS   © FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png FEMS Microbiology Ecology Oxford University Press

Ecological mechanisms underlying soil bacterial responses to rainfall along a steep natural precipitation gradient

Loading next page...
 
/lp/ou_press/ecological-mechanisms-underlying-soil-bacterial-responses-to-rainfall-uaoWd58jwq
Publisher
Blackwell
Copyright
© FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
0168-6496
eISSN
1574-6941
D.O.I.
10.1093/femsec/fiy001
Publisher site
See Article on Publisher Site

Abstract

Abstract Changes in the structure and function of soil microbial communities can drive substantial ecosystem feedbacks to altered precipitation. However, the ecological mechanisms underlying community responses to environmental change are not well understood. We used an 18-month soil reciprocal transplant experiment along a steep precipitation gradient to quantify how changes in rainfall affected bacterial community structure. We also conducted an enhanced dispersal treatment to ask whether higher immigration rates of taxa from the surrounding environment would accelerate community responses to climate change. Finally, we addressed how the composition of soil bacteria communities was related to the functional response of soil respiration to moisture in these treatments. Bacterial community structure (OTU abundance) and function (respiration rates) changed little in response to manipulation of either rainfall environment or dispersal rates. Although most bacteria were ecological generalists, a subset of specialist taxa, over 40% of which were Actinobacteria, tended to be more abundant in the rainfall environment that matched their original conditions. Bacteria community composition was an important predictor of the respiration response to moisture. Thus, the high compositional resistance of microbial communities dictated respiration responses to altered rainfall in this system. dispersal, drought, generalist, nestedness, reciprocal transplant, spatial turnover, specialist INTRODUCTION Soil microbes mediate enormous fluxes of carbon (C) among plant detritus, soil and the atmosphere via their roles in organic matter decomposition and soil respiration. Therefore, microbial responses to altered temperature, precipitation or nutrient inputs may drive substantial feedbacks to global climate change. We previously found that microbial respiration was a persistent function of historical climate rather than current conditions: soils from historically wetter sites exhibited higher respiration rates regardless of contemporary water availability (Hawkes et al.2017). Such functional legacies are likely created by microbial community resistance to environmental change (Hawkes and Keitt 2015), which has been observed in other long-term studies (e.g. Bond-Lamberty et al.2016). Nevertheless, high beta diversity in microbial communities across sites or through time is commonly reported (Freedman and Zak 2015; Kivlin and Hawkes 2016), consistent with microbial communities rapidly tracking environmental conditions. Here, we suggest that a better understanding of microbial community responses to the environment can be achieved by considering community structure beyond diversity or composition to include species nestedness, turnover and specialization. Most studies examining microbial responses to climate change focus on patterns in alpha and beta diversity. Soil bacterial communities exhibit extremely high species richness (alpha diversity) at the local scale and dramatic compositional turnover (beta diversity) among sites (Hughes et al.2001; van der Gucht et al. 2007). Few studies, however, have separately examined two independent components of microbial community beta diversity: spatial turnover and nestedness. Species turnover among sites is often equated with beta diversity, but nestedness is a component of beta diversity where the species pool at one local site is a subset of the species pool at a more diverse site (Baselga 2009). Partitioning the turnover and nestedness components of beta diversity is important because these factors respond to different ecological drivers. For example, community nestedness is highly sensitive to variation in local extinction rates and is often associated with environmental harshness (Ulrich, Almeida-Neto and Gotelli 2009). High species turnover, in contrast, is often driven by species sorting over environmental gradients (van der Gucht et al.2007). Whether spatial turnover or nestedness dominates community responses to climate change may depend on the interaction between two factors: habitat specialization within the microbial community and rates of dispersal. A microbial species can be considered a ‘habitat specialist’ if it is encountered more often in one type of habitat than any other (Mariadassou, Pichon and Ebert 2015). If the regional species pool is primarily composed of habitat specialists, then spatial turnover across environmental gradients should be high (Baselga 2009). By contrast, if habitat generalists predominate in the regional species pool, then nestedness may account for more of the variation in community composition among local sites, due to stochastic extinction of particular taxa at individual sites (Ulrich, Almeida-Neto and Gotelli 2009). In either scenario, high rates of dispersal can homogenize community structure across environmental gradients (Cadotte 2006). Dispersal may speed community acclimatization to altered environments if immigrants are better suited to contemporary conditions than residents (van der Gucht et al.2007; deVries and Shade 2013; Gianuca et al.2017). In this scenario, dispersal limitation would underlie microbial community resistance to change (Cline and Zak 2013; Lawrence, Bell and Barraclough 2016). Few studies have quantified the relative abundance of habitat specialists and generalists in soil microbial communities, although model simulations suggest that ecological specialization may drive community responses to environmental change (Hawkes and Keitt 2015). Evolutionary models demonstrate that large differences in the ‘optimal’ phenotype across different habitats tend to favor the evolution of ecological specialists, whereas generalists are expected in conditions that promote high phenotypic plasticity (van Tienderen 1997). High beta diversity of microbial assemblages suggests that most soil bacteria are specialists, yielding a strong positive correlation between taxon abundance and ecological specificity at broad spatial scales (Mariadassou, Pichon and Ebert 2015). Similarly, the degree of dispersal limitation within soil microbial communities is largely unquantified. Laboratory studies suggest that variation in dispersal rates can affect bacterial community structure and function (Lindström and Östman 2011), but there is little evidence for dispersal limitation in intact bacterial communities (van der Gucht et al.2007). Our goal in this study was to experimentally test microbial community responses to altered precipitation, focusing on how dispersal and habitat specialization shape these responses. To do so, we conducted a soil reciprocal transplant experiment across a 45-cm rain gradient in central Texas to expose soils from three dry and three wet sites to either to their normal ‘home’ or a changed ‘away’ precipitation regime (Fig. 1). In addition, half of the soil cores were subjected to elevated dispersal rates (i.e. inoculation) of microbes from the surrounding precipitation region. The experiment was designed to test four hypotheses: First, we expected to find differences in microbial community composition based on soil origin, driven by differences in historical climate and edaphic conditions at each unique site (Hypothesis 1). Next, we anticipated that the degree of microbial habitat specialization would structure communities both along the precipitation gradient and in response to transplantation treatments. If specialists are more abundant than generalists, spatial turnover should be high across the rainfall gradient; if generalists are more abundant, nestedness should play a larger role in generating beta diversity (Hypothesis 2). We also hypothesized that the enhanced dispersal treatment would reduce beta diversity (both turnover and nestedness) among soils of different origins (Hypothesis 3). Finally, we expected that changes in microbial community structure in response to transplantation and dispersal would drive the patterns of respiration in each community (Hypothesis 4) that were previously observed (Hawkes et al.2017). Figure 1. View largeDownload slide Schematic of the field reciprocal transplant experiment. Figure 1. View largeDownload slide Schematic of the field reciprocal transplant experiment. METHODS Study site and experimental design We worked in sites located across the Edwards Plateau in central Texas, where mean annual precipitation (MAP) decreases steeply from east to west (∼90 to 45 cm) while soil parent material and savanna grassland plant communities (dominated by Juniperus ashei, Bouteloua spp., Aristida spp. and Schizachyrium scoparium) remain relatively invariant. At these sites, mean rainfall peaks in early spring (April/May) and in the fall (September/October), but year-to-year variance is high. Mean annual temperature is approximately 20°C. Sample locations within sites were selected to minimize differences in topography (<2% slope) and the plant community (native grassland with <50% woody plant cover), and all soils are rocky, calcareous Mollisols. We conducted a reciprocal transplant experiment of soil cores from the drier west end to the wetter east end of the gradient and vice versa. To do so, we established two 6 m × 6 m common garden plots, one located at an eastern site where MAP = 89 cm (Ladybird Johnson Wildflower Center, or WFC, Lat: 30.18, Lon: –97.88) and one located 300 km to the west where MAP = 63 cm (MOR Ecolab, Lat: 29.44, Lon: –100.06) (Table S1, Supporting Information). Soils from three wetter eastern sites and three drier western sites were placed at each common garden in PVC cores (n = 6), creating a full factorial combination of home region and transplanted treatments. Moreover, six cores containing soils sampled on-site were also incubated in each garden as a control treatment (Fig. 1). Soils from all eight sites (three eastern sites, three western sites and the two common garden sites) were collected in July 2011, sieved to 2 mm to remove rocks and roots, packed into PVC cores at approximately native bulk density (150 g at ∼0.95 g cm−3) and installed in each common garden site within a week. Cores were arranged in rows, 50 cm apart, in random order with respect to treatment. The PVC cores measured 5 cm in diameter and 15 cm in depth, with 38-μm mesh bottoms to hold soil in while allowing for drainage. All cores were open at the top, allowing for passive dispersal of microbial propagules from the surrounding environment. The 38-μm pore size of the mesh core bottom meant that bacteria and fungi could also enter from surrounding soil. To examine the effects of dispersal rate on microbial community responses to the transplant treatment, half the cores at each garden (i.e. n = 3 in each treatment combination) received an ‘active dispersal’ treatment. We intended to simulate enhanced immigration rates of taxa from the surrounding precipitation region. Therefore, cores receiving the active dispersal treatment in the eastern common garden were inoculated with ‘eastern’ microbial assemblages, and cores in the western common garden were inoculated with ‘western’ assemblages. To create the inoculants, soils were freshly collected from the three sites in the surrounding precipitation region and the common garden site itself, thoroughly mixed in sterile water in equal proportions and filtered through 38-μm mesh to remove solid particles; 500 μL of the resulting microbial slurry was applied to each core in the enhanced dispersal treatment. A corresponding volume of sterile water was applied to the remaining experimental units. Enhanced dispersal cores were inoculated at establishment in June 2011, and again in November 2011, April 2012 and October 2012 to mimic naturally occurring variation in arrival rates (Jones et al.2017). With six experimental replicates for each factorial combination of common garden, soil origin and dispersal treatments, there were a total of 72 cores at each common garden site in addition to the six control cores (n = 78 per garden, 156 for the entire experiment). We measured a suite of physicochemical properties for each of the eight soil types included in our reciprocal transplant experiment. Soil C and N content were measured via combustion (Apollo TOC/TN, Teledyne Tekmar, Mason, OH, USA). Soil texture (% sand, silt and clay) was determined via the hydrometer method. Soil pH was measured in by mixing samples with distilled water in a 1:2 ratio. A summary of soil chemical data is given in Table S1, Supporting Information. Characterization of microbial community structure Cores were harvested in December 2012, after 18 months in the field, at which time they were removed from the ground intact and transported to the laboratory on ice. Soils from each core were sieved to 2-mm and homogenized prior to DNA extraction with a MoBio PowerSoil extraction kit (MoBio Laboratories, Carlsbad, CA). DNA extracts were amplified with bacteria-specific SSU primers Bact-0341 and Bact-0785 (Klindworth et al.2012), and the reverse primer was barcoded with an 11-bp treatment-specific barcode sequence. Amplicons were cleaned using the MoBio UltraClean PCR Clean-Up kit (MoBio Laboratories, Carlsbad, CA), quantified via Qubit dsDNA High Sensitivity assays (Life Technologies, Grand Island, NY) and pooled in equimolar mixtures. Pooled libraries were purified of fragments <150 bp in size using Agencourt AmPure XP magnetic beads (Beckman Coulter, Brea, CA). Libraries were sequenced on a 454 GS FLX sequencer with titanium chemistry (Roche, Branford, CT) at the University of Texas Genome Sequencing and Analysis Facility. We analyzed the resulting sequences using the QIIME pipeline v. 1.5.0 (Caporaso et al.2010). We eliminated sequences with quality scores below 25, fewer than 150 bases, homopolymer runs of more than 10 bp, with any ambiguous bases, or with an anomalous barcode. We also removed all singletons (unique to the entire dataset) and eliminated chimeras with USEARCH 6.1 (Edgar 2010). We used the UCLUST algorithm (Edgar 2010) to define operational taxonomic units (OTUs) at 97% sequence similarity. Representative sequences from each OTU were searched against the SILVA (Quast et al.2012) and RDP (Wang et al.2007) databases to identify and eliminate any non-bacterial taxa. To avoid bias in sampling effort, we rarefied to 1000 sequences, which was the smallest number of sequences in a sample. Sequences reads are archived in the NCBI Sequence Read Archive, Accession PRJNA368746. Quantifying microbial functional responses to soil moisture To determine whether the transplant or dispersal treatments affected microbial functional responses to water availability, soil respiration rate in the field common gardens was measured under ambient moisture conditions in November 2011, April 2012, October 2012 and December 2012. We capped cores with septa-embedded lids and sampled the headspace after 30 min. Headspace samples were stored in gas-tight vials until CO2 was quantified on a SRI 8610C gas chromatograph with FID and methanizer (SRI Instruments, Santa Monica, CA). We also examined soil respiration in response to controlled moisture conditions in the lab using data from Hawkes et al. (2017). Briefly, soils from the common garden cores were subsampled at harvest and assigned to high (25%) and low (10%) gravimetric soil moisture treatments in a laboratory incubation experiment, with microbial function (CO2 production) measured over 8 weeks. Headspace samples were analyzed for CO2 content as described above for field samples. Composition statistical analysis To examine microbial community structure and diversity across all experimental cores (Hypothesis 1), we used permutational ANOVAs to examine variation in bacterial community composition among treatments. Henceforth, we will use the term ‘region of origin’ to indicate the climatic province from which soils were originally sampled (wetter east vs. drier west), and ‘site of origin’ to indicate the specific site from which they were sampled (Table S1, Supporting Information). Therefore, in our models, site is nested within region of origin. Permutational ANOVAs were conducted using a combination of adonis in the vegan package (Oksanen et al.2016), with each term entered last in the model to obtain Type III sums of squares, and nested.npmanova in the BiodiversityR package (Kindt and Coe 2005) to obtain correct F-ratios and P values; 999 permutations were used per run. We also conducted planned pairwise comparisons among four treatment groups: soils that originated from the wetter eastern end of the gradient, incubating in home and away rainfall environments; and soils that originated from the drier western end of the gradient, incubating at home and away. A Bonferroni correction was applied to control for multiple comparisons. To visualize differences in community composition, we calculated non-metric multidimensional scaling (NMDS) coordinates based on Bray-Curtis dissimilarity matrices in PC-ORD (v6.2; McCune and Mefford 2011). Although statistics reported here are based on relative abundance data, all effect sizes and P values are nearly identical for presence-absence data. To examine variation in bacterial community composition in relation to soil properties, historical precipitation regime and recent drought stress, we used partial regularized discriminant analysis (RDA; Borcard, Legendre and Drapeau 1992; Legendre, Borcard and Peres-Neto 2005). Collinear variables were removed based on variance inflation factors >2.5, resulting in the following predictors in the final model: latitude, elevation, MAP, MAT, soil moisture in each core at harvest, the Palmer Drought Severity Index (PDSI) for 18 months in the home site prior to collection, PDSI for the experimental duration in the common garden sites, soil total organic carbon content (TOC), soil pH and soil % clay. Ecological specificity analysis We calculated the ecological specificity index (SI) following Mariadassou, Pichon and Ebert (2015) for the 1000 most abundant OTUs across the entire experiment. The index ranges between 0 and 1, with 0 indicating a perfect generalist (equal abundance across all sites) and 1 indicating a perfect specialist (found only in one of the eight sites). Specificity was calculated relative to the site of soil origin in which that OTU was most abundant; that is, we calculated local abundance in site X relative to average local abundance across all sites of soil origin, combined across gardens and dispersal treatments. Following these calculations, we performed a one-way ANOVA to determine whether the distribution of SIs varied among sites of soil origin. We also calculated a community-weighted mean SI for each site of soil origin, based on the presence or absence of each OTU within cores belonging to each site of soil origin. To determine whether the relative abundance of specialists or generalists shaped community responses to the transplant treatment (Hypothesis 2), we identified ecological specialists as those in the top decile of the SI distribution and ecological generalists as those in the bottom decile. Then, for each taxon, we calculated log response ratios (LRRs) for abundance in home vs. away transplant common gardens. We used a t test to ask whether the LRRs were greater for specialists vs. generalists, indicating higher abundance of specialists in their home precipitation region. Beta diversity partitioning We quantified community dissimilarity due to the turnover vs. nestedness components of beta diversity using the betapart package in R (Baselga and Orme 2012). We used the beta.pair function to generate dissimilarity matrices for the turnover and nestedness components of beta diversity across all experimental cores. Multi-response Permutation Procedure (MRPP) tests were conducted to determine whether dissimilarity associated with turnover and nestedness varied among soil sites of origin, transplant treatments and dispersal treatments (Hypothesis 3). To understand how turnover and nestedness were related to any shift in composition associated with transplantation, we also quantified these components of beta diversity across cores in home vs. away treatments within each site of soil origin. Finally, turnover and nestedness within each soil type were regressed against community mean SI indices to determine whether a higher relative abundance of ecological specialists is correlated with greater species turnover in response to the transplant treatment. Linking microbial community composition with functional responses to moisture To determine whether microbial community structure was a significant predictor of patterns of respiration in response to current moisture (Hypothesis 4), we used a regression approach to identify drivers of the respiration moisture response in the lab incubation experiment. Laboratory data were chosen because moisture conditions were more carefully controlled than in the field, but differences in CO2 fluxes between eastern- and western-origin soils were identical to those observed in situ (see ‘Results’). First we calculated ΔResp (microbial functional response to soil moisture) as the difference in mean CO2 flux in high vs. low soil moisture regimes in the laboratory incubation experiment. This calculation was performed for each of the 26 unique treatment groups (unique combination of soil site of origin × common garden location × dispersal treatment; and two common garden controls). Then we modeled ΔResp as a function of climatic, edaphic and microbial community variables. For our climatic predictors at each site of origin, we included the total rainfall experienced the year before transplant (June 2010 – 2011), henceforth referred to as ‘antecedent’ precipitation, and mean temperature at each site of origin during that time period (‘antecedent’ temperature). We excluded 30-year MAP and MAT because they were highly correlated with multiple predictor variables; see Table S3, Supporting Information. Because there was no significant main effect of common garden identity in analyses of bacteria diversity and composition, we did not include climate variables associated with the transplant period in the gardens. To explore the effect of soil physicochemical traits on respiratory responses, we used total organic C, pH and % clay in soils from each site of origin (total N was excluded because was highly correlated with soil C: R2 = 0.95). Finally, we included four predictor variables related to microbial community composition and structure: mean bacterial community richness and phylogenetic diversity in each treatment group, and mean NMS axis 1 and 2 scores in each treatment group. The latter were derived from the Bray-Curtis distance matrix of bacterial community composition. We used the regsubsets command in the R package leaps (Lumley 2017) to iterate through all possible combinations of predictors. The best models were selected by comparing BIC scores; pseudo R2 for each generalized linear model were calculated by comparing residual deviance vs. null deviance. RESULTS Diversity and community composition During the 18 months that soil cores were incubated in the two common gardens, cumulative rainfall was 54% higher in the eastern vs. western garden location (Fig. S1, Supporting Information), although both gardens experienced a 5-month severe drought at the start of the experiment. At the time of harvest, soil bacterial communities included 734 OTUs per sample on average (after rarefaction to 1000 sequences). However, alpha diversity was significantly affected by the interaction of site of soil origin and dispersal treatment (Table S2, Supporting Information): in soils originating from four sites, alpha diversity increased modestly by 2%–4% in the active dispersal treatment. All communities were generally dominated by Acidobacteria, Actinobacteria and Proteobacteria (Fig. S2a, Supporting Information). Bacterial community composition varied primarily by site of soil origin, explaining about 13% of the compositional variance despite 18 months of incubation in a new location (Table 1, Fig. 2). Site origin effects differed slightly between the two common gardens, but this explained only 2.6% of compositional variation (Table 1). The dispersal treatment had no effect on community composition, nor did region (drier west vs. wetter east) of soil origin. In planned pairwise comparison tests, eastern-origin and western-origin bacteria communities were always significantly different from each other (P < 0.01). However, there was no difference in community structure between eastern-origin soils in their home vs. away rainfall environment (F1,73 = 1.14, P = 0.294), and the same was true for western soils incubating at home vs. away (F1,73 = 1.14, P = 0.354). This indicates that bacterial community structure changed little in response to transplantation. Figure 2. View largeDownload slide NMS ordination of bacterial community composition in soils from each site, garden and dispersal treatment. (Stress = 22.44 after 56 iterations.) Only site and the interaction of site with common garden had significant effects. The dashed line separates soils originating from eastern and western regions. Average NMS scores for each site are plotted ±1SE (n = 4–6). Figure 2. View largeDownload slide NMS ordination of bacterial community composition in soils from each site, garden and dispersal treatment. (Stress = 22.44 after 56 iterations.) Only site and the interaction of site with common garden had significant effects. The dashed line separates soils originating from eastern and western regions. Average NMS scores for each site are plotted ±1SE (n = 4–6). Table 1. Results of permutational ANOVA examining bacteria community composition as a function of region of soil origin (Region), common garden location (Garden), dispersal treatment (Dispersal), soil moisture treatment (Moisture), site nested in region (Site|Region) and their interactions. Factors significant at P < 0.05 are in bold.   df  MS  F  R2  P  Region  1  1.516  1.446  0.031  0.163  Garden  1  0.457  1.426  0.009  0.918  Dispersal  1  0.322  1.029  0.007  0.993  Region*Garden  1  0.326  1.016  0.007  0.450  Region*Dispersal  1  0.300  0.959  0.006  0.557  Garden*Dispersal  1  0.307  0.993  0.006  0.998  Region*Garden*Dispersal  1  0.300  0.970  0.006  0.740  Site(Region)  6  1.048  3.446  0.128  0.001  Garden*Site(Region)  4  0.320  1.053  0.026  0.037  Dispersal*Site(Region)  4  0.313  1.030  0.026  0.183  Garden*Dispersal*Site(Region)  4  0.309  1.017  0.025  0.359  Error  117  0.305  0.724        df  MS  F  R2  P  Region  1  1.516  1.446  0.031  0.163  Garden  1  0.457  1.426  0.009  0.918  Dispersal  1  0.322  1.029  0.007  0.993  Region*Garden  1  0.326  1.016  0.007  0.450  Region*Dispersal  1  0.300  0.959  0.006  0.557  Garden*Dispersal  1  0.307  0.993  0.006  0.998  Region*Garden*Dispersal  1  0.300  0.970  0.006  0.740  Site(Region)  6  1.048  3.446  0.128  0.001  Garden*Site(Region)  4  0.320  1.053  0.026  0.037  Dispersal*Site(Region)  4  0.313  1.030  0.026  0.183  Garden*Dispersal*Site(Region)  4  0.309  1.017  0.025  0.359  Error  117  0.305  0.724      View Large In the RDA analysis, all environmental and spatial variables together explained only 10.7% of the variance in community composition (Table 2), despite wide variation in soil chemical properties (soil % C ranging from 2.4% to 10.9%; clay ranging from 18.5% to 39.7%; pH ranging from 7.4 to 8.2; Table S1, Supporting Information). Precipitation-related variables (MAP, gravimetric soil moisture and PDSI) accounted for 2.7%, soil characteristics (TOC, pH and clay) explained 2.1% and spatial effects explained only 0.1% of the variation in bacteria community composition among sites. Table 2. Results of partial RDA parsing how bacteria community composition relates to climatic (Clim), edaphic (Soil) and spatial (Space) factors. Most of the variation in the bacteria community was unexplained.   Inertia  R2  Adj R2  P  All vars  23.45  0.157  0.107  <0.001  Clim  20.19  0.054  0.027  <0.001  Soil  11.69  0.042  0.021  <0.001  Space  4.49  0.016  0.009  <0.001  Clim|Space + Soil  8.65    0.012  <0.001  Soil|Space + Clim  1.28    0.005  <0.001  Space|Clim + Soil  1.84    0.001  <0.001  Unexplained  125.73          Inertia  R2  Adj R2  P  All vars  23.45  0.157  0.107  <0.001  Clim  20.19  0.054  0.027  <0.001  Soil  11.69  0.042  0.021  <0.001  Space  4.49  0.016  0.009  <0.001  Clim|Space + Soil  8.65    0.012  <0.001  Soil|Space + Clim  1.28    0.005  <0.001  Space|Clim + Soil  1.84    0.001  <0.001  Unexplained  125.73        View Large Beta diversity partitioning Across all cores in the experiment, differences in bacteria community structure attributable to site of soil origin and transplantation treatments were generated by species turnover, not nestedness (Table 3). Neither turnover nor nestedness responded to the dispersal treatment because this treatment had no effect on community composition. Similarly, within sites of soil origin, changes in community structure in response to transplantation were dominated by species turnover (Table 4). The degree of nestedness (i.e. local extinction) across transplant treatments was generally quite small. Table 3. Dissimilarity among communities in different treatments due to species turnover and nestedness. P-values are based on MRPP tests on community dissimilarity matrices. Bold type indicates significance at P < 0.05.   Observed δ, (significance of δ)    Dissimilarity due to species turnover Expected δ 0.644  Dissimilarity due to nestedness Expected δ 0.021  Region of soil origin  0.628 (0.001)  0.021 (1.000)  Site of soil origin  0.570 (0.001)  0.022 (0.990)  Common garden site  0.642 (0.008)  0.021 (0.708)  Dispersal treatment  0.645 (0.562)  0.021 (0.634)    Observed δ, (significance of δ)    Dissimilarity due to species turnover Expected δ 0.644  Dissimilarity due to nestedness Expected δ 0.021  Region of soil origin  0.628 (0.001)  0.021 (1.000)  Site of soil origin  0.570 (0.001)  0.022 (0.990)  Common garden site  0.642 (0.008)  0.021 (0.708)  Dispersal treatment  0.645 (0.562)  0.021 (0.634)  View Large Table 4. Community-weighted mean specialization indices (SIs) for bacterial communities in each soil type; and the turnover and nestedness components of beta diversity between home and away precipitation treatments within each soil type.   Community-weighted mean SI  Turnover  Nestedness  Eastern sites  COL  0.208  0.123  0.010  ING  0.229  0.142  0.005  SCO  0.201  0.123  0.035  WFC (common garden site)  0.246  NA  NA  Western sites  DRS  0.218  0.136  0.040  SCS  0.218  0.141  0.005  SIE  0.213  0.150  0.005  MOR (common garden site)  0.247  NA  NA    Community-weighted mean SI  Turnover  Nestedness  Eastern sites  COL  0.208  0.123  0.010  ING  0.229  0.142  0.005  SCO  0.201  0.123  0.035  WFC (common garden site)  0.246  NA  NA  Western sites  DRS  0.218  0.136  0.040  SCS  0.218  0.141  0.005  SIE  0.213  0.150  0.005  MOR (common garden site)  0.247  NA  NA  View Large Ecological specialization Ecological specialization indices (SI) varied widely among taxa, ranging from 0.030 to 0.688, but the mean value of 0.184 suggests most taxa were not strong habitat specialists. However, a significant positive relationship between the specificity index and taxon rank abundances (R2 = 0.06, P < 0.001; Fig. 3) indicates numerical dominance by the most specialized taxa. There was significant variation in the degree of ecological specificity relative to site of soil origin (one-way ANOVA, F7,985 = 19.70, P < 0.001). In post-hoc tests, the highest SIs were found for taxa specialized to both home common garden sites (MOR, WFC; P < 0.001; Fig. 4a). Specificity indices were also significantly stronger for taxa specialized to soils originating from the drier western regions vs. the wetter eastern region (t = 2.051, P = 0.041, Fig. 4b). Community-weighted mean specialization indices in each site of soil origin followed a similar pattern, and were 18% higher in home common garden soils than in eastern-origin soils sampled from the wetter end of the rainfall gradient (Table 4). There was no relationship among community-weighted mean SIs and either the turnover or nestedness components of beta diversity. Figure 3. View largeDownload slide Relationship between OTU abundance and ecological specificity index (SI). Figure 3. View largeDownload slide Relationship between OTU abundance and ecological specificity index (SI). Figure 4. View largeDownload slide Indicators of bacteria specialization. Mean ecological specialization indices (+1 SE) for the most abundant bacterial taxa for (a) soils from each site of origin (arranged by MAP from left to right) and (b) in each precipitation region (wetter east vs. drier west). SI indices range from 0 (no specificity) to 1 (total specificity, i.e. the taxon is found only at one site). (c) Log response ratios (LRRs) for abundance in the home precipitation region vs. away region for specialist vs. generalist taxa. Positive LRRs indicate the taxon in question is more abundant in its home precipitation region. Figure 4. View largeDownload slide Indicators of bacteria specialization. Mean ecological specialization indices (+1 SE) for the most abundant bacterial taxa for (a) soils from each site of origin (arranged by MAP from left to right) and (b) in each precipitation region (wetter east vs. drier west). SI indices range from 0 (no specificity) to 1 (total specificity, i.e. the taxon is found only at one site). (c) Log response ratios (LRRs) for abundance in the home precipitation region vs. away region for specialist vs. generalist taxa. Positive LRRs indicate the taxon in question is more abundant in its home precipitation region. Specialist bacterial taxa (mean SI of 0.502) tended to be more abundant in home vs. away environments in comparison with generalist taxa (mean SI of 0.062; t = –1.795, P = 0.074; Fig. 4c). The taxonomic affiliation of specialists and generalists also varied significantly (χ2 = 25.729, P = 0.003): the largest fraction of generalist OTUs (28%) belonged to Proteobacteria, whereas 42% of specialist OTUs belonged to Actinobacteria (Fig. S2b, Supporting Information). Links between composition and function As previously reported (Hawkes et al.2017), rates of respiration were higher in soils that originated from sites located in the wetter eastern vs. drier western end of the gradient, in addition to current soil moisture conditions. This pattern was observed for in situ respiration rates measured in the field (Table S4a and b, Supporting Information) as well as in the laboratory microcosm experiment, where soil moisture was carefully controlled (Hawkes et al.2017). Our model selection procedures indicated that climatic, edaphic and microbial community variables could predict the respiratory response to soil moisture (ΔResp) moderately well: pseudo-R2 ranged from 0.18 to 0.33 (Table 5). As expected, the explanatory power of each model generally increased with the number of predictor variables. However, the two best (lowest BIC) models included only antecedent rainfall and NMS axis 1 scores as predictors. The magnitude of the respiratory response to soil moisture increased with the amount of rainfall soils received in the past year (Fig. 5a) and increased with NMS axis 1 scores (Fig. 5b), which were linked with the abundance of taxa in the Acidobacteria and Actinobacteria. NMS axis 1 scores, antecedent rainfall, soil pH and soil clay content were the variables most frequently included in the top-ranked models (Table 5); ΔResp decreased with soil clay content (Fig. 5c) and had a very weak negative relationship with soil pH (data not shown). Figure 5. View largeDownload slide The relationship between functional responses to soil moisture in the laboratory incubation (ΔResp) and climatic, edaphic and microbial community variables: (a) annual rainfall soils were exposed to prior to the transplant treatment; (b) bacteria community composition (as summarized by NMS 1 scores) and (c) soil clay content in each site of origin. Figure 5. View largeDownload slide The relationship between functional responses to soil moisture in the laboratory incubation (ΔResp) and climatic, edaphic and microbial community variables: (a) annual rainfall soils were exposed to prior to the transplant treatment; (b) bacteria community composition (as summarized by NMS 1 scores) and (c) soil clay content in each site of origin. Table 5. The five best-ranked models (as determined with BIC scores) predicting soil respiration responses to moisture in the laboratory incubation experiment. ‘Rain’ and ‘temp’ refer to antecedent rainfall and temperature at each site of origin in the 12 months prior to the transplant treatment; clay refers to soil % clay; NMS1 and NMS2 refer to mean scores on the first and second axis of an NMDS of bacteria community composition.   Predictors    Pseudo  Model rank  Biotic  Climatic  Edaphic  BIC  R2  1    Rain    – 109.46  0.182  2  NMS1  Rain    – 107.61  0.226  3  NMS1    Clay; pH  – 106.82  0.296  4  NMS1  Temp  Clay; pH  – 104.62  0.324  5  NMS1  Rain; temp  Clay; pH  – 101.58  0.329    Predictors    Pseudo  Model rank  Biotic  Climatic  Edaphic  BIC  R2  1    Rain    – 109.46  0.182  2  NMS1  Rain    – 107.61  0.226  3  NMS1    Clay; pH  – 106.82  0.296  4  NMS1  Temp  Clay; pH  – 104.62  0.324  5  NMS1  Rain; temp  Clay; pH  – 101.58  0.329  View Large DISCUSSION We found evidence for consistent differences in bacterial community composition among sites from across a natural precipitation gradient; this compositional variation was weakly structured by soil properties as well as historical and recent rainfall (Hypothesis 1). Although most taxa were habitat generalists, ecological specialists were locally abundant and tended to be more responsive to the transplantation treatments. The distribution of these specialist taxa likely contributed to the high spatial turnover observed across experimental cores (Hypothesis 2). Contrary to Hypothesis 3, however, the enhanced dispersal treatment had no effect on bacterial community structure or patterns of beta diversity. Finally, we found strong evidence that bacterial community composition mediated functional responses to soil moisture variation (Hypothesis 4). The characteristic respiratory patterns of soils from the wetter eastern and drier western regions of the rain gradient were maintained over the 18-month experiment, paralleling the compositional resistance of bacterial communities to environmental change. The most striking pattern we observed in our study was the high degree of resistance of bacteria community composition to novel rainfall regimes. This compositional resistance affected functional responses to water availability: eastern-origin (wetter) and western-origin (drier) soils exhibited characteristic patterns in soil respiration, which were unchanged following transplant and dispersal treatments. These functional patterns were correlated with variation in bacterial community structure among sites of soil origin. However, the relationship between climate, bacteria composition and respiration was not straightforward in our system, with MAP at each site explaining only ∼3% of the community compositional variance, but recent rainfall driving ∼18% of the functional response. There may be physiological adaptation of microbial taxa independent of composition, or indirect effects of rainfall on substrate availability for microbial growth (Moyano, Manzoni and Chenu 2013). Additionally, soil fungi, which we did not measure, likely also contributed to the observed patterns of respiration. Finally, the strong relationships between annual precipitation and microbial community composition that were observed along the gradient in 2008 (Hawkes et al.2017) may be lost over time, reduced by the disturbance required to construct cores or diminished under conditions of drought (we do not have initial microbial community data from the date of transplantation to address this issue definitively). Based on our results, the degree of environmental change relative to long-term climate may play a role in how resistant microbial communities are to perturbation, and therefore how community-level functions (like respiration) track changing environments. In our system, microbial communities at both ends of the gradient experience drought on a regular basis despite variation in MAP, potentially explaining why most bacterial taxa were not specialized to a specific rainfall regime. Similarly, in other soil reciprocal transplant experiments, fairly modest shifts in bacterial community composition were observed over 1 to 17 year study intervals when climate envelopes in home and away treatments were relatively similar (Lazzaro, Gauer and Zeyer 2011; Bond-Lamberty et al.2016). For example, reciprocal transplants of soil cores between forest stand types or soil types had no effect on community composition (Hannam, Quideau and Kishchuk 2007; Bond-Lamberty et al.2016). Bottomley et al. (2006) found that bacterial communities in meadow soils were largely unchanged following transplantation to a forest, although the reverse treatment did modify community structure. By contrast, moving soils between oak and grassland rapidly altered oak microbial communities that experienced hotter, drier conditions outside their historical climate (Waldrop and Firestone 2006). Microbial community composition also responded rapidly to transplantation along a dramatic (4.5°C) temperature gradient (Vanhala et al.2011). In summary, we suspect that more extreme climatic shifts (relative to natural climate variability) may be required to alter the relative abundance of generalist bacterial taxa. The compositional resistance of bacterial communities in both common gardens helps explain why patterns of respiration across soil types were unaltered by the transplant and dispersal treatments. Differences in community structure among common gardens, sites of soil origin and transplant treatments were driven by species turnover, not nestedness (local extinction), suggesting that changes in the identity or relative abundance of taxa drove the subtle compositional responses to the experimental treatments. Local specialization is a potential ecological mechanism for these compositional changes following transplantation. In support of this idea, OTUs with the highest ecological specificity tended to be more abundant in the garden matching their home precipitation region (drier west vs. wetter east), and were highest for the garden soils that were located at their specific site of origin. Moreover, decreases in the abundance of specialist OTUs following transplantation to a novel precipitation region (without changes in the chemical composition of their surrounding soil matrix) suggest that many of these taxa were responding directly to water availability. These results are consistent with laboratory studies on cultured microbes that have demonstrated ecological specialization related to moisture niche (Lennon et al.2012). Additionally, field studies have found strong niche conservatism in microbial responses to rapid change in soil moisture (Placella, Brodie and Firestone 2012). In our study, over 40% of ecological specialists were Actinobacteria, which are often found to be especially responsive to drought (Barnard, Osborne and Firestone 2013; Evans, Wallenstein and Burke 2014). This may explain why ecological specialization indices were higher in soils from drier western sites and in common garden (WFC and MOR) soils, where Actinobacteria comprised a relatively greater proportion of observed taxa. Despite the fact that most taxa were habitat generalists, specialists tended to be the most locally abundant OTUs. Specialization indices were calculated with respect to site of soil origin, not MAP alone, and so may also capture microbial responses to multiple non-climatic site factors. Soil characteristics are often reported as microbial drivers; however, soil organic C content, texture and pH explained only 2% of the variation in community structure. Although it is likely that unmeasured soil factors also shape bacterial communities, soil properties do not appear to be major controllers in the current system. Our results directly contrast with other studies that report major effects of pH and carbon availability on microbial community structure (Lauber et al.2008; Fierer et al.2009), perhaps because the relatively uniform geology of the Edwards Plateau constrains the amount of variation in these properties. We found no evidence that enhanced dispersal of microbes from the surrounding environment substantially affected bacterial community structure, similar to previous studies demonstrating a lack of dispersal limitation in soil bacteria (van der Gucht et al.2007). Cores subject to the enhanced dispersal treatment tended to have slightly more diverse communities based on number of taxa, but this response was not observed in all cases. Although these results are consistent with high microbial resistance, it is also possible that the microbial inoculation treatments failed due to environmental conditions: both common gardens experienced a historic drought at the start of the experiment, and bacteria in the soil inoculum may not have survived long under exceedingly hot, dry conditions. Alternatively, successful establishment of immigrating taxa—rather than immigration rate itself—may be the limiting factor for community-level shifts in response to environmental change (Fukami et al.2010). Establishment of new taxa may be impeded by biotic interactions with resident microbial species (Jones et al.2017), or by edaphic conditions unrelated to rainfall. It is possible that these issues were further exacerbated by low microbial numbers in our inoculum treatments; however, given typical microbial biomass values for these soils, each of the four repeated inoculum treatments likely delivered many thousands of viable cells. CONCLUSIONS Historical rainfall regime shapes bacterial community structure and function along a steep precipitation gradient (Hawkes et al.2017). We found that differences in community composition among sites are driven by species turnover rather than nestedness, consistent with some degree of ecological specialization based on precipitation regime. Although specialist taxa tended to exhibit more pronounced changes in abundance when transplanted to regions of higher or lower rainfall, compositional and functional responses to the transplant treatment were minimal. This suggests that the habitat generalists which dominated these microbial communities confer a high degree of ecological resistance to altered environmental conditions, with potentially important consequences for ecosystem-scale carbon cycling. SUPPLEMENTARY DATA Supplementary data are available at FEMSEC online. Acknowledgements We are grateful for assistance with fieldwork and labwork provided by Jim Dula, Nathan Christianson and Sarah Chen. Site access was kindly permitted by the Texas Parks and Wildlife Department, the Texas Ecolab Program and the Texas Historical Commission. Previous drafts of the manuscript were improved by comments from Jennifer Rocca and Stephanie Kivlin. FUNDING This material is based upon work supported by Texas Ecolab and the University of Texas funds to CVH. BGW was supported by a NSF Graduate Research Fellowship (DGE-1110007) and a NSF Doctoral Dissertation Improvement Grant (DEB-1210361). Conflict of interest. None declared. REFERENCES Barnard R, Osborne C, Firestone M. Responses of soil bacterial and fungal communities to extreme desiccation and rewetting. ISME J  2013; 7: 2229– 41. Google Scholar CrossRef Search ADS PubMed  Baselga A. Partitioning the turnover and nestedness components of beta diversity. Global Ecol Biogeogr  2009; 19: 134– 43. Google Scholar CrossRef Search ADS   Baselga A, Orme CDL. betapart: an R package for the study of beta diversity. Methods Ecol Evol  2012; 3: 808– 12. Google Scholar CrossRef Search ADS   Bond-Lamberty B, Bolton H, Fansler S et al.   Soil respiration and bacterial structure and function after 17 years of a reciprocal soil transplant experiment. PLoS One  2016; 11: e0150599. Google Scholar CrossRef Search ADS PubMed  Borcard D, Legendre P, Drapeau P. Partialling out the spatial component of ecological variation. Ecology  1992; 73: 1045– 55. Google Scholar CrossRef Search ADS   Bottomley PJ, Yarwood RR, Kageyama SA et al.   Responses of soil bacterial and fungal communities to reciprocal transfers of soil between adjacent coniferous forest and meadow vegetation in the Cascade Mountains of Oregon. Plant Soil  2006; 289: 35– 45. Google Scholar CrossRef Search ADS   Cadotte MW. Metacommunity influences on community richness at multiple spatial scales: a microcosm experiment. Ecology  2006; 87: 1008– 16. Google Scholar CrossRef Search ADS PubMed  Caporaso JG, Kuczynski J, Stombaugh J et al.   QIIME allows analysis of high-throughput community sequencing data. Nat Methods  2010; 7: 335– 6. Google Scholar CrossRef Search ADS PubMed  Cline LC, Zak DR. Dispersal limitation structures fungal community assembly in a long-term glacial chronosequence. Env Microbiol Rep  2013; 16: 1538– 48. Google Scholar CrossRef Search ADS   deVries FT, Shade A. Controls on soil microbial community stability under climate change. Front Microbiol  2013; 4: 265. Google Scholar PubMed  Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics  2010 26: 2460– 1. Google Scholar CrossRef Search ADS PubMed  Evans S, Wallenstein M, Burke I. Is bacterial moisture niche a good predictor of shifts in community composition under long-term drought? Ecology  2014; 95: 110– 22. Google Scholar CrossRef Search ADS PubMed  Fierer N, Strickland MS, Liptzin D et al.   Global patterns in belowground communities. Ecol Lett  2009; 12: 1238– 49. Google Scholar CrossRef Search ADS PubMed  Freedman Z, Zak DR. Soil bacterial communities are shaped by temporal and environmental filtering: evidence from a long-term chronosequence. Env Microbiol Rep  2015; 17: 3208– 18. Google Scholar CrossRef Search ADS   Fukami T, Dickie IA, Paula Wilkie J et al.   Assembly history dictates ecosystem functioning: evidence from wood decomposer communities. Ecol Lett  2010; 13: 675– 84. Google Scholar CrossRef Search ADS PubMed  Gianuca AT, Declerck SAJ, Lemmens P et al.   Effects of dispersal and environmental heterogeneity on the replacement and nestedness components of β-diversity. Ecology  2017; 98: 525– 33. Doi: /101002/ecy1666. Google Scholar CrossRef Search ADS PubMed  Hannam KD, Quideau SA, Kishchuk BE. The microbial communities of aspen and spruce forest floors are resistant to changes in litter inputs and microclimate. Appl Soil Ecol  2007; 35: 635– 47. Google Scholar CrossRef Search ADS   Hawkes CV, Waring BG, Rocca JD et al.   Historical climate controls soil respiration responses to current soil moisture. P Natl Acad Sci USA  2017; 114: 6322– 7. Google Scholar CrossRef Search ADS   Hawkes CV, Keitt TH. Resilience vs historical contingency in microbial responses to environmental change. Ecol Lett  2015; 18: 612– 25. Google Scholar CrossRef Search ADS PubMed  Hughes J, Hellmann J, Ricketts T et al.   Counting the uncountable: statistical approaches to estimating microbial diversity. Appl Environ Microb  2001; 67: 4399. Google Scholar CrossRef Search ADS   Jones ML, Ramoneda J, Rivett DW et al.   Biotic resistance shapes the influence of propagule pressure on invasion success in bacterial communities. Ecology  2017; 98: 1743– 9. Google Scholar CrossRef Search ADS PubMed  Kindt R, Coe R. Tree Diversity Analysis: A Manual and Software for Common Statistical Methods for Ecological and Biodiversity Studies . Nairobi, Kenya: World Agroforestry Centre, 2005. Kivlin SN, Hawkes CV. Temporal and spatial variation of soil bacteria richness, composition, and function in a Neotropical rainforest. PLoS One  2016; 11: e0159131. Available at: https://doi.org/10.1371/journal.pone.0159131. Google Scholar CrossRef Search ADS PubMed  Klindworth A, Pruesse E, Schweer T et al.   Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res  2013; 41. doi: 101093/nar/gks808. Lauber CL, Strickland MS, Bradford MA et al.   The influence of soil properties on the structure of bacterial and fungal communities across land-use types. Soil Biol Biochem  2008; 40: 2407– 15. Google Scholar CrossRef Search ADS   Lawrence D, Bell T, Barraclough TG. The effect of immigration on the adaptation of microbial communities to warming. Am Nat  2016; 187: 236– 48. Google Scholar CrossRef Search ADS PubMed  Lazzaro A, Gauer A, Zeyer J. Field-scale transplantation experiment to investigate structures of soil bacterial communities at pioneering sites. Appl Environ Microb  2011; 77: 8241– 8. Google Scholar CrossRef Search ADS   Legendre P, Borcard D, Peres-Neto PR. Analyzing beta diversity: partitioning the spatial variation of community composition data. Ecol Monogr  2005; 75: 435– 50. Google Scholar CrossRef Search ADS   Lennon J, Aanderud Z, Lehmkuhl B et al.   Mapping the niche space of soil microorganisms using taxonomy and traits. Ecology  2012; 93: 1867– 79. Google Scholar CrossRef Search ADS PubMed  Lindström ES, Östman Ö. The importance of dispersal for bacterial community composition and functioning. PLoS One  2011; 6: e25883. Google Scholar CrossRef Search ADS PubMed  Lumley T. leaps: regression subset selection, including exhaustive search. 2017. https://cranr-projectorg. Version: 3.0. Mariadassou M, Pichon S, Ebert D. Microbial ecosystems are dominated by specialist taxa. Ecol Lett  2015; 18: 974– 82. Google Scholar CrossRef Search ADS PubMed  McCune B, Mefford MJ. PC-ORD: Multivariate Analysis of Ecological Data . Gleneden Beach, Oregon, USA: MjM Software, 2011. Moyano FE, Manzoni S, Chenu C. Responses of soil heterotrophic respiration to moisture availability: an exploration of processes and models. Soil Biol Biochem  2013; 59: 72– 85. Google Scholar CrossRef Search ADS   Oksanen J, Blanchet FG, Friendly M et al.   Package ‘vegan.’ 2016. https://cranr-projectorg. Version: 2.4-1. Placella SA, Brodie EL, Firestone MK. Rainfall-induced carbon dioxide pulses result from sequential resuscitation of phylogenetically clustered microbial groups. P Nat Acad Sci USA  2012; 109: 10931– 6. Google Scholar CrossRef Search ADS   Quast C, Pruesse E, Yilmaz P et al.   The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res  2012; 41: D590– 6. Google Scholar CrossRef Search ADS PubMed  Ulrich W, Almeida-Neto M, Gotelli NJ. A consumer's guide to nestedness analysis. Oikos  2009; 118: 3– 17. Google Scholar CrossRef Search ADS   van der Gucht K, Cottenie K, Muylaert K et al.   The power of species sorting: local factors drive bacterial community composition over a wide range of spatial scales. P Natl Acad Sci USA  2007; 104: 20404– 9. Google Scholar CrossRef Search ADS   van Tienderen PH. Generalists, specialists, and the evolution of phenotypic plasticity in sympatric populations of distinct species. Evolution  1997; 51: 1272– 380. Google Scholar CrossRef Search ADS   Vanhala P, Karhu K, Tuomi M et al.   Transplantation of organic surface horizons of boreal soils into warmer regions alters microbiology but not the temperature sensitivity of decomposition. Glob Change Biol  2011; 17: 538– 50. Google Scholar CrossRef Search ADS   Waldrop MP, Firestone MK. Response of microbial community composition and function to soil climate change. Microb Ecol  2006; 52: 716– 24. Google Scholar CrossRef Search ADS PubMed  Wang Q, Garrity GM, Tiedje JM et al.   Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy. Appl Environ Microb  2007; 73: 5261– 6267. Google Scholar CrossRef Search ADS   © FEMS 2018. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com

Journal

FEMS Microbiology EcologyOxford University Press

Published: Feb 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 12 million articles from more than
10,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Unlimited reading

Read as many articles as you need. Full articles with original layout, charts and figures. Read online, from anywhere.

Stay up to date

Keep up with your field with Personalized Recommendations and Follow Journals to get automatic updates.

Organize your research

It’s easy to organize your research with our built-in tools.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

Monthly Plan

  • Read unlimited articles
  • Personalized recommendations
  • No expiration
  • Print 20 pages per month
  • 20% off on PDF purchases
  • Organize your research
  • Get updates on your journals and topic searches

$49/month

Start Free Trial

14-day Free Trial

Best Deal — 39% off

Annual Plan

  • All the features of the Professional Plan, but for 39% off!
  • Billed annually
  • No expiration
  • For the normal price of 10 articles elsewhere, you get one full year of unlimited access to articles.

$588

$360/year

billed annually
Start Free Trial

14-day Free Trial