Linking the community structure of arbuscular mycorrhizal fungi and plants: a story of interdependence?

Arbuscular mycorrhizal fungi (AMF) are crucial to plants and vice versa, but little is known about the factors linking the community structure of the two groups. We investigated the

association between AMF and the plant community structure in the nearest neighborhood of Festuca brevipila in a semiarid grassland with steep environmental gradients, using high-throughput

sequencing of the Glomeromycotina (former Glomeromycota). We focused on the Passenger, Driver and Habitat hypotheses: (i) plant communities drive AMF (passenger); (ii) AMF communities drive

the plants (driver); (iii) the environment shapes both communities causing covariation. The null hypothesis is that the two assemblages are independent and this study offers a spatially

explicit novel test of it in the field at multiple, small scales. The AMF community consisted of 71 operational taxonomic units, the plant community of 47 species. Spatial distance and

spatial variation in the environment were the main determinants of the AMF community. The structure of the plant community around the focal plant was a poor predictor of AMF communities,

also in terms of phylogenetic community structure. Some evidence supports the passenger hypothesis, but the relative roles of the factors structuring the two groups clearly differed, leading

to an apparent decoupling of the two assemblages at the relatively small scale of this study. Community phylogenetic structure in AMF suggests an important role of within-assemblage

interactions.

Arbuscular mycorrhizal fungi (AMF) are one of the most important symbiont groups for plants, forming relationships with the majority of land plants and having a significant role in the

acquisition of phosphorus (Smith and Read, 2008). Yet, despite some important progress in recent years, especially in relation to interactions with other soil biota or how AMF respond to

management (Alguacil et al., 2014; Caravaca and Ruess, 2014; Leifheit et al., 2015; Knegt et al., 2016), there are many aspects of the assembly processes regulating the community ecology of

these organisms that are poorly understood: a key challenge remains disentangling the relative contribution of dispersal limitation, environmental filtering and biotic interaction on AMF

community structure (Vályi et al., 2016). The cryptic nature of the group and the complexity of the three-way interaction between plants, AMF and the environment complicate the study of the

factors that regulate AMF community structure. Dispersal limitation remains one of the most complex aspects of AMF ecology (Zobel and Öpik, 2014): as for example reviewed in Vályi et al.

(2016), AMF can disperse via local mycelium spread but also spores, hyphal fragments and colonized root fragments, and the importance of these mechanisms could be scale-dependent, although

direct evidence is missing. Still, large AMF spores and hyphal fragments are mostly spread via zoochory, which implies limited dispersal capability and seems reflected by small-scale

patterns in community structure (Mummey and Rillig, 2008; Dumbrell et al., 2010a; Horn et al., 2014). The effects of dispersal limitations are entangled with those of environmental

gradients, biotic interactions within the AMF assemblage, and between AMF and plants (for example, Mummey and Rillig, 2008; Dumbrell et al., 2010a; Horn et al., 2014; Martinez-Garcia et al.,

2015; García de León et al., 2016a, 2016b).

The study of AMF in grasslands is of particular importance as grassland ecosystems cover a significant proportion of the earth’s surface, harbor the majority of herbaceous plant diversity

(Shantz, 1954) and it is in grasslands that AMF reach their highest abundance and diversity (Treseder and Cross, 2006; Kivlin et al., 2011). Studies on plant biodiversity in grassland

ecosystems at small scales have revealed connections between species richness of AMF and plants (Hiiesalu et al., 2014) and host plant effects on AMF community composition (Vályi et al.,

2015). Still, effects can be very localized: AMF can form extended hyphal networks but spatial autocorrelation in their distribution is typically found at submeter scales (Mummey and Rillig,

2008), with a potential role for biotic interactions (Vályi et al., 2016). To date, only a few studies have taken this fact into account and applied a sufficiently fine-grained sampling

design for a solid statistical analysis of the patterns generated by local processes (Dumbrell et al., 2010b; Horn et al., 2014).

AMF and plants form two sets of communities associated with each other but assembled through different processes that take place at different spatial and temporal scales (Zobel and Öpik,

2014). The plant set can drive the fungal set or vice versa (Figure 1), but which group is driving might depend on successional stage, which is linked to differences in dispersal processes

between plants and AMF. Zobel and Öpik (2014) have used the concept of difference in dispersal between AMF and plants to revisit the Driver and Passenger hypotheses originally proposed by

Hart et al. (2001). Zobel and Öpik (2014) also formulated the Habitat hypothesis to distinguish a situation where AMF and plant communities covary but are not directly causally linked, as

opposed to the null hypothesis of no covariation (‘independence’). For example, during primary succession, plants typically arrive before AMF and then act as a potential filter to AMF: AMF

are Passengers as they are following plants. However, dispersal limitation in an established AMF assemblage can cause the AMF assemblage to more strongly determine which plants will

establish during secondary succession: the AMF assemblage becomes the Driver (Zobel and Öpik, 2014). Zobel and Öpik (2014) further predict that the Habitat hypothesis would be most common in

regions with a stable community (for example, climax vegetation) where environmental variation within regions will cause a non-mechanistic covariation between AMF and plant communities. The

general null hypothesis is that plants and AMF may vary independently of each other, which could possibly happen at very broad or global scales, where plants are more disperal limited than

AMF seem to be (Kivlin et al., 2011; Öpik et al., 2013; Davison et al., 2015). Accordingly, Vályi et al. (2016) have recently proposed that the host effect is minimal at regional and global

scales.

Autocorrelation (Semivariogram) and trends in environmental variables create (arrow a) spatial structure and environmental gradients. Variation in the environment generates variation in

plants and AMF (arrows b). AMF and plants can thus be structured by changes in habitat conditions, which can then simply lead to covariation between the two assemblages (Habitat hypothesis).

Alternatively, AMF could either drive the plant assemblage (Driver hypothesis, arrow c) or be driven by the plant assemblage (Passenger hypothesis, arrow d). In all cases, the driving

factors/assemblage (b–d) have a spatial structure that will be, at least partially, reflected by spatial structure in the driven assemblage. This spatial dependence calls for a spatially

explicit approach to the testing of the three hypotheses. Spatial scale and successional stage have also been hypothesized to be the major factors in determining which among the Habitat,

Driver and Passenger hypotheses apply to real systems. In addition to all these factors, AMF can also be structured by interactions within the assemblage, independently of plants, which has

been hypothesized to happen at local scale and that could create very patchy distribution. All data are simulated.

There are studies that have touched upon components of these hypotheses. For example, AMF taxa are generally found to be able to colonize any AM (as opposed to non-AM) plant species

(Klironomos, 2000), still there may be a bias towards easily cultivable species (Ohsowski et al., 2014) and ‘specificity’ might be quantitative rather than qualitative (Vályi et al., 2015).

Therefore, AM fungal communities and plant communities may still be directly causally correlated despite the perceived generalism of the AM symbiosis. A thorough account of the studies

supporting the various hypotheses is given in Zobel and Öpik (2014) and we are aware of only two recent, observational studies that have addressed the subject (Martinez-Garcia et al., 2015;

García de León et al., 2016a). However, a problematic aspect of observational field studies remains to tease apart cause and effect in the correlations between the two organism groups in the

presence of spatial structure in the environment (Figure 1). To solve this problem, we applied a spatially explicit design to sample AMF and plant communities along a replicated steep but

short (≈15 m) soil environmental gradient (Horn et al., 2014). We could therefore control for spatial patterns and environmental effects when testing for the effects of plants on AMF

communities and vice versa. We used a standardized focal plant of high abundance to investigate environmental, plant and AMF community variation at sufficiently small scales. We also took

into account the phylogenetic community structure of both plant and AMF assemblages to allow community relationships to occur at levels other than species/operational taxonomic unit (OTU)

between and within the groups.

Our main aim was to collect for the first time multiple scales and high spatial resolution data to test the general null hypothesis that plant community structure, including phylogenetic

structure, is independent of AMF community structure and vice versa. If the hypothesis were rejected, given the scales included in the study, we aimed to collect support for one or more of

the three alternative hypotheses (Figure 1), with the overall goal of shedding light on the mutual relationships between plant and AMF communities.

Sampling was conducted in a nature protection area located in north-eastern Germany (Brandenburg, 52°27.778'N, 14°29.349'E), a Natura 2000 biodiversity hotspot that contains over 200

different plant species and combines floral elements of steppes and coastal habitats. Given the high diversity of plants (Ristow et al., 2011) and AMF (Horn et al., 2014), the area is very

suitable for this study. We sampled by a hierarchical nesting of plots in April 2011: twelve 3 × 3 m2 plots were sampled at the four corners of three 15 × 15 m2 larger plots (henceforth

called ‘macroplots’) located on the slope of a hillside (Supplementary Figure S1). The distances between the macroplots ranged from 20 to 500 m (Supplementary Figure S2), leading to overall

intersample distances from a few cm to 3 m (within a plot) and up to 500 m between macroplots. The uphill–downhill axes of the three macroplots were characterized by a steep textural

gradient from sandy-loamy (uphill) to highly sandy (downhill) soils (Supplementary Figure S3). Soil parameters varied significantly and to a large extent (for example, almost 3 units of pH)

along the texture gradient (Horn et al., 2015).

We assessed the local AM fungal community in the roots and surrounding soil of Festuca brevipila plants plus the neighboring plant species around these Festuca plants. F. brevipila is one of

the most abundant species in sampled plots (Ristow et al., 2011; Horn et al., 2015). Soil cores (5 cm radius, 15 cm deep) were taken from five F. brevipila plants per plot, resulting in 60

(5 plants x 12 plots) sampling locations. Each sample position was random within the plot (minimum distance of 30 cm between any two samples in the same plot; Supplementary Figure S1). Plant

presence/absence was assessed in the surrounding area in a radius of 15 cm around each soil core to target local interactions present in the rhizosphere of our focal plant (neighborhood

plant community structure). This scale is consistent with the minimal observed spatial autocorrelation of AM fungi (30–100 cm; Mummey and Rillig, 2008).

Soil cores, including roots and plant material, were stored at −20 °C before analysis. Each soil core was thoroughly homogenized and subsampled for soil chemical analyses (Supplementary

Information, part a). We measured water content, pH, carbon, nitrogen and phosphorus content of the soil, which are known to affect AMF community variation (Camenzind et al., 2014; Horn et

al., 2014). Additionally, dehydrogenase activity was assessed as a proxy for microbial activity. Roots were washed in Millipore water before analysis.

We extracted genomic DNA twice from each core, once from 150 mg of washed, fine-ground F. brevipila roots and once from 250 mg of soil material, which was sieved through a 2 mm mesh. We used

the PowerSoil DNA Isolation Kit (MoBio Laboratories Inc., Carlsbad, CA, USA) following the procedure in the manufacturer’s manual. We then created 454 pyrosequencing amplicon pools for the

AMF using a nested PCR design, using the AMF-specific primer set SSUmAf and LSUmAr for the first and SSUmCf and LSUmBr for the second, nested PCR (Krüger et al., 2009). The amplified region

spans genes for the SSU (small ribosomal subunit), the complete ITS (internal transcribed spacer) region and a part of the LSU (large ribosomal subunit). Subsequently, amplicons of ~600 bp

in length were created from the AMF-specific PCR fragments using general fungal primers located in the LSU gene modified with 454 adapters and sample-specific barcode sequences

(Supplementary Information, part b). The 454 sequencing was carried out on a Roche GS FLX+ system with titanium chemistry at the Göttingen Genomics Laboratory at the Georg-August University

of Göttingen (Göttingen, Germany).

Sequences were denoized using the PyroNoise approach (Quince et al., 2009) implemented in Mothur (Schloss et al., 2009). The denoizing approach removes bad quality sequences, creates

sequence clusters and removes chimera sequences. After denoizing and preclustering, sequences from roots and soil were clustered into OTUs using CROP (Clustering 16S rRNA for OTU Prediction;

Hao et al., 2011), which uses a Bayesian clustering algorithm. This approach addresses species delineation uncertainty better than hierarchical clustering methods because of its flexible

cutoff, thereby creating significantly less artifact OTUs than fixed cutoff clustering approaches (Hao et al., 2011). We checked the final OTU sequences against chimeras using the Mothur

implementation of the uchime algorithm and the Krüger et al. (2012) SSU-ITS-LSU alignment, as well as the slayer algorithm against the sequences themselves. Default settings were used for

both algorithms.

Owing to the nature of pyrosequencing, we found differences in read numbers for every sampling location, so we resampled the read numbers to equal amounts of 500 reads per sample using a

bootstrap approach with 10 000 iterations per sample (Efron, 1979; Wehner et al., 2014). Samples with considerably lower read numbers than the estimated resampling threshold (