Mycorrhiza DOI 10.1007/s00572-014-0624-1

ORIGINAL PAPER

Mycorrhizal fungal communities in coastal sand dunes and heaths investigated by pyrosequencing analyses Synnøve Botnen & Håvard Kauserud & Tor Carlsen & Rakel Blaalid & Klaus Høiland

Received: 17 October 2014 / Accepted: 23 December 2014 # Springer-Verlag Berlin Heidelberg 2015

Abstract Maritime sand dunes and coastal ericaceous heaths are unstable and dynamic habitats for mycorrhizal fungi. Creeping willow (Salix repens) is an important host plant in these habitats in parts of Europe. In this study, we wanted to assess which mycorrhizal fungi are associated with S. repens in four different coastal vegetation types in Southern Norway, three types from sand dunes and one from heaths. Moreover, we investigated which ecological factors are important for the fungal community structure in these vegetation types. Mycorrhizal fungi on S. repens root samples were identified by 454 pyrosequencing of tag-encoded internal transcribed spacer 1 (ITS1) amplicons. Significantly higher fungal richness was observed in hummock dunes and dune slacks compared to eroded dune vegetation. The compositional variation was mainly accounted for by location (plot) and vegetation type and was significantly correlated to content of carbon, nitrogen and phosphorus in soil. The investigated maritime sand dunes and coastal ericaceous heaths hosted mycorrhizal taxa mainly associated with Helotiales, Sebacinales, Thelephorales and Agaricales. Keywords Mycorrhiza . Salix repens . Ericaceous heaths . Sand dunes

Introduction Maritime sand dunes and coastal ericaceous heaths are among the most dynamic terrestrial ecosystems in temperate Northern Electronic supplementary material The online version of this article (doi:10.1007/s00572-014-0624-1) contains supplementary material, which is available to authorized users. S. Botnen : H. Kauserud : T. Carlsen : R. Blaalid : K. Høiland (*) Section for Genetics and Evolutionary Biology, Department of Biosciences, University of Oslo, Oslo, Norway e-mail: [email protected]

Europe (Gimingham 1972; Ranwell 1972). The sand dunes harbour some characteristic macrofungi, including Phallus hadriani and Psathyrella ammophila, which are saprotrophic, and Inocybe serotina, Inocybe dunensis and Laccaria maritima, which form ectomycorrhiza (ECM) (Andersson 1950; Høiland 1975; Watling and Rotheroe 1989). Based on fruit body surveys, the ericaceous heaths seem to be rather species poor (Høiland 1981), which could be because the tannic acids in raw humus with Calluna are believed to suppress ECM fungi (Gimingham 1972; Bending and Read 1996). Creeping willow (Salix repens) is often the only ECM plant in northern European maritime sand dunes and coastal ericaceous heaths (Gimingham 1972; Høiland and Elven 1980; Watling and Rotheroe 1989; Arnolds and Kuyper 1995; van der Heijden and Vosatka 1999; van der Heijden et al. 1999). In Norwegian sand dune ecosystems, its ecology is wide, but it is most important in the following three vegetation types (Høiland 1978; Høiland and Elven 1980): (1) dune slacks, temporarily flooded, flat depressions with a more or less dense layer of S. repens and plant species typical for moist areas; (2) hummock dunes, dome-shaped dunes with dense S. repens and species typical for dry areas; and (3) eroded dunes, demolished (by wind or man) hummock dunes with impoverished S. repens and no or few other species. As Watling 2005 pointed out, S. repens can grow in environments which are basically unfavourable for macromycetes, such as high salt concentrations, high temperatures due to high solar radiation and low water availability. It produces a shrub layer enriched with litter in which different individuals or even species may be interconnected by common mycorrhizal networks (Read 1989). In addition to ECM fungi, S. repens may host arbuscular mycorrhiza (AM), depending on the relative availability of nitrogen and phosphorus in the soil (Read 1989; van der Heijden 2001; Watling 2005). Mycorrhizal associations in S. repens have been thoroughly studied in the Netherlands (van der Heijden and Vosatka

Mycorrhiza

1999; van der Heijden et al. 1999; van der Heijden 2001; van der Heijden and Kuyper 2001a, b, 2003). These works emphasize the difference in habitat preference of various ECM morphotypes and AM fungi, and it has been suggested that the mycorrhizal symbionts contribute to its broad ecological amplitude. The mentioned studies on mycorrhizal associations of S. repens have mainly been based on traditional morphotyping. High-throughput sequencing techniques, such as 454 pyrosequencing, enable more in-depth analyses of belowground fungal communities (Lindahl et al. 2013). Soil samples from the rhizosphere of S. repens in Dutch sand dunes have recently been investigated by 454 pyrosequencing by Geml et al. (2014). The aim of this study was (1) to record mycorrhizal fungal groups and their diversity in roots of S. repens by means of environmental DNA sequencing, (2) to investigate the changes in richness of mycorrhizal fungi in sand dunes and ericaceous heaths, (3) to investigate the changes in fungal composition regarding the variable ecology in these dynamic systems and (4) to reveal whether the mycorrhizal communities change temporally within a growing season.

Material and methods Field work Field work was carried out on the Lista peninsula in SW Norway. In October 2008, nine study plots including S. repens were established in the sand dune areas in the S and SE parts of the Lista encompassing (1) dune slacks, (2) hummock dunes and (3) eroded dunes (Høiland 1978; Høiland and Elven 1980) and three in an area with grazed (sheep) ericaceous heaths in the NW part of the peninsula (Fig. 1). A stratified sampling was done according to the four above-mentioned vegetation types, where three plots were established in each vegetation type (Fig. 1). Each plot included an area delimited by a 3-m-diameter circle. The midpoint geographic coordinates were obtained for each plot. Two sampling sites of 0.25 m2, A and B, were placed within each plot. For an attempt to insure independence, the distance between the samples was always more than 1 m. For details concerning the 12 plots, see Online Resource 1. In 2009, all plots were visited three times, in April, July and October, in order to characterize the fungal diversity in spring, summer and autumn. From sites A and B, an approximately 10 cm of continuous S. repens root was collected at each time point. Hence, each plot comprised six root samples, i.e. 72 root samples in the total data set. For further analyses, the data from the two sampling sites, A and B, were combined, reducing the number to three seasonal samples per plot, altogether 36 samples in the total data set. Each vegetation type (heaths, dune slacks, hummock dunes and eroded dunes) is therefore

represented by three plots, each with three seasonal samples, altogether nine samples per vegetation type. Within 24 h, each root was carefully rinsed. Five to ten root fragments with healthy ECM root tips were detached from the cleaned root under a dissection microscope and pooled into a 2-mL Eppendorf tube containing 2 % cetyltrimethylammonium bromide (CTAB) buffer. The samples were stored at −20 °C until DNA extraction. In July 2009, soil samples from each site, A and B, were collected and mixed according to plot and 100–500 g was dried in indoor temperature and stored in paper bags before analysis of total carbon (C), nitrogen (N) and phosphorus (P). Chemical and DNA analyses Total C and N were measured by rapid combustion of 8– 10 mg of the sample in pure oxygen using a Thermo Finnigan EA 1112 Series Flash Elemental Analyzer. For measuring total P, 10–14-mg of the sample was mixed with 10 mL 1 % K2S2O8 for 30 min in 121 °C and then analyzed using the colorimetric method by Murphy and Riley (1962) using Bran + Luebbe AutoAnalyzer III, method No. G-297-03 (Multitest MT). DNA was extracted from root samples following the CTAB miniprep method described by Murray and Thompson (1980) combined with the purification steps from E.Z.N.A soil DNA Kit (Omega Bio-Tek, Doraville, GA, USA). Samples were prepared for 454 pyrosequencing by performing nested PCR amplification using the fungus-specific primers internal transcribed spacer 1F (ITS1F) and ITS4 (White et al. 1990; Gardes and Bruns 1993) in the first step and fusion primers including ITS5 and ITS2 (White et al. 1990) in the nested step. Fusion primers were constructed by adding five different unique 10-bp tags and 454 pyrosequencing adaptors A and B to ITS5 and ITS2, respectively. The same tags were added to both forward and reverse primers. PCR was performed on MJ Thermal Cycler PTC-200 in 20-μL reactions containing 2 μL template DNA and 18 μL reaction mix. Final concentrations were 0.10 mM dNTP mix, 0.125 μM of each primer and 0.5 units polymerase (Phusion Hot Start II, Finnzymes, Vantaa). The PCR amplification program was as follows: 30 s at 98 °C, followed by 20 cycles of 10 s at 98 °C, 20 s at 50 °C and 20 s at 72 °C and a final extension step at 72 °C for 7 min before storage at −20 °C. The nested PCR was run with the same reaction concentrations and amplification program, but with a 50× diluted PCR mix as a template. PCR products were pooled into eight equimolar amplicon libraries using SequalPrep™ Normalization Plate (96) Kit following the manufacturer’s protocol (Invitrogen, CA, USA) and cleaned with Wizard® SV Gel and PCR Clean-Up System (Promega, Madison, USA); 454 titanium sequencing of the tagged amplicons was performed at the Norwegian HighThroughput Sequencing Centre (http://www.sequencing.uio.

Mycorrhiza Fig. 1 Sketch map of the investigated area, the Lista peninsula, Norway, Vest-Agder County, Farsund Municipality. H1, H2 and H3 are plots from ericaceous heaths; DS1, DS2 and DS3 are plots from dune slacks; SRD1, SRD2 and SRD3 are plots from hummock dunes; SRDE1, SRDE2 and SRDE3 are plots from eroded dunes

no/) using a 454 plate divided into 16 compartments. The raw data were submitted to the European Nucleotide Archive SRA0307041. Bioinformatics analyses Tag switching has been shown to be a source of error in pyrosequencing (Carlsen et al. 2012). To remove reads with nonmatching 3′ and 5′ tags, we used a python script available on https://github.com/karinlag/karinlag_utils/tree/original. Further quality filtering and de-noising were conducted using Qiime (Caporaso et al. 2010). Reads of 500 bp, with a quality score of less than 25 in a sliding window of 50 bp, with more than one ambiguous base (N) and with homopolymers >10 bp, were discarded. One-base pair mismatch was allowed in both the forward and reverse primer. De-noising was conducted using default settings in Qiime. The remaining reads were clustered into operational taxonomic units (OTUs) at a 97 % similarity level using the UCLUST algorithm (Edgar 2010) as implemented in Qiime. In most comparable studies, a 97 % clustering level has been used. It is an inherent problem

in this type of analysis that some OTUs include multiple species and that some species are split into multiple OTUs (see Blaalid et al. 2013). This is because there is no general cut-off level that works well for all taxa. While the clustering level directly affects the number of OTUs (i.e. alpha and gamma diversity), it does, to a less extent, affect the compositional turnover (i.e. beta diversity), as demonstrated by Lekberg et al. (2014) for AM fungi. OTUs containing only one sequence were removed, as they are likely sequencing errors. Numerous studies have shown that singletons to a large degree represent PCR/sequencing errors and that the inclusion of singletons may lead to severe overestimation of the richness (see e.g. Kunin et al. 2010). Sequences were considered as chimeric and removed from the data set if they were identified by the Perseus algorithm in mothur (Schloss et al. 2009) and had a BLASTn hit with 85 % identity and >85 % coverage against known AM or ECM fungi for

Mycorrhiza

further analyses. Fungi were considered ECM if they were listed by Tedersoo et al. (2010). The sample OTU matrix was rarefied to 400 reads per sample. To test if any of the OTUs were significantly correlated with the different vegetation types or seasons, G tests of independence were carried out on the obtained OTU matrix. To be considered in the G test, an OTU had to be present in at least five samples; apart from this, the tests were conducted with default settings in Qiime. Downstream statistical analyses were conducted in the statistical environment R (R Development Core Team 2010). All species richness and multivariate analyses were performed using the vegan package (Oksanen et al. 2012). Accumulation curves were calculated from the rarefied OTU matrix curves using the approach described in Ugland et al. (2003). Extrapolated species richness was calculated using the estimators first-order Jackknife, second-order Jackknife and bootstrap. To test whether the OTU richness (here defined as a number of OTUs per sample) was influenced by vegetation type, a general linear model (GLM) with quasi-Poisson distribution was performed. Samples from the same plot, but from different seasons, were merged to exclude temporal dependence. The model included 12 samples and 8 degrees of freedom. A total of four models were made, and they were all similar, except that they had different reference levels: one for each of the different vegetation types. This was done to check all combinations of vegetation types against each other. For further analyses, the non-rarefied sample OTU matrix was transformed to a present/absent matrix. A global nonmetric multidimensional scaling (GNMDS) and a detrended correspondence analysis (DCA) were conducted in parallel on the OTU matrix. The GNMDS was conducted using the following settings, as suggested by Liu et al. (2008): dissimilarity measure (Bray-Curtis), maximum iterations (200), dimensions (2) and initial configuration (100). The DCA was conducted using default settings. The ordinations were screened for outliers and possible artefacts like the tongue effect in the DCA and arch effect in the GNMDS. To test for correlations between the axes/dimensions in the ordinations, Kendall’s τ coefficient tests were conducted and the tests were used as an indicator that reliable gradients had been recovered. The environmental variables (N, C and P) were standardized using zero-skewness standardization (Sokal and Rohlf 1995) to improve the homogeneity (Økland et al. 2001) and bring the variables to the same scale. The environmental variables were then fitted to the GNMDS with the envfit function in vegan. Correlation tests between the dimensions in the GNMDS and the environmental variables were conducted, using Kendall’s τ coefficient. To estimate the variation attributed to different variables, constrained ordinations (CCA in R) were preformed, using a forward selection approach on the variables: sample site, vegetation type, season and the environmental variables. A separate GNMDS was conducted on a

sample OTU matrix excluding samples from the heath. Furthermore, a GNMDS, with the settings described above, was conducted on a present/absent OTU matrix including only common OTUs that occurred in at least five samples. OTU (species) scores were calculated from plot scores and OTU abundance by weighted averaging regression, using the function wascores(). This was done to get an indication on whether most common OTUs were similar when it came to which samples they were in.

Results Data properties After removing sequences with inconsistent tags, quality filtering and de-noising, 118,687 of the original 140,855 reads remained for further analyses. These clustered into 1190 OTUs, of which 351 contained only one read (singletons). Nineteen OTUs were considered to be chimeric, and 646 had top BLAST hits to what were considered nonmycorrhizal fungi, and were removed from further analyses since we chose to focus only on mycorrhizal fungi in this study. One sample (from the ericaceous heath) was excluded because of a low number of reads. A total of 170 OTUs (Online Resource 2) considered to represent ECM or AM fungi, appearing in 35 samples, were kept for further analyses in the non-rarefied Hellinger transformed OTU matrix. In the rarefied OTU matrix (400 sequences per sample), 148 OTUs in 34 samples were kept. Fungal richness On average, each OTU occurred in 5.18 samples (range 1–25) and the mean number of OTUs in each sample was 20.6 (range 6–38). The accumulation curve of OTU richness versus sampling effort did not reach asymptote but levelled off, and the different species richness estimates (first-order Jackknife, second-order Jackknife and bootstrap) indicated that 75.1– 87.6 % of the diversity was captured (Fig. 2). The average number of OTUs in each vegetation type was 38. In the heaths, the average was 37.6 OTUs; in dune slacks, 41.3 OTUs; in hummock dunes, 46.3 OTUs; and in eroded dunes, 27 OTUs. Although few replicates were included, the GLM analysis revealed a significantly higher number of OTUs in hummock dunes and dune slacks than in eroded dunes (p