# Soil fungal community development in a high Arctic glacier foreland follows a directional replacement model, with a mid-successional diversity maximum.

Ke Dong, Binu Tripathi, Itumeleng Moroenyane, Woosung Kim, Nan Li, Haiyan Chu, Jonathan Adams

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ke Dong, Binu Tripathi, Itumeleng Moroenyane, Woosung Kim, Nan Li, et al., Soil fungal community development in a high Arctic glacier foreland follows a directional replacement model, with a mid-successional diversity maximum., Scientific Reports, 2016, 6, pp.26360. 10.1038/srep26360. pasteur-01351106

HAL Id: pasteur-01351106

https://riip.hal.science/pasteur-01351106

Submitted on 2 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# Soil fungal community development in a high Arctic glacier foreland follows a directional replacement model, with a mid-successional diversity maximum

Ke Dong<fnref n="1" />, Binu Tripathi<fnref n="1" />,², Itumeleng Moroenyane<fnref n="1" />,³,⁴, Woosung Kim<fnref n="1" />, Nan Li⁵, Haiyan Chu⁶ & Jonathan Adams<fnref n="1" />

Directional replacement and directional non-replacement models are two alternative paradigms for community development in primary successional environments. The first model emphasizes turnover in species between early and late successional niches. The second emphasizes accumulation of additional diversity over time. To test whether the development of soil fungal communities in the foreland of an Arctic glacier conforms to either of these models, we collected samples from the Midtre Lovénbreen Glacier, Svalbard, along a soil successional series spanning >80 years. Soil DNA was extracted, and fungal ITS1 region was amplified and sequenced on an Illumina Miseq. There was a progressive change in community composition in the soil fungal community, with greatest fungal OTU richness in the Mid Stage (50–80 years). A nestedness analysis showed that the Early Stage (20–50 years) and the Late Stage (>80 years) fungal communities were nested within the Mid Stage communities. These results imply that fungal community development in this glacier succession follows a directional replacement model. Soil development processes may initially be important in facilitating arrival of additional fungal species, to give a mid-successional diversity maximum that contains both early- and late-successional fungi. Competition may then decrease the overall diversity due to the loss of early successional species.

The processes by which communities and ecosystems assemble themselves in new or disturbed environments have long fascinated ecologists<fnref n="1" />. One of the most frequently studied successional environments is the foreland area uncovered by retreating glaciers²,³. Different models of community assembly have been tested and observed in primary succession⁴⁻⁶. During the 20th Century, two different models describing primary vegetation succession (directional replacement and directional non-replacement) were proposed and modified⁷,⁸. These two models have also been applied to invertebrate and fungi⁹,¹⁰. The dominant model of primary vegetation succession is essentially the directional replacement model, which emphasizes the replacement of primary colonizers by later stage species, and with a dominant role of competition in driving this community change¹¹. On the other hand, invertebrate succession in glacier forelands has been found to follow a directional non-replacement model - with the progressive addition and persistence of taxa, and little or no loss of species over time¹⁰.

A number of studies have examined fungal communities from forelands of retreating glaciers, and their results revealed that time since deglaciation has large effects on fungal community composition⁹,¹²⁻¹⁵. A recent study of the Midtre Lovénbreen glacier on Svalbard found that primary succession of root-associated fungi of

Received: 13 January 2016

Accepted: 29 April 2016

Published: 31 May 2016

<footnote n="1">¹Department of Biological Sciences, College of Natural Sciences, Seoul National University, Seoul, 08826, South Korea. ²Arctic Research Centre, Korea Polar Research Institute, Incheon, 21990, South Korea. ³Division of Life Sciences, School of Science, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, 999077, Hong Kong. ⁴Institut National de la Recherche Scientifique, Centre INRS-Institut Armand-Frappier, 531 boulevard de Prairies, Laval, Quebec, H7V 1B7, Canada. ⁵School of the Environment, Florida A&M University, Tallahassee, FL, 32307, USA. ⁶State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, East Beijing Road 71, Nanjing 210008 China. Correspondence and requests for materials should be addressed to J.A. (email: geograph.ecol@gmail.com)</footnote>
the ectomycorrhizal forb *Bistorta vivipara* (L.) Delabre (Polygonaceae) exhibited a 'directional non-replacement model' pattern of community change in following its host plant⁹. However, the root-associated fungal community alone is not necessarily representative of soil fungi in general. Root-associated fungi are always strongly influenced by their host plant, and so this community is particularly strongly tied to the species composition, ground cover and activity of vegetation⁹,¹⁶.

Soil fungal communities in successional systems are complex and they interact with many factors, beyond just the composition of the aboveground vegetation¹⁷⁻²⁰. For instance, the buildup of organic nutrients like soil carbon and nitrogen derived from primary plant colonizers is considered to be one important aspect of soil fungal community development as the density of primary plant colonizers increases in a successional system²¹. Stochastic factors depending on colonization of soil organisms can also strongly influence soil fungal communities in successional systems. For example, Jumponen² noted the importance of airborne spore deposition in the fungal communities. While the aerially deposited, dormant spore bank is present across the succession, it is masked in the older substrate areas by a larger active mycelial biomass. Therefore, the whole fungal community of the bulk soil of a glacial foreland would not be expected to follow strictly the same successional patterns as root-associated fungi.

We were interested in understanding the fungal communities of the Midtre Lovénbreen glacier foreland on Svalbard, against the background of the previously published study of root-associated fungi at the same site by Davey *et al.*⁹, to test whether soil fungal community variation along the chronosequence also follows the directional non-replacement model. However, our primary hypothesis was that the directional replacement model would instead be the best description of overall soil fungal community change, based on the following mechanisms: 1) changes in soil conditions such as pH and nutrient concentrations, as the glacier foreland soil develops over several decades, should favour different sets of fungal species in terms of physiological requirements. This will affect their viability in competition, causing turnover of species. 2) During the early stages of glacier foreland soil development, soil resources are patchily distributed, and the soil fungal community is dominated by opportunistic lottery competition²² which favours rapid colonization and population growth on resource patches relatively free of other competing species. As the ecosystem develops over time, soil resources will build up and become less patchily distributed: lottery competition becomes less important, and steady competition between species and their competitive hierarchy under this form of competition begins to play an increasing role in the community²³. This, we predicted, would eliminate fungal species that are generally weak competitors when placed under conditions of steady competition, causing loss of these species and dominance of communities by other species (which may have been slow to arrive and colonize) that are strong competitors. We also hypothesized that the soil fungal succession would relate strongly to soil total organic carbon (TOC) and total nitrogen (TN), as is the case with glacial succession generally for other groups of organisms and in glacier forelands elsewhere²¹,²⁴.

## Results

**Fungal community analysis.** A total of 161,200 high quality ITS1 region sequences were assigned to 774 OTUs at ≥97% similarity level from 31 samples after the removal of low quality, chimeric, and rare sequences. On average, 165 OTUs (±7.1 SE) were found in each sample. The individual highest OTU abundance was found in the Mid Stage with 406 OTUs, while the lowest OTU abundance was found in the Early Stage with 134 OTUs. The fungal communities in the deglaciated foreland area were dominated by Ascomycota which accounted for 79.2% of the classifiable reads, followed by Basidiomycota (5.5%), Zygomycota (4.1%) and Chytridiomycota (0.2%) or as basal fungal lineages. Basidiomycota and Chytridiomycota showed significant lower abundance in the Early Stage than Mid Stage and Late Stage (p < 0.05). A total of 4,636 sequences belonged to known groups of EcM fungi, representing 2.9% of the subsampled sequences. 11 genera were identified belonging to possible EcM fungi from our soil samples, with *Inocybe* as the most dominant genus (36.9% of total EcM sequences), followed by *Cortinarius* (28.8%).

Six OTUs had an abundance of more than 2% of the total sequences, and four of these (OTU00026, OTU00034, OTU00037, and OTU00046, all classified as Ascomycota) showed a significant or marginally significant difference between the Early Stage and the Late Stage (p ≤ 0.05). A total of 77 OTUs (18.5% of all generated reads) were potential indicator OTUs. Though the database cannot classify all of the sequences at a high taxonomic resolution, at a broad level most of these OTUs (63.6% of OTUs) were classified as Ascomycota. Basidiomycota as indicators were found only in the Late Stage. Among the indicator OTUs, 40 (12.6% of reads) were from the Early Stage; 18 (2.8% of reads) were from the Mid Stage and 21 (3.1% of reads) were from the Late Stage (Supplementary Table 1).

With increasing age since deglaciation, the fungal diversity showed a significant quadratic trend in relation to OTU richness (R = 0.449, p < 0.05) and a marginally significant quadratic regression by Shannon index (R = 0.426, p = 0.06) (Fig. 1). This revealed that OTU richness/diversity increased from the Early Stage of glacial succession to the Mid Stage, and then decreased from the Mid to Late Stage. Rarefaction curves of OTU richness also strongly supported the view that the Mid Stage had greater fungal diversity, with OTU richness per individual sample significantly greater than corresponding Early Stage and Late Stage samples (Supplementary Fig. S2). ANOSIM performed on Bray-Curtis community matrix revealed that samples from three different successional age categories varied significantly from each other (Table 1).

Nestedness analysis showed that the fungal communities across the whole time span of succession follow a nested structure (P ≪ 0.0001). A packed matrix order of all samples was generated, in which the nestedness of each sample was categorized from high to low, and the lower ones are nested in the higher ones (Supplementary Table 4). The Spearman's rho test indicated that this packed matrix order is highly correlated to the OTU richness of each sample (Fig. 2). Moreover, the OTU richness was found to be higher in the Mid Stage than in the Early Stage and Late Stage. All these results together indicate that those samples with lower OTU richness, which
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Fungal diversity measured along the chronosequence. (a) OTU richness and (b) Shannon index were used. To assess the relationship of years since deglaciation with diversity/richness indices of the total fungal community, we fitted three regression models (linear, quadratic, and cubic). Model selection was carried out based on adjusted R² and root mean square error. Only the best fit quadratic regression was shown.</figcaption></figure>

<table><caption>**Table 1.** R value and significance level (p) from 2-way ANOSIM test for differences in the fungal community structure between stages. ES, Early Stage; MS, Mid Stage; LS, Late Stage. Significance codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘’ 1.</caption><thead><tr><td>Groups</td><td>R statistic</td></tr></thead><tbody><tr><td>ES, MS</td><td>0.226**</td></tr><tr><td>ES, LS</td><td>0.581***</td></tr><tr><td>MS, LS</td><td>0.274*</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Soil fungal OTU richness re-ordered along the rank order of nestedness. Significance level is shown at ***p < 0.001.</figcaption></figure>

are mostly from Early Stage and Late Stage, tend to have an OTU composition that is a subset of Mid Stage communities.

**Physicochemical characteristics of the glacier foreland.** With age since deglaciation, concentrations of important nutrients for microbial activity such as TOC and TN increased steadily, while pH decreased (Fig. 3). The TOC in the soil increased significantly from a very low content of 0.03% in the newly exposed soils, which are close to the edge of glacier, to over 0.69% in the oldest foreland soils around 80 years old (R = 0.66, p < 0.0001), while the content of TN increased from 0.004% to 0.055% (R = 0.67, p < 0.0001). In contrast, the pH - which ranged from 7.8 to 9.4 in the whole sample set - decreased along the chronosequence (R = 0.38, p < 0.05). No significant time-related patterns were found in TP, C/N ratio, soil moisture, and sand or clay-silt proportion in relation to age since deglaciation. Cluster analysis of environmental variables and age since deglaciation in the retreating glacier area showed high correlation between TOC and TN (Spearman ρ² ≥ 0.6; Supplementary Fig. S3), suggesting TOC and TN are covariables.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Three environmental parameters are shown along the age since deglaciation of soil in the glacier of Midtre Lovénbreen: (a) TN; (b) TOC; (c) pH. TOC, total organic carbon; TN, total nitrogen. The values of TOC and TN are shown in percentages.</figcaption></figure>

CCA analysis was used to reveal which environmental properties best predicted variation in fungal communities in the deglaciated foreland area (Fig. 4). Together with two axes on the biplot, in an accumulative variance for the interaction between communities and variables, a total of 41.1% of variation was explained. Axis 1 explained 24.0% of the variation in the data, while axis 2 explained 17.1%. Axis 1 and 2 together separated out the three stages, with Early Stage soils associated with low values of CCA axis 1 and 2. TOC, pH, C/N ratio, soil moisture, the proportion of silt and clay, age since deglaciation, TN and TP together accounted for 30.6% of the total fungal community variation in combination. Among these selected environmental factors, TOC (pseudo-F = 2.2, p < 0.001) and deglaciated age (pseudo-F = 1.6, p < 0.005) were significant contributors to fungal community variability, and a forward test indicated that the most important factor was TOC.

## Discussion

We expected to see a directional replacement trend in community composition during succession, under the hypothesis that competition between species would begin to play an increasing role as nutrients accumulated within the initially exposed low-nutrient soil from the glacier, and as total population densities increased. The result showed a significant peak in OTU richness/diversity in the Mid Stage successional soils, and a progressive change in community composition with increasing age (Fig. 4). From this observation, together with the fact that there is nesting of the relatively OTU-poor Early and Late Stage communities into Mid Stage OTU-rich communities, it appears that a directional replacement model is the most consistent model for changes in the soil fungal community across the whole successional time span of this glacier foreland. By contrast, other studies in high resistance environments typical of the High Arctic have found directional non-replacement succession in arbuscular and root-associated fungal communities²⁵,²⁶. For example, Davey *et al*. studied a relatively restricted system - the plant roots of one particular species of plant (*Bistoria*), in a system dominated by mutualistic interaction between the fungus and host plant⁹. In this sense, it is not surprising that the root-associated fungal community
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Canonical Correspondence Analysis ordination plot of fungal community composition based on ITS1 gene OTUs and a vector overlay of the environmental variables. Environmental variables were shown in arrows and the significant ones were presented in red. Different colors of dots denote different age since deglaciation: green, 20–50 years; blue, 50–80 years; yellow, over 80 years. 1, age since deglaciation; 2, soil moisture; 3, C/N ratio; 4, proportion of silt and clay; TN, total nitrogen; TOC, total organic carbon; TP, total phosphorus.</figcaption></figure>

and the general soil fungal community may be following different trends. For example, the trade-offs for the host plant in terms of accumulating a diverse EcM community, as opposed to eliminating part of it, may need to be considered as part of the process²⁷. This consideration might go some way to explaining why EcM diversity instead showed turnover in another more fertile environment, an alpine glacier foreland⁹. Accumulating extra EcM species as they arrive in a relatively resource poor environment may improve the flexibility of nutrient capture by the host plant²⁷. In contrast, in a more nutrient-rich alpine glacier foreland, the same plant species may be more selective about its fungal community, eliminating some species which are no longer helpful in resource capture⁹. Additionally, from the perspective of fungal-fungal interaction, the competition between fungal species might be intense in a more nutrient rich environment (such as an alpine glacier foreland), which may also play a role in eliminating some species.

In this study, Ascomycota indicator OTUs occurred most frequently in the Early Stage of succession, while Basidiomycota indicator OTUs were more frequent in the Late Stage. Four of the six most abundant OTUs in the whole dataset were classified as Ascomycota and each showed a decreasing trend with deglaciated age. These results, and the change in the whole community composition (Fig. 4), suggest a replacement of Ascomycota by Basidiomycota along the chronosequence. Elsewhere, similar turnover in the fungal community, with Ascomycota replaced by Basidiomycota during succession, has been reported by Zumsteg *et al.*²¹ in the Damma Glacier foreland in the Central Alps. It is known that many groups of Ascomycota are abundant in extreme environments²⁸⁻³⁰ for example *Lecanoromycetes*, to which most of the lichen-forming fungal species belong²⁹, and *Dothidiomycetes* which is the ‘pioneer’ in unfavorable environments and enhances plant growth by improving nutrient and water acquisition²⁸. By contrast, Basidiomycota colonize more heavily vegetated soils²¹. The trend of Ascomycota being replaced by Basidiomycota along the succession therefore may imply: i) the accumulation of nutrients, ii) the maturation of the ecosystem. We were interested in finding what specific fungal subgroups of the Ascomycota and Basidiomycota gave rise to the turnover during succession, but unfortunately low resolution of our ITS1 sequences in the available databases results prevented interpretation at levels below phylum, in many cases. Lack of taxonomic resolution also hampered our identification of EcM fungi. However, the relative abundance of EcM sequences showed a linear regression with the deglaciated age (R = 0.43, p < 0.05, Supplementary Fig. 1) and was significantly lower in the early stage than middle stage and late stage (P < 0.05).

TOC, TN and pH, which showed high correlation to the deglaciated age in this study, have been considered important soil properties in shaping fungal communities²⁴,³¹. However, only TOC, besides the deglaciated age, was selected as a significant influencing variable by the forward selection in CCA analysis. TOC and TN were classified as covariables by a cluster analysis. As these two factors emerge as being so close in terms of predictive power, it is difficult to tell whether TOC or TN is actually more important. However, considering that fungi are heterotrophic and that some fungi can fix nitrogen³², it seems less unlikely that nitrogen limitation is the strongest structuring factor for the community because they are not wholly dependent on soil nitrogen. We therefore would suggest that TOC is a more important environmental variable in primary glacier succession.

In the succession we observed here, fungal OTU richness/diversity increased from the Early Stage towards the Mid Stage, and declined from the Mid Stage to the Late Stage. We suggest that the later decline in diversity might
partly occur because of increasing importance of steady interference competition (as opposed to very variable lottery competition²² on highly localized resource patches) over time as both fungal population densities and TOC and TN build up in the soil. The initial increase in diversity from the Early Stage to Mid Stage demands other explanations. In the Early Stage, diversity may be low because fewer viable niches are available in the very resource-poor soil, before significant build-up of TOC and TN. In these circumstances, the ability to colonize very patchy localized resources within the soil matrix, the ability to continue to grow slowly under low resource availability, and the physiological adaptations to higher soil pH, may all be necessary for establishment and survival. These would generally be regarded in ecology as conditions in which interspecific competitive interactions were weak and/or unpredictable, although the use of the term 'competition intensity' is slippery and elusive²³.

As resources in the soil accumulate with soil development³³,³⁴ and the pH becomes closer to neutral, the physiological specializations of the pioneer species may no longer be generally necessary, and a greater range of fungal species will be able to survive³⁵. Furthermore, there has now been time for groups that are slow dispersers but strong competitors to arrive³⁶. Nevertheless, we suggest that by the Mid Stage the resource concentrations are still low enough, or patchy enough, that many of the earlier pioneer species can still coexist along with the new arrivals (as is shown by the nestedness analysis, where the 'Early Stage' community is shown to be contained within the 'Mid Stage community'). A key factor may be that the fine scale mosaic of environments to some extent shields the pioneer fungal species against intense interference competition. This would be analogous to the situation in plant communities many mid successional environments, where 'pioneer' plant species coexist for a time in a fine scale mosaic of low soil nutrient conditions alongside faster-growing later successional species that require higher nutrient concentrations³⁷.

By the Late Stage, as a result of several decades of plant growth and decomposition, high resource concentrations are predictable enough throughout the soil that population densities of fungi are higher, with intense and stable interference competition and less opportunity for 'lottery' colonization of localized resources within the soil matrix. Fungi that can most effectively exploit these more abundant, stable resources are able to dominate, at the expense of the 'pioneer' fungi that had persisted into the Mid Stage. The loss of these early-stage fungi causes a decrease in overall diversity. This would be analogous to the displacement of 'stress tolerant' pioneer vascular plants in primary successional systems as conditions are ameliorated¹¹,³⁸.

A second and slightly different explanation may be that rather than the changing intensity of 'competition' being important in producing the pattern of community change and diversity, there is simply a distinct set of 'early' successional niches and another set of 'late' successional niches³⁹,⁴⁰. The 'mid-stage' diversity maximum may thus reflect a stage of intermediate soil development, pH and nutrient concentrations, in which both 'early' and 'late' niches are viable in competition. The result of this overlap is greater diversity. As soil develops beyond the niche limits of early successional species, they are lost and only late successional species survive (under lower diversity).

It would be very interesting to study other successional systems along chronosequences, to understand how widespread the turnover patterns we observed here are in fungal ecology. Additionally, field-based or laboratory-based microcosm systems might be able to add far more detailed observations, and tests of mechanisms that can come from experimental manipulation of conditions. Combining the metagenetic approach used here with metagenetic sampling of small soil metazoans, bacteria and archaea, and also adding soil metagenomes, could potentially give far more detailed and sophisticated understanding of the mechanisms in community and ecosystem development of primary successional systems.

## Methods

**Sampling and site description.** Our sampling site was the foreland of the retreating Midtre Lovénbreen Glacier, Svalbard (Supplementary Table 2). The glacier began to retreat in the 1920s, leaving an area of moraine which now covers around 10 km². Progressive time stages of the retreat have been mapped based on high resolution aerial photography and satellite imagery⁴¹. Interpolation between these isochrones allows estimation of the age since deglaciation of each point in the foreland. Moreau *et al.*⁴¹ described the vegetation in this deglaciated area as strongly influenced by age since deglaciation, and further from successional equilibrium than tundra outside the glacial foreland zone. In the foreland, pioneer species including *Salix polaris* are found as a sparse vegetation in younger areas, while other species such as *Carex nardina* and *Polygonum viviparum* also occur in the older foreland soils.

We took three transects, running from the current snout of the glacier to the edge of the moraine that has formed since the 1920s (Supplementary Figs S4–S6). In total, 39 samples were collected from three transects from the glacier front to the edge of the foreland moraine to ensure a representative coverage of the foreland (15 samples were taken from transect 1, 15 from transect 2, and 9 from transect 3). Because of the differing distances covered during retreat in different parts of the foreland, Transect 1 was sampled every 50 meters, and the Transect 2 was sampled every 80 meters, reflecting the longer length of the middle part of the foreland due to more rapid glacial recession. Transect 3 was sampled near transect 2 to ensure that enough samples were collected from each soil stage. At each sampling site, five subsamples were taken within a 1 m² quadrat, taking a scoop of 100 g of soil at 0–10 cm depth at each corner and at the central point of the square. The five subsamples were pooled together to make one sample of 500 g for that 1 m² quadrat. Sampling was carried out over 3 days at the end of July and beginning of August 2014.

**DNA extraction and soil properties.** On the same day as sampling, soil DNA extraction and soil moisture measurement were carried out at the field station laboratory on Svalbard. The remainder of each soil sample was packed and stored at –20 °C. The soil DNA was extracted using the MoBio Power Soil DNA isolation kit (MoBio Laboratories, Carlsbad, CA, USA) following manufacturer's instructions, and random empty vials were
chosen and run to serve as controls. The isolated DNA was stored at -80 °C until the PCR stage. Soil texture, TOC, TN, total phosphorus (TP), and pH were measured at National Instrumentation Center for Environmental Management, Korea by following the standard protocol of the Soil Science Society of America (SSSA).

**PCR amplification and sequencing of ITS1 region.** Each sequenced sample is prepared according to the Illumina 16S Metagenomic Sequencing Library protocols to amplify the ITS1 region. The DNA quality is measured by PicoGreen and Nanodrop. Input gDNA (10 ng) is PCR amplified. The barcoded fusion primer sequences used for amplifications were as follows: ITS1-F (5'-CTTGGTCATTTAGAGGAAGTAA-3') and ITS2 (5'-GCTGCGTTCTTCATCGATGC-3'). The final purified product is then quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantificatoin kits for Illumina Sequecing platforms) and the quality is identified by the LabChip GX HT DNA High Sensitivity Kit. (PerkinElmer, Massachusetts, USA). And then the paired-end (2 × 300 bp) sequencing was performed by the Macrogen Inc. (Seoul, Korea) using the MiSeq™ platform (Illumina, San Diego, USA).

**Data analysis.** Eight samples from the glacial foreland area were excluded from further analysis owing to either insufficient DNA for sequencing, or low read numbers (less than 2000 reads) generated by sequencing. Six of these eight excluded samples were collected from soils near the snout of the glacier (which had a deglaciated age less than 20 years). The remaining 31 samples were divided into three groups in relation to the isochrones of three main moraines⁴¹. The samples collected from areas deglaciated between 20 years and 50 years ago were categorized as the Early Stage. Those more than 50 years and less than 80 years old were assigned to the Mid Stage, while those more than 80 years old were assigned to the Late Stage.

The sequence data obtained was processed in Mothur⁴². Sequences with any ambiguous bases, sequences with more than 8 homopolymers and sequences with lengths less than 200 bp were removed using the screen. seqs command in Mothur. Putative chimeric sequences were detected and removed via the Chimera Uchime algorithm contained within Mothur⁴³ in de novo mode, which first splits sequences into groups and then checks each sequence within a group using the more abundant groups as reference. Rare sequences (≤10 reads) were removed to avoid the risk of including spurious reads generated by sequencing errors⁴⁴. A final even subsampling depth of 5,200 reads was applied to each sample. Operational taxonomic units (OTUs) were assigned using the QIIME implementation of UCLUST⁴⁵, with a threshold of 97% pairwise identity. Taxonomic classification of each quality sequence was done against the UNITE version 6.0 database⁴⁶ using the classify command in Mothur at an 80% cutoff with 1000 iterations. EcM fungi were determined by matching taxonomy assignments with established EcM lineages as determined by recent phylogenetic and stable isotope data⁴⁷. The Miseq sequence data used in this study are deposited in the MG-RAST server⁴⁸ under project ID 15403 (http://metagenomics.anl.gov/metagenomics.cgi?page=MetagenomeProject&project=15403).

**Statistical analysis.** To assess the relationship between fungal richness/diversity and environmental variables, as well as with age since deglaciation, the richness of OTUs and other diversity indices were calculated using the Mothur platform⁴² for samples standardized to 5,200 reads (size decided by default, set to the size of smallest sample size), as diversity is directly correlated with the number of sequences collected. Bray-Curtis distance was calculated to analyze fungal community similarity. To reduce the contribution of highly abundant OTUs in relation to less abundant ones in the calculation of Bray-Curtis measure, abundance data of OTUs were square root transformed. The pairwise differences in community composition of the four groups were calculated by analysis of similarity (ANOSIM) in relation to Bray-Curtis distance. To assess the best fitting model of correlations between age since deglaciation and diversity/environmental variables, linear, polynomial (quadratic) and power (cubic) models were tried out using SigmaPlot v 10.0 (Systat Software, San Jose, CA). R version 2.15.2 and R package 'labdsv' was used to calculate the indicator value and the significance of OTUs⁴⁹. Only OTUs with greater than 10% frequency, an indicator value >0.3⁵⁰ and a p value <0.05 were selected as indicator OTUs. Cluster analysis of all measured environmental variables was carried out using the Varclus procedure in the Hmisc package in the R platform, version 2.15.2.

**Nestedness analysis and Canonical correspondence analysis.** The nestedness analysis was performed by BINMATNEST with default input parameters⁵¹ to test the nested structure of the fungal communities. The significance of nestedness was tested using default input parameters and null model 3 which calculates the p-value for row and column totals following Geel *et al.*⁵². The correlation of OTU richness and packed matrix order (indicating nestedness categorized from high to low) of each sample was evaluated by Spearman's rho test⁵³. A constrained analysis was conducted using CANOCO⁵⁴ to assess the effect of environmental variables on the fungal community. Canonical correspondence analysis (CCA)⁵⁵ was used since the length of the extracted gradient was larger than 3 SD units. Forward selection was used to select significant explanatory variables with 999 permutations and only significant (p < 0.05) variables were included in the models.

## References

1. Schimel, J. In *Ecological Studies: Analysis and Synthesis* Volume 113 (eds Stuart Chapin, F. III & Körner, C.) 239–254 (Springer, 1995).

2. Jumponen, A. Soil fungal community assembly in a primary successional glacier forefront ecosystem as inferred from rDNA sequence analyses. *New Phytol.* **158**, 569–578 (2003).

3. Zumsteg, A., Bäath, E., Stierli, B., Zeyer, J. & Frey, B. Bacterial and fungal community responses to reciprocal soil transfer along a temperature and soil moisture gradient in a glacier forefield. *Soil Biol. Biochem.* **61**, 121–132 (2013).

4. Kaufmann, R. & C. Raffl. In *Global Mountain Biodiversity: a Global Assessment* (eds Körner, C. & Spehn, E.) 177–190 (London Parthenon, 2002).

5. Matthews, J. A. In *Ecosystems of the world* (ed. Lawrence, W.) Ch. 2, 17–38 (Elsevier, 1999).

6. Robbins, J. A. & Matthews, J. A. Regional variation in successional trajectories and rates of vegetation change on glacier forelands in south-central Norway. *Arct. Antarct. Alp. Res.* **42**, 351–361 (2010).
7. Matthews, J. A. *The Ecology of Recently-deglaciated Terrain: A Geocological Approach to Glacier Forelands* (Cambridge University Press, 1992).

8. Svoboda, J. & Henry, G. Succession in marginal arctic environments. *Arct. Alp. Res.* 37–384 (1987).

9. Davey, M. *et al.* Primary succession of Bistorta vivipara (L.) Delabre (Polygonaceae) root-associated fungi mirrors plant succession in two glacial chronosequences. *Environ. Microbial.* 17(8), 2777–2790 (2015).

10. Vater, A. E. & Matthews, J. A. Testing the ‘addition and persistence model’ of invertebrate succession in a subalpine glacier-foreland chronosequence: Fåbergstølsbreen, southern Norway. *Holocene* 23, 1151–1162 (2013).

11. Walker, L. R. & Del Moral, R. *Primary succession and ecosystem rehabilitation* (Cambridge University Press, 2003).

12. Blaalid, R. *et al.* Changes in the root-associated fungal communities along a primary succession gradient analysed by 454 pyrosequencing. *Mol. Ecol.* 21, 1897–1908 (2012).

13. Brown, S. P. & Jumpponen, A. Contrasting primary successional trajectories of fungi and bacteria in retreating glacier soils. *Mol. Ecol.* 23(2), 481–497 (2014).

14. Brown, S. P. & Jumpponen, A. Phylogenetic diversity analyses reveal disparity between fungal and bacterial communities during microbial primary succession. *Soil Biol. Biochem.* 89, 52–60 (2015).

15. Mühlmann, O. & Peintner, U. Ectomycorrhiza of Kobresia myosuroides at a primary successional glacier forefront. *Mycorrhiza* 18, 355–362 (2008).

16. Mundra, S. *et al.* Arctic fungal communities associated with roots of Bistorta vivipara do not respond to the same fine-scale edaphic gradients as the aboveground vegetation. *New Phytol.* 205, 1587–1597 (2015).

17. Lauber, C. L., Strickland, M. S., Bradford, M. A. & Fierer, N. The influence of soil properties on the structure of bacterial and fungal communities across land-use types. *Soil Biol. Biochem.* 40, 2407–2415 (2008).

18. Nielsen, U. N., Osler, G. H., Campbell, C. D., Burslem, D. F. & van der Wal, R. The influence of vegetation type, soil properties and precipitation on the composition of soil mite and microbial communities at the landscape scale. *J. Biogeogr.* 37, 1317–1328 (2010).

19. Smith, A., Marin-Spiotta, E., de Graaff, M. & Balser, T. Microbial community structure varies across soil organic matter aggregate pools during tropical land cover change. *Soil Biol. Biochem.* 77, 292–303 (2014).

20. Yarwood, S., Wick, A., Williams, M. & Daniels, W. L. Parent Material and Vegetation Influence Soil Microbial Community Structure Following 30-Years of Rock Weathering and Pedogenesis. *Microb. Ecol.* 69, 383–394 (2015).

21. Zumsteg, A. *et al.* Bacterial, archaeal and fungal succession in the forefield of a receding glacier. *Microb. Ecol.* 63, 552–564 (2012).

22. Kennedy, P. Ectomycorrhizal fungi and interspecific competition: species interactions, community structure, coexistence mechanisms, and future research directions. *New Phytol.* 187(4), 895–910 (2010).

23. Welden, C. W. & Slauson, W. L. The intensity of competition versus its importance: an overlooked distinction and some implications. *Q. Rev. Biol.* 61(1), 23–44 (1986).

24. Thormann, M. N., Currah, R. S. & Bayley, S. E. Succession of microfungal assemblages in decomposing peatland plants. *Plant Soil* 250, 323–333 (2003).

25. Nara, K., Nakaya, H., Wu, B., Zhou, Z. & Hogetsu, T. Underground primary succession of ectomycorrhizal fungi in a volcanic desert on Mount Fuji. *New Phytol.* 159, 743–756 (2003).

26. Wu, B., Hogetsu, T., Isobe, K. & Ishii, R. Community structure of arbuscular mycorrhizal fungi in a primary successional volcanic desert on the southeast slope of Mount Fuji. *Mycorrhiza* 17, 495–506 (2007).

27. Pohjanen, J. *et al.* Interaction with ectomycorrhizal fungi and endophytic Methylobacterium affects nutrient uptake and growth of pine seedlings *in vitro*. *Tree physiol.* 34, 993–1005 (2014).

28. Jumpponen, A. R. I. & Trappe, J. M. Dark septate endophytes: a review of facultative biotrophic root-colonizing fungi. *New Phytol.* 140(2), 295–310 (1998).

29. Miadlikowska, J. *et al.* New insights into classification and evolution of the Lecanoromycetes (Pezizomycotina, Ascomycota) from phylogenetic analyses of three ribosomal RNA-and two protein-coding genes. *Mycol.* 98(6), 1088–1103 (2006).

30. Selbmann, L. *et al.* Biodiversity, evolution and adaptation of fungi in extreme environments. *Plant Biosyst.* 147, 237–246 (2013).

31. Rousk, J. *et al.* Soil bacterial and fungal communities across a pH gradient in an arable soil. *ISME J.* 4, 1340–1351 (2010).

32. Miransari, M. Arbuscular mycorrhizal fungi and nitrogen uptake. *Arch. Microbial.* 193, 77–81 (2011).

33. He, L. & Tang, Y. Soil development along primary succession sequences on moraines of Hailuogou Glacier, Gongga Mountain, Sichuan, China. *Catena* 72, 259–269 (2008).

34. Kabala, C. & Zapart, J. Initial soil development and carbon accumulation on moraines of the rapidly retreating Werenskiold Glacier, SW Spitsbergen, Svalbard archipelago. *Geoderma* 175, 9–20 (2012).

35. Fujiyoshi, M. *et al.* Successional changes in ectomycorrhizal fungi associated with the polar willow Salix polaris in a deglaciated area in the High Arctic, Svalbard. *Polar Biol.* 34, 667–673 (2011).

36. Tilman, D. In *Praktische Zahnmedizin Odonto-Stomatologie Pratique Practical Dental Medicine* Vol. 99 (eds Schulze, E., Mooney, H.) 327–344 (Springer Berlin Heidelberg, 1994).

37. Jones, G. A. & Henry, G. H. Primary plant succession on recently deglaciated terrain in the Canadian High Arctic. *J. Biogeogr.* 30(2), 277–296 (2003).

38. Grime, J. Primary strategies in plants. *Transactions of the Botanical Society of Edinburgh* 43(2), 151–160 (1979).

39. McKane, R. B. *et al.* Resource-based niches provide a basis for plant species diversity and dominance in arctic tundra. *Nature* 415, 68–71 (2002).

40. Taylor, D. L. *et al.* A first comprehensive census of fungi in soil reveals both hyperdiversity and fine-scale niche partitioning. *Ecol. Monogr.* 84, 3–20 (2014).

41. Moreau, M., Laffly, D., Joly, D. & Brossard, T. Analysis of plant colonization on an arctic moraine since the end of the Little Ice Age using remotely sensed data and a Bayesian approach. *Remote Sens. Environ* 99, 244–253 (2005).

42. Schloss, P. D. *et al.* Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. *Appl. Environ. Microb.* 75, 7537–7541 (2009).

43. Edgar, R. C., Haas, B. J., Clemente, J. C., Quince, C. & Knight, R. UCHIME improves sensitivity and speed of chimera detection. *Bioinformatics* 27, 2194–2200 (2011).

44. Brown, S. P. *et al.* Scraping the bottom of the barrel: are rare high throughput sequences artifacts? *Fungal Ecol.* 13, 221–225 (2015).

45. Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. *Bioinformatics* 26, 2460–2461 (2010).

46. Abarenkov, K. *et al.* The UNITE database for molecular identification of fungi–recent updates and future perspectives. *New Phytol.* 186, 281–285 (2010).

47. Tedersoo, L., May, T. W. & Smith, M. E. Ectomycorrhizal lifestyle in fungi: global diversity, distribution, and evolution of phylogenetic lineages. *Mycorrhiza* 20, 217–263 (2010).

48. Meyer, F. *et al.* The metagenomics RAST server–a public resource for the automatic phylogenetic and functional analysis of metagenomes. *BMC bioinformatics* 9, 386 (2008).

49. Roberts, D. labdsv: Ordination and multivariate analysis for ecology. R package version 1 (2007) Available at: https://cran.r-project.org/web/packages/labdsv/labdsv.pdf. (Accessed: 17th April 2016).

50. Logares, R. *et al.* Biogeography of bacterial communities exposed to progressive long-term environmental change. *ISME J.* 7, 937–948 (2013).

51. Rodriguez-Gironés, M. A. & Santamaría, L. A new algorithm to calculate the nestedness temperature of presence–absence matrices. *J. Biogeogr.* 33, 924–935 (2006).

52. Geel, M., Ceustermans, A., Hemelrijck, W., Lievens, B. & Honnay, O. Decrease in diversity and changes in community composition of arbuscular mycorrhizal fungi in roots of apple trees with increasing orchard management intensity across a regional scale. *Mol. Ecol.* 24, 941–952 (2015).
53. McDonald, J. H. In *Handbook of biological statistics* 2nd edn (eds McDonald, J. H.) 221–224 (Sparky House Publishing, 2009).

54. Šmilauer, P. & Lepš, J. *Multivariate analysis of ecological data using CANOCO* 5 (Cambridge university press, 2014).

55. Ter Braak, C. J. Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis. *Ecology* 67, 1167–1179 (1986).

## Acknowledgements

This work was supported by a grant of the Polar Academic Program (PAP) funded by Korea Polar Research Institute (KOPRI) under the number 0409-20140091.

## Author Contributions

K.D. and W.K. contributed to field sampling. K.D., B.T., I.M. and N.L. did data analysis. K.D., H.C. and J.A. wrote the main manuscript text. All of the authors reviewed and contributed to the manuscript.

## Additional Information

Supplementary information accompanies this paper at http://www.nature.com/srep

Competing financial interests: The authors declare no competing financial interests.

**How to cite this article:** Dong, K. et al. Soil fungal community development in a high Arctic glacier foreland follows a directional replacement model, with a mid-successional diversity maximum. *Sci. Rep.* **6**, 26360; doi: 10.1038/srep26360 (2016).

<figure><img src="image_8.png" /></figure>

This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/