Soil microbial community variation correlates most strongly with plant species identity, followed by soil chemistry, spatial location and plant genus

Soil ecologists have debated the relative importance of dispersal limitation ("everything" is not "everywhere") and ecological factors in determining the structure of soil microbial communities. The relative explanatory power of spatial and ecological factors, including plant species identity and even plant relatedness, for different fractions of the soil microbial community (i.e. bacterial and fungal communities) are poorly understood. We sampled field soils in a northern California field site, and find that soil microbial community variation correlates most strongly with plant species identity, followed by soil chemistry, spatial location and plant genus.


Introduction
The study of microbial community structure has long centred around the debate between the hypothesis that 'everything is everywhere' (Bass Becking 1934), and the alternative that dispersal limitation influences microbial community structure (O'Malley 2007). Because the soil microbial community is diverse, heterogeneous and difficult to characterize (Singh et al. 2004), understanding the environmental correlates of soil microbial communities has lagged behind community ecology studies in other systems (Fierer and Ladau 2012;. Thus the relative explanatory powers of ecological factors (e.g. plant species identity, soil chemistry and plant relatedness) and dispersal limitation (i.e. spatial location) for determining soil microbial communities are not well characterized (reviewed in Berg and Smalla 2009). By comparing bacterial and fungal community datasets, we explore whether ecological and spatial factors structure soil microbial communities, and, if so, how bacterial and fungal communities differ.
Plant -microbial relationships are often plant speciesspecific (reviews in Ehrenfeld et al. 2005;Hardoim et al. 2008;Berg and Smalla 2009). For example, legume species are often associated with particular strains or species of rhizobia bacteria (Long 2001). Plants exude chemicals from their roots that can foster beneficial microbes in the rhizosphere (reviewed in Ehrenfeld et al. 2005;Bais et al. 2006;Compant et al. 2010). For example, Arabidopsis thaliana exudes malic acid in the presence of a pathogen, and malic acid attracts a beneficial bacterium, Bacillus subtilis, protecting the roots from the pathogen (Rudrappa et al. 2008). However, plant species-specific relationships are not universally found for soil microbial communities (Ehrenfeld et al. 2005), and data sets comparing bacterial and fungal communities for plant species-specificity are rare (but see, e.g. Ushio et al. 2008).
Soil chemistry also influences soil microbial community composition, diversity and activity (Ehrenfeld et al. 2005). For example, differences in pH explained 70 % of the variance in bacterial community diversity across soils sampled from North and South America (Fierer and Jackson 2006). Greater soil fertility (NPK) increases bacterial biomass and activity in grassland mesocosms and these effects interact with plant species identity in some systems (Innes et al. 2004) and not in others (Bardgett et al. 1999), thus effects of soil chemistry are often system-dependent. Whether or how these soil chemistry effects may interact with plant species identity or plant relatedness is not generally known.
Spatial location could also influence soil microbial composition, especially if microbial taxa are dispersal limited, resulting in spatial autocorrelation of microbial communities (e.g. Ettema and Wardle 2002;Fierer and Ladau 2012). Early workers argued that microbes do not experience significant dispersal limitation, and thus 'everything is everywhere: but the environment selects' (Bass Becking 1934, reviewed in O'Malley 2007. However, more recent work has found strong spatial autocorrelation in many microbial communities, such that spatially close communities are more similar than expected by chance, and this similarity decays with distance (Ettema and Wardle 2002). Such a pattern of spatial autocorrelation suggests that dispersal limitation does influence patterns of distribution for microbes (O'Malley 2007;Rout and Callaway 2012). The relative importance of dispersal limitation and ecological factors is still an area of ongoing research (Berg and Smalla 2009).
If closely related plants are similar in root morphology (e.g. Comas and Eissenstat 2009;Valverde-Barrantes et al. 2014), root exudates, or other drivers of microbial associations, plant relatedness might also be expected to correlate with soil microbial community structure. Consistent with this hypothesis, plant clades, in some cases, have similar responses to soil biota (e.g. Brandt et al. 2009;Liu et al. 2012;Reinhart et al. 2012;Anacker et al. 2014). Closely related plants also often share characteristics of root morphology, such as the absence of many fine roots (Valverde-Barrantes et al. 2014), a trait which is correlated with higher mycorrhizal colonization (Manjunath and Habte 1991;Brundrett 2002). In addition, there is evidence, especially in Orchidaceae, of closely related plants sharing similar soil fungal associates (Jacquemyn et al. 2011;Martos et al. 2012) and close relatives sometimes respond similarly to arbuscular mycorrhizal fungi (e.g. Reinhart et al. 2012;Lugo et al. 2015). However, such relationships are by no means universal, and there are cases were plant relatedness does not correlate with effects of soil biota (e.g. Reinhart and Anacker 2014). Previous studies have often focused on whole-soil effects (i.e. plant-soil feedbacks, e.g. Burns and Strauss 2011) or have focused on a narrow fraction of the soil community (e.g. mycorrhizae, Reinhart et al. 2012;Reinhart and Anacker 2014); therefore, there is little comparative evidence about the relative influence of plant relatedness on different components of the soil microbiome (e.g. bacterial and fungal communities).
We quantified the relative contributions of plant species, soil chemistry, spatial location and plant relatedness on soil microbial communities at Bodega Bay Marine Reserve, CA, USA. (i) We predicted that plant species identity may influence soil microbial community composition, if species differ in factors, such as root exudates, which influence microbial community composition. (ii) Soil chemistry (e.g. N, pH) may correlate with microbial communities, for example, if soil microbe species differ in their nutrient preferences or soil chemistry niche axes. (iii) Spatial separation may influence microbial community composition across samples, if microbes are dispersal limited. (iv) Finally, plant congeners may have similar soil microbial communities, if close relatives have similar traits or habitat preferences, which influence soil microbial community composition.

Methods
To explore microbial community composition, we sampled field soils and used DNA fragment analysis to characterize differences among soil samples. For each of the 14 angiosperm species, we sampled soils from 6 unique collection locations in the field, at least 100 m apart, at Bodega Bay Marine Reserve, CA, USA [see Supporting Information]. These species are Cirsium occidentale, C. quercetorum, Fragaria chiloensis, F. vesca, Gilia capitata ssp. chamissonis, G. millefoliata, Plantago erecta, P. subnuda, Rumex crassus, R. occidentalis, Sanicula arctopoides, S. crassicaulis, Trifolium fucatum and T. gracilentum. Soil samples were collected from the root zone of each focal plant species to 15 cm depth, roots were removed and soils were mixed within a core. Soils therefore represent rhizosphere and non-rhizosphere soil within the active root zone of the selected plants species, which we expect to be under the influence of the plant and therefore represent microbial communities cultivated by the influence of the representative plant species. Soil samples were stored on dry ice in the field and moved to 280 8C within ,12 h. To quantify the role of spatial proximity, we mapped collection locations using GPS coordinates.

Soil bacterial and fungal communities
Terminal-restriction fragment length polymorphism (TRFLP) was used to characterize the soil bacterial and fungal communities following established methods (as in, for example, Burke et al. 2008). Terminal-restriction fragment length polymorphism is a DNA fragmentation technique that quantifies differences in the soil microbial community among samples but which cannot identify microbes to specific taxa. Terminal-restriction fragment length polymorphism estimates of microbial community structure have been found to be generally robust and concur with alternative approaches such as 454 pyrosequencing in documenting microbial community differences (Pilloni et al. 2012). See details of the TRFLP methodology in the Supporting Information. In brief, for each soil sample, DNA was extracted from two 500-mg subsamples of collected soil, amplified using PCR to target bacterial or fungal-specific target regions and restriction enzymes were used to digest the PCR product. We then averaged the subsample extractions within a soil sample for analysis. Terminal-restriction fragment length polymorphism results in size fragments that represent operational taxonomic units (OTUs) in our analyses.

Soil chemistry
To characterize soil chemical properties, soil samples taken from each soil collection (n ¼ 84) were analysed for NO 3 (ppm), Olsen P (ppm), K (ppm), K (meq/100 g), Na (ppm), Ca (meq/100 g), Mg (meq/100 g), cationexchange capacity, organic matter and pH by the Division of Agriculture and Natural Resources Analytical Laboratory at the University of California, Davis.

Phylogeny estimation
We sampled seven congeneric pairs and six soil replicates per species, and two subsamples per soil sample, to assess heterogeneity in soil microbial communities. We tested for effects of plant phylogeny, as well as tested just for effects of plant genus on microbial community composition (see below).
To estimate the phylogeny among our sampled plant species, we downloaded all available DNA regions for matK, ITS and trnL-trnF from genBank for the 14 sampled plant species [see Supporting Information]. Each species was represented in the data set by at least 1 DNA region, matK was present for 7 species, ITS for 13 and trnL-trnF for 11 [see Supporting Information]. We aligned the sequences separately for each DNA region with MEGA version 5.2.2 (Tamura et al. 2011) using MUSCLE (Edgar 2004a, b) with the default gap open penalty of 2400 and then concatenated the alignments. We conducted a maximum likelihood search using Garli version 0.951 (Zwickl 2006) using a GTR + G + I model of evolution. To incorporate prior knowledge about angiosperm phylogeny, we enforced a constraint tree with the known relationships among genera within the Asterids (Stevens 2009). Bootstrap analysis with 100 replicates was conducted without constraints. Thus bootstrap values reflect repeatability of clades, based on the raw DNA data. We calibrated the branch lengths using r8s with the NPRS method (Sanderson 2003) and dates fixed at the Fabaceae, Rosaceae, Asteraceae, Apiaceae, Polygonaceae, Lamiales, Caryophyllales, coreeudicots, eurosid1 and easterid2 based on published fossil calibration dates (Wikström et al. 2001).

Data analysis
Soil community diversity. We ordinated the community data using non-metric multidimensional scaling. We used the 'metaMDS' function in the 'vegan' package (Oksanen et al. 2010), using 20 random starts and a Bray-Curtis distance matrix (Faith et al. 1987). We conducted ordinations with and without data transformation, and with and without removing rare taxa (taxa present in ≤7 % of subsamples) prior to ordination. These strategies did not reduce ordination stress, and we therefore present ordinations on untransformed whole community data. In addition, we conducted visual inspection of a graph of stress versus the number of dimensions in the ordination, and chose the number of axes that reduced stress below 0.20, where the graph appeared to asymptote (k ¼ 4). NMDS ordination resulted in a reduced number of AoB PLANTS www.aobplants.oxfordjournals.org dimensions that describe differences in soil communities, and with which we could conduct further analyses.
To describe the amount of community variance potentially represented by the ordination, we conducted a post hoc partial Mantel correlation between a Euclidean distance matrix from the ordination and the total community data matrix using Bray -Curtis distances following recommendations in McCune and Grace (2002). We used multi.mantel in the phytools package with 1000 permutations (Revell 2012) to conduct this analysis and present R 2 values converted to percentages.
Do plant species identity, soil chemistry, spatial location and phylogenetic relatedness predict similarity in rhizosphere microbial communities? We used two statistical tools to explore the correlations in this data set: phylogenetic eigenvector regression (PVR) and a variance partitioning analysis. The PVR approach allows us to test for statistical significance of effects of our predictors on soil microbial community structure. Because many predictors may covary in this observational data set, the variance partitioning approach allows us to ask about the relative explanatory power of combinations of predictors and is purely descriptive.
First, we fit four PVR models, each with one of the ordination axes as the response variable; predictors in each model included the eigenvectors that represented plant species, soil chemistry, spatial location and plant phylogeny, as described below. To acquire phylogenetic eigenvectors, we first decomposed the phylogeny with 14 plant taxa using PVRdecomp in the PVR library (Diniz-Filho et al. 1998, 2012. We examined a graph of the eigenvalues against the vector number for the decomposition and used a visual break in the data to determine which eigenvectors to consider for the model selection, which included the first seven. Next, we decomposed spatial distance into a set of two eigenvectors (Bellier et al. 2007). To select a minimal model for each ordination axis, we fit a full linear model with the first seven phylogenetic eigenvectors and two spatial eigenvectors using the lm function. We then choose the subset of phylogenetic eigenvectors that were statistically significant (P , 0.05) for the reduced model (Diniz-Filho et al. 2012). This approach never resulted in more than three phylogenetic predictors. Because identical values for eigenvectors were associated with each species, species identity could not simultaneously be included in this model. We compared this full model (e.g. fungal community axis 1 plant phylogeny + spatial location) to models dropping effects from the model using a likelihood ratio test to determine the statistical significance of each predictor. Because 'phylogeny' and 'spatial location' consisted of multiple vectors, this approach tested the overall significance of these predictors of microbial community structure. To select minimal soil chemical models, we used stepAIC on a model with all soil chemical vectors and included only statistically significant predictors in the final model. This reduced variance inflation factors to below 2.5 for all predictors. Because phylogeny co-varied with soil chemistry (Spearman rank correlation range: 20.44 to 0.59) and models with both phylogeny and soil chemistry resulted in high multicollinearity (VIF . 100), phylogeny and soil chemistry effects were tested separately. Finally, we report the adjusted R 2 for models with each of four classes of predictors (i.e. plant species, soil chemistry, spatial location and plant phylogeny) individually, to estimate the maximum relative explanatory power of each, adjusted for the number of predictors in the model. Because our sampling was concentrated on congeners, we also conducted a PVR analysis with plant genus as a predictor, rather than plant phylogeny, to ask whether plant genus alone could predict differences in microbial community composition.
Standard diagnostic plots were examined and residuals were normally distributed for all PVR models. We used adjusted R 2 value comparisons among models to ensure that models were not over-fit. Variance inflation factors were ,2.5 for all predictors in all final models (median VIF ¼ 1.03), which is generally considered very low multicollinearity (Fox and Monette 1992;O'Brien 2007).
Finally, some of these predictors may covary, such as if closely related species prefer the same habitats and are thus spatially close together, and multiple predictors combined could influence microbial community composition. To describe the combined effects of plant phylogeny, soil chemistry and spatial location on soil community structure, we used variance partitioning methods (Desdevises et al. 2003;Peres-Neto et al. 2006;Gonçalves-Souza et al. 2014). See the Supporting Information for more details about this analysis. All analyses were conducted in the R statistical package (version 3.0.2) (R Development Core Team 2008). Data are available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.9bf7f.

Soil chemistry
Soils from the Bodega Bay site were primarily loamy sand with a somewhat acidic average pH of 6.07 + 0.05 SE (range ¼ 5.01 -8.03). Organic matter was generally low with a mean 4.2 % + 0.2 SE (range ¼ 0.6-15.5 %). Salinity was highly variable and ranged from 12 to 390 ppm with a mean of 91.8 + 6.7 SE at this coastal site.

Bacterial communities
NMDS ordination resulted in a stress of 0.12 with four ordination axes ( Fig. 1A and B), and the total ordination (axes 1-4 combined) explained 82 % of the bacterial community composition in a partial Mantel test. Bacterial communities in the rhizosphere exhibited plant speciesspecific patterns, and plant species explained 15, 15 and 38 % of the variance on ordination axes 2, 3 and 4, respectively (Table 1). Variation in pH explained 19 % of the variance in bacterial community composition on axis 2 (Table 1). Variation in Mg and Ca across samples explained 16 % of the variation on axis 3 (Table 1). Variation in K, Na and Ca explained 11 % of the variance on axis 4 (Table 1). Spatial location explained 19 and 16 % of the variation on ordination axes 3 and 4 (Table 1). Plant phylogeny was a significant but poor predictor of bacterial community composition, explaining 1, 5, 9 and 8 % of the variation on ordination axes 1-4 (Table 1). An alternative analysis approach with plant genus instead of plant phylogeny explained 12 % of the variance on ordination axis 3 [see Supporting Information- Table S2], and this pattern was primarily driven by the Gillia and Sanicula congeners, which were very similar in their bacterial community composition on axis 3 (Fig. 2). In the variance partitioning analysis, soil chemistry alone explained more of the variance in bacterial community composition than other factors, and combinations of predictors (e.g. spatial location, soil chemistry and plant phylogeny combined), explained additional variation [see Supporting Information- Fig. S2].

Fungal communities
The NMDS ordination for the fungal data set had a stress of 0.15 ( Fig. 1C and D), and the total ordination (axes 1-4 combined) explained 71 % of the fungal community composition in a partial Mantel test. Plant species identity explained 56, 31 and 20 % of the variation on ordination axes 1, 2 and 3, respectively (Table 1). Phosphorus, Ca and Mg explained 45 % of the variance on fungal axis 1 (Table 1). Variation among soil samples in pH explained  (Table 1). Spatial location also explained a significant amount of the variance in fungal community composition, explaining 43, 9 and 10 % of the variance on ordination axes 1, 3 and 4 (Table 1). Plant phylogeny explained 29 % of the variance on ordination axis 2 (Table 1). The alternative analysis with plant genus, rather than plant phylogeny, suggested that plant genus explained 16 and 30 % of the variance on ordination axes 1 and 2, respectively [see Supporting Information- Table S2]. For example, Rumex species fungal rhizosphere communities were very close in ordination space (Fig. 1C, blue points) as were Gillia congeners (Fig. 1C, grey points), and Sanicula congeners (Fig. 1C, pink points). This signal of plant genus was especially strong on ordination axis 2 (Fig. 2, Supporting Information- Table S2). Combinations of the predictors explained larger percentages of the fungal community composition than individual predictors alone [see Supporting Information- Fig. S2].  2. Soil bacterial communities were significantly, but weakly, correlated with plant phylogeny on ordination axis 3 (Table 1). Plant phylogeny explained 29 % of the variance in soil fungal communities on ordination axis 2 (Table 1).

Discussion
Plant species identity was the single best statistical predictor of soil microbial community composition Consistent with a large body of literature (reviews in Brundrett 2002;Ehrenfeld et al. 2005; Berg and Smalla 2009), we found a signal of plant species identity on the soil bacterial and fungal communities, and plant species was the single best statistical predictor of soil microbial community variation (Table 1). Variance in bacterial community composition was higher than for fungal communities, and no combination of predictors could explain a large amount of this variance [see Supporting Information- Fig. S2]. This large variance could reflect a high degree of fine-scale heterogeneity within the root zone to which bacterial communities respond. For example, the presence of fungal biomass was found to influence the community structure of important bacterial functional groups (e.g. denitrifying bacterial) at fine spatial scales within a forest system . This high heterogeneity of microsites for bacterial within the root zone, coupled with the high degree of species diversity within the bacterial domain, could be responsible for the weaker predictors of bacterial communities in our study. In addition, our examination of bacterial community structure using general 16S rDNA primers could miss the potential relationship between plant species and specific bacterial functional groups. Analysis of the nitrogen-fixing functional group, for example, could reveal plant effects on that functional group, because we expect some plants to have co-evolved relationships with these nitrogen-fixing bacteria (e.g. rhizobia-plant relationships; Long 2001). The weaker signal of plant species identity on bacteria could therefore be due to our focus on overall bacterial community structure and not functional group distribution, which could yield different results. Plant species explained a significant amount of the variance in fungal community composition ( Table 1). Effects of plant species identity on soil microbial communities might be the result of 'active' selection by plants for soil microbes (e.g. via root exudates), a 'passive' by-product of nutrient uptake, root structure (Singh et al. 2004) or shared habitat preferences by plants and microbes. It is increasingly accepted that soil fungi, and especially mutualistic organisms such as mycorrhizal fungi (Brundrett 2002), may form plant species-specific relationships that can affect plant growth and community composition. Plant speciesspecific associations with soil microbes likely contribute to the strong plant-soil feedbacks observed in many systems (Kulmatiski et al. 2008), which have important implications for coexistence and community assembly (Bever et al. 2010).

Soil chemistry was correlated with soil microbial community composition
Soil chemistry had the second strongest correlation with differences in microbial communities (Table 1). Previous studies have found strong correlations between soil chemistry and soil microbial communities at large spatial scales (e.g. Lauber et al. 2009), consistent with our results at smaller spatial scales. Variation in pH is often correlated with bacterial communities at the continent scale (Lauber et al. 2009), and we found pH to correlate with bacterial community composition. Fungal communities, including AMF and EMF communities, also frequently differ across soils with differences in pH (Kluber et al. 2012), also consistent with the correlation between fungal ordination axis 2 and pH. Further, many fungi can affect C and N availability and cycling through their saprotrophic capabilities (e.g. production of extracellular enzymes for liberation of organic nutrients) (Leake et al. 2002), though we found no strong evidence for correlations between the fungal community and N.
Spatial location was the third best predictor of microbial community composition and was correlated with soil chemistry and plant relatedness We found significant associations between soil microbial communities and spatial location (Table 1). Thus 'everything' is not 'everywhere' in this data set, as many other workers have also found for soil microbial communities (reviews in Ettema and Wardle 2002;Rout and Callaway 2012). This is consistent with the hypothesis that dispersal limitation influences microbial distributions even at fairly restricted spatial scales [see Supporting Information- Fig. S1]. The correlations we observed between soil chemistry, plant relatedness and spatial locations on fungal community composition suggest that these factors co-varied in the field, making disentangling the relative importance of these predictors difficult [see Supporting Information- Fig. S2].
Plant congeners had similar soil fungal, but not bacterial, communities Plant relatedness explained less of the variation in soil microbial communities than plant species, soil chemistry or spatial location (Table 1). Plant congeners had similar soil fungal (but not bacterial) communities, even when they occurred in divergent habitats. For example, congeners in Sanicula and Rumex had similar fungal communities ( Fig. 2) but they are found in different habitats and are spatially distant (.1 km) from one another [see Supporting Information- Fig. S1]. Thus, spatial location and similarity in habitat preferences between close relatives are not the sole drivers of the correlation between plant AoB PLANTS www.aobplants.oxfordjournals.org relatedness and soil fungal communities. Similar root exudates or root morphology between close relatives (Comas and Eissenstat 2009;Valverde-Barrantes et al. 2014) may be other explanations for shared fungal community composition.
DNA fragmentation tools, like TRFLP, do not identify taxa without additional work. Because our soil sampling procedure may have broken some plant roots, endophytic fungi could be contributing to the pattern we observed. Our general ITS fungal primers preferentially amplify nonarbuscular mycorrhizal fungi (D. J. Burke, pers. obs.;Taylor et al. 2014), suggesting that arbuscular mycorrhizal fungi are not driving the correlation that we observed between plant phylogeny and fungal rhizosphere community. Given the species-specific nature of some mycorrhizal associations (Brundrett 2002), we may be underestimating the importance of host plant relatedness in influencing fungal rhizosphere communities.
Meta-analyses suggest that plant relatedness does not explain much of the overall variance in plant-soil feedbacks, where plant performance is compared between conspecific soils and heterospecific soils (Mehrabi and Tuck 2015); however, phylogenetic relatedness in plant-soil feedbacks did improve predictions of plant abundance for 57 species within an old field community (Anacker et al. 2014). Some fractions of the soil community may be similar between closely related plants. For example, arbuscular mycorrhizal effects are similar between closely related plants (Reinhart et al. 2012;Lugo et al. 2015) and some mycorrhizal associates are specialized to particular plant clades (Brundrett 2002). Our results suggest that plant relatedness correlates with the fungal fraction of the soil community, but only weakly with the bacterial community. Thus the relative importance of fungal and bacterial components of the soil microbiota to plant performance may influence whether close relatives have similar plant -soil feedback effects on one another (Burns and Strauss 2011).

Conclusions
Plant species identity was the best predictor of bacterial and fungal community composition in the soil, followed by soil chemistry, spatial location and plant relatedness. Because 'everything' is not 'everywhere' in microbial communities, spatial location can confound our attempts to predict microbial community structure. Future work should use experimental manipulations to disentangle covarying predictors, such as soil chemistry, spatial location and plant relatedness in the field. Understanding the feedbacks between soil microbial community structure and plant communities may benefit from considering the evolutionary relationships among plants.

Sources of Funding
Our work was funded by the National Science Foundation (United States) DEB (11-20387) and by the Holden Arboretum Trust and the Corning Institution for Education and Research.

Supporting Information
The following additional information is available in the online version of this article - Table S1. Plant species sampled to examine rhizosphere soil microbial communities at Bodega Bay, CA, USA. Figure S1. Map of the spatial sampling locations for the 14 species and 6 replicate soil collections per species at Bodega Bay, CA, USA. Figure S2. Variance partitioning analysis on the bacterial and fungal data sets. Table S2. Linear PVR models were used to test the effects of plant species, spatial location and plant genus on soil microbial community composition.
Phylogenetic conservatism in plant-soil feedback and its implications for plant abundance. Ecology Letters 17:1613 -1621.