Anthropogenic land use shapes the composition and phylogenetic structure of soil arbuscular mycorrhizal fungal communities

Abstract

Arbuscular mycorrhizal (AM) fungi play an important role in ecosystems, but little is known about how soil AM fungal community composition varies in relation to habitat type and land-use intensity. We molecularly characterized AM fungal communities in soil samples (n = 88) from structurally open (permanent grassland, intensive and sustainable agriculture) and forested habitats (primeval forest and spruce plantation). The habitats harboured significantly different AM fungal communities, and there was a broad difference in fungal community composition between forested and open habitats, the latter being characterized by higher average AM fungal richness. Within both open and forest habitats, intensive land use significantly influenced community composition. There was a broad difference in the phylogenetic structure of AM fungal communities between mechanically disturbed and nondisturbed habitats. Taxa from Glomeraceae served as indicator species for the nondisturbed habitats, while taxa from Archaeosporaceae, Claroideoglomeraceae and Diversisporaceae were indicators for the disturbed habitats. The distribution of these indicator taxa among habitat types in the MaarjAM global database of AM fungal diversity was in accordance with their local indicator status.

By analysing AM fungal DNA in soil from different habitats, its shown that forested and open habitats harbour different fungal communities, while disturbance changes the phylogenetic structure of fungal communities.

AMF, 454 pyrosequencing, habitat types, soil biodiversity, disturbance

Introduction

Arbuscular mycorrhizal (AM) fungi (phylum Glomeromycota) colonize the roots of most terrestrial plants, facilitating mineral nutrient uptake from soil in exchange for plant-assimilated carbon (Smith & Read, 2008). Functional aspects of plant–AM fungal interactions have been a major focus of research during recent decades, and it has been shown that the taxonomic identity of AM fungal taxa can influence plant growth and performance (van der Heijden et al., 1998; Klironomos, 2003). Plant growth has also been found to respond differently to communities of AM fungi originating from different ecosystems (Johnson, 1993; Corkidi et al., 2002; Moora et al., 2004; Johnson et al., 2008; Uibopuu et al., 2009, 2012; Martinez & Johnson, 2010; Williams et al., 2011). Reflecting the important role played by local AM fungal communities in determining plant growth, there is now increasing interest in describing and explaining the distribution of AM fungal diversity in human-dominated landscapes (van der Heijden & Scheublin, 2007).

Global meta-analyses have shown that AM fungal community diversity and composition differ among broad habitat types such as forest, grassland and arable fields (Öpik et al., 2006; Kivlin et al., 2011). Several studies have also recorded variation of AM fungal community composition among natural vegetation types within the same region (Dumbrell et al., 2010; Lekberg et al., 2011; Meadow & Zabinski, 2012). However, most attention has been focussed on the effect of land-use intensity on AM fungal communities and specifically those found in agricultural ecosystems. High land-use intensity has been linked with low AM fungal diversity, based on studies of fungal spores (Boddington & Dodd, 2000; Oehl et al., 2003, 2004, 2010; Bainard et al., 2012) and the molecular diversity of AM fungi in soil (Lumini et al., 2010, 2011; Verbruggen et al., 2012), although no changes in diversity have been observed in some cases (Jansa et al., 2002; Mathimaran et al., 2007). AM fungal community composition in soil has also been shown to change along gradients of land-use intensity (studies of spores: Jansa et al., 2002; Oehl et al., 2004, 2010; using molecular techniques: Lumini et al., 2010; Miras-Avalos et al., 2011). Similarly, AM fungal communities associating with plant roots tend to exhibit low diversity in agricultural ecosystems with high land-use intensity (Helgason et al., 1998; Daniell et al., 2001; Hijri et al., 2006; Alguacil et al., 2008; Lumini et al., 2011; Schnoor et al., 2011; Bainard et al., 2012; although see Galvan et al., 2009) and changes in community composition along gradients of land-use intensity (Jansa et al., 2003; Alguacil et al., 2008; Li et al., 2010; Miras-Avalos et al., 2011; Schnoor et al., 2011; Bainard et al., 2012). At the same time, we lack a good understanding of the impact of land use beyond agricultural ecosystems. No comparative analyses have been conducted to determine how AM fungal communities vary among ecosystems with divergent land-use type and intensity in conditions where climate and geology do not vary.

Variation of AM fungal assemblages in nature is typically investigated using taxon-based approaches that do not consider the phylogenetic relationships of the taxa involved (e.g. studies cited above). Functional and life-history traits of AM fungi are considered to be phylogenetically conserved (Hart & Reader, 2002; Powell et al., 2009) meaning that phylogeny represents a viable proxy of life-history strategy in AM fungi (Chagnon et al., 2013). Variation in the phylogenetic composition of AM fungal communities could therefore provide an interesting insight into the question of whether the functional characteristics of AM fungal communities are associated with specific habitats or land-use types. Variation in the phylogenetic structure of AM fungal communities has been addressed in experimental systems (Maherali & Klironomos, 2007, 2012; Bainard et al., 2014) or single ecosystems (Busby et al., 2013; Lekberg et al., 2013; Moebius-Clune et al., 2013; Horn et al., 2014; Saks et al., 2014), but, as far as we are aware, it has not been studied across habitats and land-use gradients. Here, we focus on phylogenetic beta diversity – a measure of phylogenetic distance between communities (Graham & Fine, 2008), which allows community structure and the associated traits of species in a community to be evaluated in relation to differences in land-use type and intensity.

In this study, we addressed the soil AM fungal communities inhabiting an extensive gradient of habitat types and land-use intensities, from conventional arable fields, where large quantities of fertilizers and pesticides are used, through sustainable arable fields, permanent grasslands and commercial forest plantations, up to primeval forest under strict conservation. All habitats were located within a relatively small geographic area and experienced similar climatic and edaphic conditions. In addition to comparing specific habitats, we compared broad habitat types. First, we considered the distinction between forested and open habitats because earlier meta-analyses have identified differences in the AM fungal communities inhabiting them (Öpik et al., 2006; Kivlin et al., 2011). Second, we considered mechanically disturbed (regularly ploughed arable fields and spruce plantation where soil is disturbed during repetitive thinning and timber removal) and nondisturbed habitats, because mechanical disturbance may be responsible for changes in AM fungal community composition (Miras-Avalos et al., 2011; Schnoor et al., 2011). We identified soil AM fungal communities on the basis of 454 sequencing of fungal DNA, which has the benefit of incorporating both the extraradical hyphae and spores in soil (Hempel et al., 2007; Davison et al., 2012). We hypothesized that AM fungal diversity in soil will decrease with increasing land-use intensity. We also expected AM fungal community composition and phylogenetic structure to reflect broad ecosystem distinctions, particularly that of structurally open (grassland and agricultural) vs. forested habitats. We expect such information to show how human activity has changed current AM fungal communities.

Materials and methods

Study sites

The study was conducted in Estonia, in a variety of habitats, representing on one hand a distinction between forested and structurally open habitat (including agricultural habitats) and on the other a gradient of land use intensity. The following habitats were sampled, in order of most to least disturbed: open habitats – conventional intensively managed arable land (intensive, three sites); sustainable arable land (sustainable, two sites); permanent urban grassland (grassland, two sites); forested habitats – managed forest – 40-year-old spruce plantations on former arable land (plantation, two sites); and primeval forest – old-growth herb rich mixed spruce forest (forest, three sites) in a protected area (see Supporting Information, Table S1 and Table 1 for more information). AM fungal samples from primeval forest are published in Davison et al. (2012). All ecosystems were situated on a Gleyic luvisol soil with slightly calcaric moraine parent material (Reintam et al., 1987). Soil conditions were similar among study sites, except that the topsoil was more acid under coniferous tree canopies (Table 1).

Table 1

Location and soil characteristics of the study sites and details of sampling; number of samples (i.e. samples yielding < 9 hits), number of Glomeromycota sequences and number of AM fungal molecular taxa (VT) per study site

Habitat	Geographic coordinates		Sampling date			pHKCl				P (mg kg−1)					K (mg kg−1)						Ca (mg kg−1)							Mg (mg kg−1)								N (%)									Org. C (%)										Number of successful samples											Number of Glomeromycota sequences												Number of VT
Open
Intensive 1	58°15′03″N 26°17′46″E		06/2009			6.9				89					229						2303							362								0.2									2.1										9											2551												21
Intensive 2	58°21′45″N 26°31′06″E		08/2009			6.2				114					220						1445							135								0.12									1.48										9											3650												26
Intensive 3	58°12′57″N 27°14′44″E		08/2009			6.4				104					153						1654							176								0.13									1.55										9											3027												18
Sustainable 1	58°14′39″N 26°13′31″E		08/2009			6.3				97					83						1420							120								0.13									1.64										1											2081												16
Sustainable 2	58°15′19″N 26°21′36″E		08/2009			7.2				46					145						2682							184								0.14									1.68										9											8503												28
Grassland 1	58°22′N 26°43′E		09/2008			7.1				134					338						6031							341								0.3									4.08										7											2047												20
Grassland 2	58°22′37″N 26°43′31″E		09/2009			7.1				115					278						6688							321								0.5									6.01										9											13815												35
Forested
Plantation 1	58°14′59″N 27°17′44″E		09/2009			4.1				37					75						984							146								0.23									3.66										5											1164												9
Plantation 2	58°14′59″N 27°17′50″E		09/2009			4.2				48					68						969							181								0.24									3.6										4											276												5
Forest 1	58°16′48″N 27°19′28″E		05/2009			5.2				13					113						2289							242								0.42									6.46										8											2277												21
Forest 2	58°16′48″N 27°19′28″E		06/2009			4.4				15					47						1205							146								0.26									4.69										9											3743												20
Forest 3	58°16′48″N 27°19′28″E		09/2009			4.2				27					142						2340							210								0.678									12.69										9											1928												22

Sites under regular soil disturbance.

Sample collection

Each habitat type was represented by two or three sites. At each site, nine soil samples were collected from nine points on a regularly spaced 10 × 10 m sampling grid in either 2008 or 2009. Although sampling of different sites was not synchronous (Table 1), there were no systematic differences in timing between habitats and earlier analysis of the forest samples indicated that the timing of sampling did not influence AM fungal soil community composition (Davison et al., 2012). Each soil sample consisted of 10 g of soil collected from the 2- to 10-cm topsoil where most of the roots of herbaceous plant species are located. Soil samples were dried with silica gel and stored airtight at room temperature before further analyses. We also collected topsoil (2–10 cm) samples for chemical analysis (see Data S1) from each sampling location and pooled these within each site (altogether about 500 g per site).

Molecular analyses

DNA was extracted from dried soil. Glomeromycota sequences were amplified from soil DNA extracts using the SSU rRNA gene primers NS31 and AML2 (Simon et al., 1992; Lee et al., 2008), linked to 454-sequencing adaptors A and B, respectively, and following the 454-sequencing approach of Davison et al. (2012) and Öpik et al. (2013). This gene provided us with a larger comparative sequence dataset than would be available for any other genomic region (Öpik et al., 2010, 2014). For more details on molecular analysis, see Data S1.

Sequence reads were included in subsequent analyses only if they carried the correct bar code and forward primer sequences and were ≥ 170 bp long (excluding the bar code and primer sequence). Potential chimeras were detected and removed from the data using uchime (Edgar et al., 2011) in reference database mode (MaarjAM, Öpik et al., 2010) and the default settings. After stripping the bar code and primer sequences, we used the MaarjAM database of published Glomeromycota SSU rRNA gene sequences for taxonomic assignment of the obtained reads (status 31.03.2013). The MaarjAM database contains representative sequences covering the NS31/AML2 amplicon from published environmental Glomeromycota sequence groups and known taxa. As of April 2013, it contained a total of 6064 SSU rRNA gene sequence records that had been classified on the basis of phylogenetic analysis into manually curated OTUs with sequence similarity threshold ≥ 97% (frequently 99%) or so-called virtual taxa (VT cf. Öpik et al., 2009). VT nomenclature allows easy comparison of data from individual datasets and consistent taxon naming in time. More details about VT principles as implemented in the MaarjAM database can be found in Öpik et al. (2014). Such preclustered identified sequences permit blast-based taxonomy assignment of sequences without the need to use OTU clustering algorithms with set similarity thresholds and correct for inconsistencies in the taxonomic identification of INSD accessions. The taxonomic nomenclature used in MaarjAM follows recent changes in Glomeromycota taxonomy where possible (discussed in Öpik et al., 2013). For more details on bioinformatics, see Data S1.

A set of representative sequence reads has been deposited in the EMBL nucleotide collection (accession numbers HF954560–HF954914; the set consists of 1–2 randomly selected reads representing each VT from each site, dependent on availability).

Statistical methods

Sampling efficacy was assessed with rarefaction analysis of data subsets, using the specaccum() function from r package vegan (Oksanen et al., 2013). Linear mixed-effects models (LME) (Pinheiro et al., 2013) were used to test for differences in AM fungal VT richness in different habitats. As dependent variables, we used both observed richness and a rarefied estimate of richness (rarefied richness) to correct for differences in sampling intensity. We used rarefaction to the median number of sequence reads per sample, as analysis using analogous microbial datasets has shown this to be an optimal approach (de Cárcer et al., 2011). The LME models also contained habitat type as a fixed predictor and site as a random effect (random intercept).

For analyses of AM fungal community structure, we used quantitative data, where the proportions of different VT reads were used as a proxy for the relative abundance of AM fungal taxa in a sample (cf. Öpik et al., 2009 for methodological discussion). Similarity between communities was measured using Bray–Curtis dissimilarity. Variation in soil AM fungal community composition was visualized using nonmetric multidimensional scaling [three-dimensional nmds with 50 iterations; implemented using function metaMDS() from vegan]. Nested permanova [function nested.npmanova() from BiodiversityR (Kindt & Coe, 2005)] was used to make comparisons between the AM fungal communities present in different habitats, broad habitat types (open vs. forested) and disturbance regimes (mechanically disturbed vs. nondisturbed). Differences between habitat types in multivariate dispersion (Anderson, 2006) were tested using an LME model. The LME model included as a dependent variable the Euclidean distance between each sample and the respective multivariate group (habitat) centroid. The model also contained habitat type as a fixed predictor and site as a random effect (random intercept).

To identify the AM fungal taxa associated with particular land-use regimes and broad ecosystem distinctions (structurally open vs. forested, mechanically disturbed vs. nondisturbed), we used indicator species analysis (Dufrene & Legendre, 1997), as implemented by the function indval() from the r package labdsv (Roberts, 2012) and considered only those VT with an indicator value of at least 0.25. To account for the nested nature of the data, we tested the significance of indicator taxa using permutation (n = 999) of entire sites between habitat types. We also studied available information about the global distribution of indicator taxa with the help of the MaarjAM database (Öpik et al., 2010). For VT that were identified as significant indicators in this dataset, we calculated the proportion of records in MaarjAM derived from different habitat categories defined in that database (anthropogenic, forest, grassland and successional), where sufficient records were available (9 VT with < 10 records were not considered). These empirical proportions were then compared to proportions calculated under a null model constructed by permuting (n = 999) the habitat variable. Empirical proportions that were more extreme than the 2.5% or 97.5% quantiles of the null distribution were taken as evidence of significantly nonrandom VT-habitat association.

Mantel tests were used to detect correlations between community composition and soil variables. As a single measure of each soil variable was available for each site, AM fungal community composition was pooled for each site: The mean proportional abundance of VT at each site was calculated. A Mantel procedure was then carried out on pairs of distance matrices comprising the AM fungal community data (using Bray–Curtis dissimilarity) and one of the soil variables (using Euclidean distance).

Phylogenetic analysis

A phylogenetic tree containing the type sequences of all VT in the MaarjAM database (Öpik et al., 2013) was used to depict phylogenetic relationships between VT recorded in this study. Phylogenetic similarity between communities of AM fungi and their correlation with soil edaphic variables were investigated using nmds, permanova, a multivariate dispersion model and Mantel tests as described above. However, in place of a distance matrix based on taxonomy, these calculations used pairwise between-sample phylogenetic distance, estimated from the abundance-weighted mean pairwise phylogenetic distance (mpd) between taxa present in different communities (using function comdist() from r package picante (Kembel et al., 2010).

Results

A total of 45 061 454 reads carried a correct bar code, the correct NS31 primer sequence, were ≥ 170 nucleotides in length and received a hit (similarity ≥ 97%) against a virtual taxon (VT; i.e. a molecular operational taxonomic unit) from the MaarjAM database of published Glomeromycota SSU rRNA gene sequences (Öpik et al., 2010). Samples yielding < 9 hits and singleton VT were removed, leaving a data matrix consisting of 88 samples (i.e. 82% of collected samples) and 72 Glomeromycota VT (Table S2). In particular, a number of spruce plantation and sustainable arable land samples did not yield amplified Glomeromycota DNA, resulting in low replication from these habitats.

Richness of AM fungal communities

Rarefaction suggested that some undetected taxa remained in all habitat types, with sustainable arable land particularly lacking an asymptotic rarefaction curve (Fig. 1).

Fig. 1

Rarefaction analysis of soil AM fungal samples (VT) from different habitats.

The mean richness of AM fungal taxa per sample differed significantly between habitats (observed richness F_5,7 = 19.1, P < 0.001, rarefied richness F_5,7 = 26.4, P < 0.001). Observed richness in spruce plantation was significantly lower than richness in sustainable agriculture and grassland (Tukey, P < 0.05, Fig. 2a), while rarefied richness was significantly lower in spruce plantation than in all other habitats (Tukey, P < 0.05, Fig. 2b).

Fig. 2

Observed (a) and rarefied (b) arbuscular mycorrhizal fungal (AMF) taxon richness (VT) per soil sample in different habitats. Solid lines indicate medians; boxes and whiskers indicate quartiles and ranges, respectively; points indicate outlying values. The mean richness of AM fungal taxa differed significantly between habitats (observed richness F_5,7 = 19.1, P < 0.001, rarefied richness F_5,7 = 26.4, P < 0.001). Letters indicate significant differences between habitats (P < 0.05 Tukey's HSD test).

There was a strong positive correlation between VT richness per site and topsoil pH (R = 0.7, P = 0.019) and soil Ca content (R = 0.6, P = 0.02) across the sampled sites. Significant correlation between other measured soil properties (P, N, K, Mg and organic C) and VT richness was not detected.

AM fungal community composition

permanova and NMSD ordination (Fig. 3a) revealed significantly different taxonomic composition in AM fungal communities from different habitats (permanova pseudo-F = 3.5, P = 0.001, R² = 0.26), from structurally open and forested habitats (permanova pseudo-F = 3.8, P = 0.001, R² = 0.11) and from mechanically disturbed and nondisturbed habitats (permanova pseudo-F = 3.5, P = 0.01, R² = 0.09). Fungal communities from agricultural land were located separately on the ordination biplot, and the sustainable and intensive agricultural samples formed separate groups within this cluster. Forested communities also formed a distinct cluster, within which primeval forest differed from spruce plantation. Samples from permanent grassland were located between arable land and forest communities (Fig. 3a). The multivariate dispersion of samples differed significantly between groups: Mean dispersion was significantly greater among intensive agricultural samples (0.60) than among plantation (0.26) or sustainable agricultural samples (0.42); meanwhile, dispersion was significantly lower in plantation samples than in all other habitat types besides sustainable agriculture (0.53–0.60; Tukey's HSD test; Fig. 3a).

Fig. 3

nmds plots of variation in soil AM fungal communities (a) taxonomic composition and (b) phylogenetic composition. The taxonomic analysis was based on the Bray–Curtis dissimilarity between samples (stress = 0.18); the phylogenetic composition analysis used pairwise abundance-weighted phylogenetic distances (stress = 0.16). Ellipses indicate one SD around group centroids [forested and open habitats in (a); undisturbed and disturbed habitats in (b)].

Primeval forest and permanent grassland contained several indicator taxa that were either absent or sparse in other habitats: Indicator species analysis identified seven and six significant indicator taxa for these habitat types, respectively (all Glomeraceae, Table 2, indicator value > 0.25). Sustainable arable land (two Diversisporaceae taxa), intensive arable land (one Archaeosporaceae taxon) and spruce plantation (one Diversisporaceae and one Claroideoglomeraceae taxon) were characterized by fewer indicator taxa (Table 2, indicator value > 0.25). We also identified indicator taxa for broad ecosystem types. Structurally open ecosystems were characterized by nine and forested ecosystems by five indicator taxa (Table 2, indicator value > 0.25). Altogether, the analysis resulted in 24 AM fungal taxa with significant indicator values > 0.25 with respect to either habitat type or broad ecosystem type.

Table 2

Results of indicator species analysis showing AM fungal VT characteristic (indicator value > 0.25) of particular habitat types (X) and/or broad ecosystem distinctions – open vs. forested (shaded cells)

VT	Genus/group		Species			Open			Forested
Intensive				Sustainable					Grassland						Plantation							Forest
VTX00245 (−)	Archaeosporaceae		Archaeospora trappei
VTX00008	Archaeosporaceae		Archaeospora sp.			X
VTX00009	Archaeosporaceae		Archaeospora sp.
VTX00056 (ns)	Claroideoglomeraceae		Claroideoglomus sp.
VTX00193 (ns)	Claroideoglomeraceae		Claroideoglomus lamellosumg																		X
VTX00060 (−)	Diversisporaceae		Diversispora celata, eburnean							X
VTX00061	Diversisporaceae		Diversispora sp., Glomus versiforme																		X
VTX00062 (ns)	Diversisporaceae		Diversispora sp.							X
VTX00306	Diversisporaceae		Diversispora sp.							X
VTX00067 (+)	Glomeraceae		Glomus mosseae
VTX00074 (+)	Glomeraceae		Glomus sp.																									X
VTX00113 (−)	Glomeraceae		Glomus fasciculatum, intraradices																									X
VTX00125 (+)	Glomeraceae		Glomus sp.												X
VTX00130 (+)	Glomeraceae		Glomus sp.												X
VTX00137	Glomeraceae		Glomus sp.												X
VTX00140 (+)	Glomeraceae		Glomus sp.												X
VTX00142	Glomeraceae		Glomus sp.																									X
VTX00156 (ns)	Glomeraceae		Glomus sp.												X
VTX00160 (+)	Glomeraceae		Glomus sp.																									X
VTX00166 (+)	Glomeraceae		Glomus sp.																									X
VTX00191 (+)	Glomeraceae		Glomus sp.																									X
VTX00199 (+)	Glomeraceae		Glomus hoi, macrocarpum																									X
VTX00222 (+)	Glomeraceae		Glomus indicum												X
VTX00281 (+)	Paraglomeraceae		Paraglomus laccatum

VT names are followed by (+) when the indicator status from this study is in accordance with the corresponding habitat association of records in the MaarjAM database (determined using randomization; see Methods and Table S2). Inconsistency between the two data categorizations and nonsignificant results in the MaarjAM analysis are indicated by (−) and (ns), respectively. Remaining VT had < 10 records in the MaarjAM database and were not included in the randomisation analysis.

We investigated how these 24 indicator taxa are distributed among habitats in the MaarjAM database of global AM fungal diversity. Among the 24 indicators, 18 were represented by > 10 records in the MaarjAM database and were included in randomisation analysis. In 11 cases of 18, the local indicator status of the AM fungal taxa was matched by a corresponding significant positive habitat association in the MaarjAM database; in three cases, a negative habitat association in the MaarjAM database indicated inconsistency between the classifications; and in four cases, the habitat association in the MaarjAM database was neither significantly positive nor negative (Table 2, Fig. S1).

Differences in mean community composition at sites were significantly correlated with soil pH (R = 0.42; P < 0.01) and P content (R = 0.35; P = 0.01) in Mantel tests, with marginally nonsignificant correlation for organic C (R = 0.20; P = 0.06).

Phylogenetic community structure

The phylogenetic placement of VT occurring in different habitats was determined using the phylogenetic tree from Öpik et al. (2013) which contains all known Glomeromycota VT (MaarjAM database type sequences; Fig. 4). This indicated that Acaulosporaceae and some clades within Glomeraceae were underrepresented in all studied habitat types. It was also apparent that the phylogenetic profile of taxa detected in plantation soils differed considerably from all others, with a particular underrepresentation of Glomeraceae.

Fig. 4

Glomeromycota phylogenetic tree with taxa recorded in different habitats highlighted at the tips. The tree contains type SSU rRNA gene sequences of VT from the MaarjAM database and is taken from Öpik et al. (2013). Circles at the tree tips indicate the presence of VT in the habitat types listed above the tips; the size of circles is proportional to the mean relative abundance of taxa in the respective habitat.

nmds and permanova based on mean phylogenetic distance (mpd) between taxa in pairs of samples indicated a significant effect of habitat (permanova pseudo-F = 4.3, P = 0.002, R² = 0.24) and mechanical disturbance (permanova pseudo-F = 2.7, P = 0.009, R² = 0.16) on phylogenetic composition (Fig. 3b), while the effect of broad habitat type (open vs. forested) was marginally nonsignificant (permanovaF = 2.3, P = 0.11, R² = 0.06). The permanova habitat results broadly reflected differences in multivariate dispersion: Sample-to-sample turnover in phylogenetic composition was significantly lower in plantation (0.20) and forest (0.31) than in the agricultural habitats (0.51–0.57), while that of grassland (0.36) was significantly lower than that of intensive agriculture (Tukey's HSD test; Fig. 3b).

Mean phylogenetic distances between sites were not significantly correlated with any measured soil variables (Mantel tests).

Discussion

We found that habitats with different types and intensities of land use harbour significantly different AM fungal communities. The taxonomic composition of fungal communities varied to reflect both differences between forested vs. structurally open habitats and mechanically disturbed vs. nondisturbed habitats. At the same time, differences between habitats in the phylogenetic composition of AM fungal communities reflected most clearly the degree of mechanical disturbance occurring in the habitat.

Richness of AM fungal communities

Structurally open ecosystems tended to exhibit higher AM fungal taxon richness than forested habitats. Permanent grasslands were characterized by a diverse AM fungal community, as has been reported by earlier studies (Vandenkoornhuyse et al., 2002a, b; Saito et al., 2004; Horn et al., 2014). When different types of arable land were compared, the impact of land-use intensity was not so evident – sustainable fields exhibited slightly higher AM fungal richness than intensive fields, but the difference remained statistically nonsignificant. This result contrasts with the findings of many previous studies (see citations in Introduction), although similar results have also been reported (Jansa et al., 2002; Mathimaran et al., 2007).

Relatively low AM fungal richness was recorded in forest ecosystems. The same pattern was reported in a review by Öpik et al. (2006). Although the total number of AM fungal taxa inhabiting plant roots in boreonemoral coniferous forest may be fairly high (Öpik et al., 2009; Davison et al., 2011; Saks et al., 2014), the diversity per soil sample tends to be moderate (Davison et al., 2012). There might be several mechanisms underlying this pattern. The pH of the topsoil in the forest ecosystems addressed in this study was lower than that in the open ecosystems. Indeed, forest soils under coniferous tree species, especially Norway spruce, are generally more acidic than those under deciduous trees (Augusto et al., 2003), and invading conifers can decrease soil pH in former grasslands (Cumming & Kelly, 2007). Second, forested habitats were predominated by ectomycorrhizal Norway spruce. Previous work has shown that ectomycorrhizal trees can suppress AM fungal diversity in soil (Tyndall, 2005) and root AM fungal colonization (Becklin et al., 2012). In addition, Norway spruce reduces the diversity and abundance of the predominantly AM herbaceous forest field layer (Koorem & Moora, 2010) due to the inhibiting effect of spruce litter on seed emergence and seedling establishment (Koorem et al., 2011). The lower AM fungal taxon richness in forested compared to open habitats may thus be the result of the complex influence of coniferous trees through their associated ectomycorrhizal fungi, through acidification of the topsoil due to needle litter or through suppression of predominantly AM plant species in the understory.

Within forested habitats, spruce plantation exhibited considerably lower AM fungal taxon richness than primeval forest, which is consistent with a generally adverse effect of intensive land use on AM fungi. Although forest management is known to influence the structure of microbial communities in soil (Mummey et al., 2010), information about its impact on AM fungal diversity is very scarce. Öpik et al. (2008) did not record an effect of past logging on AM fungal richness inhabiting plant roots in a boreonemoral forest, while Bennett et al. (2013) demonstrated that AM fungal network structure changes with forest age. Even-aged Norway spruce stands alter the structure (Moora et al., 2009) and reduce the diversity (Koorem & Moora, 2010) of the herbaceous field layer, resulting in significantly different understory plant community in plantations compared to unmanaged old-growth forest (Moora et al., 2007). This may have an effect on the abundance and diversity of the soil AM fungal communities and may have resulted in the less diverse AM fungal communities in plantation.

We are not aware of other comparisons of AM fungal communities in natural forest and plantation in the temperate zone.

AM fungal community composition

Earlier global analyses have revealed varying composition of AM fungal communities between open and forested habitat types (Öpik et al., 2006; Kivlin et al., 2011). Here, we also observed such a difference in soil AM fungal communities at the local scale. There is also evidence that AM fungal communities exhibit different structure depending upon habitat soil characteristics (Dumbrell et al., 2010; Lekberg et al., 2011; Meadow & Zabinski, 2012). In common with these studies, we found that community composition varied along a topsoil pH gradient, but we cannot differentiate the effects of vegetation which also influences topsoil pH. We also noted a significant change in AM fungal community composition along a soil P gradient. Previously, Bainard et al. (2014) found that although phosphate flux was negatively correlated with AM fungal diversity in semiarid prairie agricultural soils, it did not significantly influence the composition of AM fungal communities.

There was also significant within-broad-ecosystem-type variation in AM fungal community composition, reflecting land-use intensity. Soil AM fungal communities differed between intensive and sustainable arable land. A similar pattern has been recorded by Lumini et al. (2010) and Miras-Avalos et al. (2011). AM fungal communities also differed between primeval forest and spruce plantation. Again, existing information about AM fungal community patterns in forests with different land-use intensity is scarce. Öpik et al. (2008) did not find differences between formerly logged and old-growth stands in boreonemoral coniferous forest. Meanwhile, Haug and colleagues found differences between some pristine tropical forest and afforestation sites (Haug et al., 2010), but not between others (Haug et al., 2013).

We used the MaarjAM database to investigate the global habitat associations of indicator taxa recorded in this study. This revealed that 61% of the local indicator VT was also significantly overrepresented in the corresponding habitats in the global data. This suggests that at least a proportion of AM fungal taxa exhibit relatively narrow ecological requirements that generate consistent variation in abundance at a large geographical scale. In addition to life-history traits (Chagnon et al., 2013), information about the ecology of AM fungi can be an invaluable component of future AM fungal community analyses.

Phylogenetic community structure

In contrast to variation in taxonomic composition, differences in mean phylogenetic distance between AM fungal communities were not significantly correlated with any measured soil variables in this study; neither did phylogenetic community structure differ significantly between forested and agricultural habitats. Changes in phylogenetic community structure were most apparent in the comparison of mechanically nondisturbed and disturbed habitats. Nondisturbed habitats, notably grassland and primeval forest, were predominated by Glomeraceae, although the taxonomic identity of the characteristic taxa differed; meanwhile, Glomeraceae were considerably less represented among the indicator taxa for other, less stable habitats. This pattern was also apparent from the habitat associations of local indicator taxa in the global MaarjAM database (Öpik et al., 2010). Phylogenetic similarity is expected to indicate a degree of functional similarity (Hart & Reader, 2002; Maherali & Klironomos, 2007; Chagnon et al., 2013). For instance, root colonization in Glomeraceae takes place primarily from hyphal fragments rather than from spores as in Gigasporaceae (Hart & Reader, 2002), and their capacity for dispersal and recolonization after disturbance might be consequently low (Hart & Reader, 2004). In the context of primary succession, Sikes et al. (2012) also reported higher representation of Glomeraceae in later successional stages. Previous studies have shown a significant impact of soil disturbance on AM fungal community composition (Helgason et al., 1998; Fitzsimons et al., 2008; Miras-Avalos et al., 2011; Schnoor et al., 2011; Wetzel et al., 2014). The current study indicates that the disturbance regime might be an important determinant of AM fungal community phylogenetic structure and hence also functional structure.

Spruce plantation was considered a disturbed habitat, and it indeed differed from primeval forest and permanent grassland in the representation of Glomeraceae. Although the level of mechanical soil disturbance in plantation is considerably lower than that occurring in arable fields, some disturbance occurs due to regular thinning and the subsequent transport of logged trees. This could have contributed to the difference in phylogenetic composition between plantation on one hand and primeval forest and permanent grassland on the other.

Our study reveals significant variation in AM fungal community diversity and taxonomic and phylogenetic composition under different land-use regimes, within both forested and open habitats. Taxonomic composition of AM fungal communities tends to be relatively more influenced by broad habitat type, while phylogenetic composition tends to respond relatively more to land-use intensity (i.e. the frequency of soil disturbance).

Acknowledgements

This research was funded by grants from the Estonian Science Foundation (9050, 9157), targeted financing (IUT 20-28), EU LIFE+ project 08 ENV/EE000258, European Union 7th framework project SCALES (FP7-226852), the European Regional Development Fund (Centre of Excellence FIBIR) and Environmental Protection and Technology R&D programme project ERMAS.

Author notes

Editor: Ian C. Anderson

Estonian Science Foundation

9050

9157

IUT 20-28

EU LIFE+

ENV/EE000258

European Union

FP7-226852

European Regional Development Fund (Centre of Excellence FIBIR)

Environmental Protection and Technology R&D