Multiple‐trophic patterns of primary succession following retreat of a high‐elevation glacier

Howmultiple, interacting components of complex soil communities assemble within forefields of receding glaciers is still largely unknown, especially at high-elevation sites (>5000 m a.s.l.) where plant succession is very slow. To examine succession of soil communities across different trophic levels, we investigated four major soil groups: bacteria, fungi, nematodes, and other non-fungal non-nematode microbial eukaryotes at the Puca Glacier in the Peruvian Andes spanning 9-, 24-, and 89-year-old deglaciated soils. This is the first study of microbial communities, other than bacteria, at a high-elevation chronosequence in the Andes Mountains. In addition, we characterized soil biogeochemical properties (e.g., C, N, moisture, and pH) and rates of microbial enzyme activities associated with C, N, and P acquisition. We found significantly correlated increases in estimated richness and high species turnover in all soil groups along the chronosequence. These shifts in soil communities were significantly correlated with microbial enzyme activities and measures of C, N, moisture, and pH. Stoichiometric comparisons of enzyme activities showed phosphorus (P) and carbon (C) limitation of microbial activity across the entire chronosequence with no hint of nitrogen (N) limitation. Taken together, the observed shifts in soil communities and biogeochemistry indicate coordinated increases in trophic complexity and ecosystem functioning during the initial 90 yr of microbial succession along the post-glacial chronosequence of the Puca Glacier.


INTRODUCTION
Records of natural glacial retreats are common and date back to the early 19th century. However, as humans exert more influence on climate warming, glaciers are receding at faster and faster rates (Marshall 2014). Although there have been many studies of the succession of bacterial communities following glacier retreat (reviewed in Nemergut et al. 2007, Bradley et al. 2014, Kim et al. 2017, how different components (e.g., microbial eukaryotes and microinvertebrates) of soil communities assemble within forefields of receding glaciers is still largely unstudied. For example, to our knowledge there have been no studies of trophic groups higher than bacteria at any receding glacier forefield in the entire Andes Range. Bacteria are generally better adapted to nutrient deficient environments than fungi and protists, because many bacteria can fix N 2 and CO 2 (via photoautotrophy and chemoautotrophy), and therefore, bacteria should dominate early successional communities (Schaaf et al. 2011. Given the wider metabolic breadth of the bacteria as a whole, they should be less tied to the appearance of plants during primary succession. However, successional trajectories of these early communities may also be influenced by local and regional factors. For instance, bacterial communities in the forefield of the Lyman Glacier in the North Cascade Mountains were linked to the presence of plants more than fungal communities (Brown and Jumpponen 2014), but the opposite pattern was observed at Bl aisen Glacier in Norway (Blaalid et al. 2012). Likewise, fungal but not bacterial components of the soil community at the Lyman Glacier chronosequence responded to soil age (Brown and Jumpponen 2014), while the opposite was found in the forefields of lava flows of Mt. Hekla in Iceland (Cutler et al. 2014). A recent study at the forefront of Hailuogou Glacier in China suggested that differences can be attributed to a combination of many factors including site specificity and random effects (Jiang et al. 2018). We suggest that another factor playing a role in community assembly following glacial retreat could include biotic interactions with soil organisms at higher trophic levels.
The deglaciated forefield of the Puca Glacier in the Peruvian Andes has been characterized as having high abundances of cyanobacteria (e.g., Nostoc, Pseudanabaena, and Leptolyngbya) (Schmidt et al. 2008). In the early stages (spanning 0-80 yr) of the chronosequence, the diversity and evenness of bacterial communities were positively correlated with the age of soils (Nemergut et al. 2007) as were their activities such as rates of N-fixation (Schmidt et al. 2008). Later studies at this site established that microbial and plant colonization was limited by nutrients (mainly P; Castle et al. 2017) and that soil nutrients were tied closely to microbial community succession (Knelman et al. 2014) suggesting that microbes and plants show similar successional patterns at this site.
Although previous studies at the Puca Glacier chronosequence provided insights into the succession of plant and bacterial communities, how communities of other soil microbiota such as fungi and microinvertebrates (e.g., nematodes) assemble along this post-glacial chronosequence is unknown. As many microinvertebrates feed on bacteria and fungi, the successional patterns of their communities could be tied to the patterns of their primary food sources. Nematodes are a major component of soil communities and are the most abundant and diverse microinvertebrates on Earth (Lambshead 2004, Hodda et al. 2009). Through their positioning as primary and intermediate consumers (Yeates et al. 1993, Bongers and, they play important roles in soil ecosystem processes such as decomposition of organic matter, mineralization of nutrients (Schratzberger et al. 2019), and plant productivity (Bernard et al. 2017). Hence, more detailed examination of nematode assemblages enables a better understanding of the complex interactions among members of the soil community.
To examine the process of succession of the entire soil community in this high-elevation glacial chronosequence, we collected replicated samples at three different spatial distances from the glacier representing 9-, 24-, and 89-year-old deglaciated soils. Our main objectives were to examine: (1) the diversity and composition of the main components of the entire soil community (e.g., bacteria, fungi, nematodes, and other nonfungal non-nematode microbial eukaryotes) along the chronosequence; (2) the potential functional roles of these groups along the chronosequence; and (3) the correlations among these soil groups and their relationship to edaphic factors. Given the relatively high pH of these soils which should be more conducive to bacterial than fungal life (Rousk et al. 2010, Zhang et al. 2016, and extremely low carbon and nitrogen content of these soils (Schmidt et al. 2008, Knelman et al. 2014, we hypothesized that succession of fungi, nematodes (especially fungivores and omnivores), and other eukaryotes would be slower than that of bacteria during the early stages of succession at this site.

Study site and sampling
The study was conducted at the forefield of Puca Glacier located in the Cordillera Vilcanota of the Peruvian Andes (13°46 0 24″ S, 71°04 0 17″ W, 5000 m a.s.l.; Knelman et al. 2014). The site receives about 100 cm of annual precipitation mostly as snow (Nemergut et al. 2007) and soil temperatures can oscillate between À10°C and over 25°C on a daily basis during the dry season (Schmidt and Vimercati 2019). Soil samples were collected on 23 March 2015 from the same v www.esajournals.org transects that were previously described by Knelman et al. (2014), which also correspond to transects 2, 3, and 4 as described by Castle et al. (2017). Briefly, we sampled three transects representing advancing stages of succession having been exposed by the retreating Puca Glacier for approximately 9, 24, and 89 yr. The 9-year-old location was identified by direct field observations of multiple expeditions between 2005 and 2012, and the age of older locations was based on aerial and satellite photographs taken between 1931 and 2003 and ground truthing the locations based on landmarks on the ground as described in more detail in previous publications (Nemergut et al. 2007, Schmidt et al. 2008, Knelman et al. 2014. In contrast to the presence of plants (25-50% plant cover) at the 89year-old site, there were no plants in the youngest soils and only biological soil crusts and some mosses at the 24-year-old site (Knelman et al. 2014). At each successional stage, four composite samples (consisting of 30-50 g subsamples) from the top 5 cm of soil within~400 cm 2 area were collected, placed into sterilized Ziplock bags, and gently homogenized. The four composite samples at each successional stage were at least 5 m from each other along the transect. Collection spoons were surface sterilized with ethanol wipes between samples to eliminate cross contamination. Samples were frozen overnight and shipped to the laboratory at the University of Colorado where there were immediately subsampled for specific biological and biochemical analyses as described below and subsequently processed or frozen till later.

DNA extraction and PCR
For the analysis of bacteria, fungi, and nonfungal non-nematode eukaryotes, total DNA was extracted from~0.3 g of soil subsamples using PowerSoil DNA Isolation kit according to the manufacturer's instructions. Because 0.3 g of soil is not enough to accurately represent nematodes, nematodes were extracted from~50 g of soil using a mobility-dependent method (Whitehead tray) over a period of 24 h ) and counted to a trophic group level (Yeates et al. 1993) under an inverted microscope. Extracted nematodes were reduced to 0.5 mL, transferred to PowerSoil bead-beating tubes, and then processed for DNA as described above.

Soil biogeochemistry and microbial activity
The detailed methods were described elsewhere . Briefly, soil moisture (%) was measured by drying 5 g of soil for 48 h at 60°C. Dissolved organic carbon (DOC), total dissolved nitrogen (TDN), and inorganic N (NH 4 + , NO 3 À ) were measured using soil extracts (5 g soil in 25 mL 0.5 mol/L K 2 SO 4 ). After mixing, soil slurries were centrifuged for 3 min at 20,124 g and immediately filtered through 0.3 µm glass fiber filters. A Shimadzu total organic carbon analyzer was used to measure DOC and TDN, and inorganic N (IN) was quantified with Lachat QuikChem 8500 Flow Injection Analyzer (Lachat Instruments, Loveland, Colorado, USA) and Synergy 2 Multi-Detection Microplate Reader (BioTek Instruments, Vinooski, Vermont, USA). Total N (%) was measured by combustion of~50 mg air-dried and ground soil using a Thermo Finnigan Flash EA 1112 Series CHN analyzer (Thermo Fisher Scientific, Waltham, Massachusetts, USA).
Microbial biomass C and N (µg/g dry soil) were estimated by treating 5 g of soil with 2 mL alcohol-free chloroform. After 24 h of incubation, the soil slurries were vented for 1 h to release chloroform and then measured with K 2 SO 4 extraction method as described above.
Activities (nmolÁh À1 Ág À1 dry soil) of seven extracellular enzymes involved in N (leucine aminopeptidase-LAP and N-acetylglucosamine -NAG), C (a-glucosidase-AG, b-glucosidase-BG, b-xylase-BXYL, cellobiosidase-CBH), and v www.esajournals.org P (phosphatase-PHOS) acquisition were measured in soil slurries made of 1 g of soil mixed together with 125 mL 50 mmol/L sodium bicarbonate buffer adjusted to pH 8.0 and homogenized at 3000 rpm for 1 min using Ultra-Turrax homogenizer. The sample slurries were transferred to 96-well plates with controls and substrates and incubated at 13°C for 22 h (Weintraub et al. 2007, King et al. 2010, Castle et al. 2017. The activity of enzymes was measured using Synergy HT Multi-Detection Microplate reader (BioTek, Winooski, Vermont, USA). Stoichiometric analyses using the seven enzymes were done to estimate which nutrients were most limiting to microbial activity at the time the soils were collected as discussed elsewhere (Hill et al. 2012. To derive enzyme ratios (C/N, N/P, and C/P) indicative of which soil nutrients are most limiting to microbes (Sinsabaugh et al. 2009), mean value of the sum of enzymes associated with C (AG+BG+CHB+BXLY) and N (LAP+NAG) were used.
Sequencing data processing QIIME 2 , Bolyen et al. 2019) was used to process all sequencing reads. First, the raw sequences were demultiplexed to each sample based on unique barcode sequence, and then, the forward and reverse sequences were joined to reduce sequence errors. The primer sequences of the 16S/18S/ITS region were further trimmed from the joined sequences. The average length of final reads for 16S and ITS was 250 and 130 bp for 18S. The dada2 (Callahan et al. 2016) pipeline was used to check chimeras and generate amplicon sequence variants (ASVs) tables at 100% similarity (Callahan et al. 2017). The taxonomy was assigned to the ASVs using blast (Camacho et al. 2009). The inhouse curated Silva_111 database (Pruesse et al. 2007) was used to assign taxonomy for bacteria and eukaryotes (including nematodes), while UNITE version8 (20190202) for fungi (Nilsson et al. 2019).

Community analyses
The ASV tables were filtered depending on the sample type and specific sequenced region. The archaea and mitochondrial ASVs were removed from the 16S ASV table, for the 18S non-fungal eukaryotic ASV table, fungi, nematode, and plant ASVs were removed and these are designated from now on as "other eukaryotes." For nematode 18S ASV and fungal ITS ASV tables, only nematode and only fungal ASVs were retained. The unassigned ASVs were removed from all of the ASV tables prior to analyses. The ASV tables generated in QIIME 2 were transferred into R V.3.5.3 (R Development Core Team 2019) for all downstream analyses. The number of reads of bacterial, fungal, and non-fungal eukaryotic ASVs was standardized to per 1 g dry soil, and of nematode ASVs to per 50 g dry soil. Alpha diversity indices (Chao1 and Shannon) were calculated using phyloseq package (McMurdie and Holmes 2013) and beta diversity (based on Bray-Curtis distance dissimilarity matrices) was calculated using vegan (Oksanen et al. 2011).
Because successional community development may result in changes of species functional traits, we assigned putative functional profiles to our fungal ASVs using FUNGuild (Nguyen et al. 2016) and to nematode ASVs using the Outline for Soil Ecologists (Yeates et al. 1993).

Statistical analyses
Venn diagrams of the unique and standardized ASVs for each soil group were generated using R package eulerr (Larsson 2020). The presence of significant differences in alpha diversity indices (Chao1 and Shannon) among communities from soils of different successional stages was tested using ANOVA with age/distance along the chronosequence as the main factor followed by Fisher's least significant difference (LSD). The presence of significant differences in community composition (beta diversity) was tested with PERMANOVA with 10,000 permutations in R package vegan (Oksanen et al. 2011). Because these tests were based on unrarefied data (differences in abundance and diversity likely reflect real biology and ecology in this system), the tests were repeated using Aitchison distance, which has some advantages for compositional data (Aitchison et al. 2000). Details of the Aitchison distance and tests are included in the Appendix S1. To determine whether alpha diversity patterns were similar across all soil communities, we compared the correlation of richness estimator Chao1 among the different communities (e.g., bacteria vs. fungi, bacteria vs. v www.esajournals.org nematodes, etc.) using a model-II regression implemented in R package lmodel2 (Legendre and Oksanen 2018). Because we ran six comparisons, we adjusted the P values with the False Discovery Rate correction (Benjamini and Hochberg 1995) with the R package stats (R Development Core Team 2019). The abundance of most notable taxa (e.g., bacterial genera) and functional traits (e.g., fungal saprotrophs and bacterial-feeding nematodes) was tested with ANOVA as above. All values were considered significant at P < 0.05 unless otherwise stated.
Differences in soil chemistry along the chronosequence were tested with ANOVA as above. Relationship of soil communities with geochemical characteristics was tested using envfit function in R package vegan, and significant soil characteristics were reported via PCoA ordination plots at P < 0.05. GDM (Ferrier et al. 2007) was used to compare how bacterial, fungal, nematode, and other eukaryote communities responded to soil geochemical gradients. Predictor data for GDM were highly co-correlated and were clustered using k-means clustering of pairwise unsigned correlation distance (1 À |r|), using k = 5 which was determined using the elbow method. Clusters were (1) PHOS, (2) inorganic N, (3) microbial N and microbial C, (4) DON, DOC, and LAP, and (5) pH, total organic N, total organic C, AG, BG, BXYL, CBH, NAG, and successional age. GDM was run with PHOS, inorganic N, microbial C, DOC, and age as predictor variables, explaining compositional dissimilarity (Aitchison distance using the isometric log ratio transformation; Egozcue et al. 2003) of community data for each of the four biological groupings mentioned above. We used GDM v1.3.11 (Manion et al. 2016).
To better understand potential roles of specific taxa at each successional stage, we used the multipatt function in the indicspecies package (De Caceres and Legendre 2009).

Community diversity
There was a general pattern of increasing alpha diversity along the chronosequence for all soil groups. The bacterial community ASV table amounted to a total of 4843 ASVs with 48,510 total sequences. The total number of observed ASVs increased with increasing age of soils with only 207 ASVs being shared among all three stages of the successional gradient (Appendix S1: Fig. S1A). Diversity measures (Chao1 and Shannon) showed that bacterial communities were significantly more diverse in the 89-year-old soils than in the younger soils, and the communities in the 9-and 24-year-old soils were not significantly different from one another (Table 1).
The number of nematode ASVs was at least an order of magnitude lower than of all investigated soil groups with a total of 67 ASVs represented by 77,980 sequences. Only four ASVs were present in the 9-year-old soils, compared with 62 ASVs in 89-year-old soils, and only two nematode ASVs (both Aphelenchoidid fungal feeders) were shared among the three sites (Appendix S1: Fig. S1C). Both Chao1 and Shannon diversity measures showed that the most diverse nematode community was present in the 89-year-old soils (Table 1).
There were 941 other eukaryotic (non-fungal, non-nematode) ASVs with 219,106 sequences. Similar to the bacterial and nematode communities, the other eukaryotic community was characterized by the significantly higher richness in the 89-year-old soils (Table 1; Appendix S1: Fig. S1D). However, the Shannon diversity index was highest in the 24-year-old soils ( Table 1).
The Chao1 richness estimates correlated significantly between all pairs of communities (Fig. 1). The correlation between alpha diversities of fungi and nematodes was particularly strong (Fig. 1D), and the relationship between bacteria and non-fungal, non-nematode eukaryotes richness was a close second (Fig. 1C).

Community composition, structure, and function
The community composition and structure (beta diversity) were significantly different among the three locations along the chronosequence for all four soil community groups (Table 2). Successional age along the chronosequence explained~40% of the variance of the soil communities with communities from the same location closely clustering together and clearly v www.esajournals.org separating from communities of other locations ( Fig. 2A-D) and this pattern remained unchanged when examined using Aitchison distance (Appendix S1: Table S1, Fig. S2).
Proteobacteria, Bacteroidetes, and Acidobacteria were the most common and abundant phyla in all our samples (Fig. 3A). Interestingly, despite high soil pH (~8), the most abundant ASVs belonged to a presumptive photosynthetic Chloracidobacterium (Acidobacteria) distributed fairly uniformly across all successional stages. Among the other most abundant photosynthetic taxa were cyanobacteria (e.g., Microcoleus, Nostoc, and Phormidium), which unlike Chloracidobacterium, were generally more abundant in the youngest soils ( Fig. 3A; Appendix S1: Table S2). The most  abundant Proteobacteria in all soils were Sphingomonadales (e.g., Sphingomonas), but other Alphaproteobacteria, especially Rhizobiales (including some known to be associated with plants) increased with soil age (Appendix S1: Table S2). Likewise, the more abundant Burkhoderiales (Betaproteobacteria) and Pseudomonadales (Gammaproteobacteria) (Appendix S1: Table S2) had their highest abundance in the 89-year-old soils. Verrucomicrobia, represented predominantly  by Chthoniobacter ( Fig. 3A; Appendix S1: Table S2), generally exhibited a similar pattern, however within Chthoniobacter many individual ASVs deviated from that general pattern (e.g., some ASVs were more abundant in the youngest soils).
In the 9-year-old soils, only three fungal-feeding (Aphelenchoides) and one bacterial-feeding (Tridentulus) nematode ASVs were observed (Figs. 3C, 4B). In the 24-year-old soils, the nematode community was slightly more complex with a single occurrence of an omnivore (Mesodorylaimus), two bacterial-, and mostly fungal-feeding species (9 in total) being present. In contrast, the 89-year-old soils contained complex nematode communities with all known trophic groups represented (Fig. 4B). This diversity at the trophic level was mirrored by the species diversity (Fig. 3C). For example, predatory nematodes in the family Nygolaimidae and plant parasites in the family Anguinidae were only present in the oldest soils (Appendix S1: Table S2). Bacterial feeders in the oldest soils were dominated by Plectus, Teratocephalus, and Prismatolaimus. With the exception of Plectus, all bacterial-feeding taxa were present only in the oldest soils.
Within the other eukaryotic community, the abundance of metazoan ASVs (mostly rotifers and arachnids) were highest within soils of the intermediate age (Fig. 3B) but the pattern was not significant. Collembola, Tardigrada (Hybsibidae), and Turbellaria were virtually absent within the youngest soils (Appendix S1: Table S2). Interestingly, the pattern of nearly complete absence in younger soils was also common to all types of algae including golden (Chrysophyceae), green (Chlorophyceae and Ulvophyceae), and red (Rhodophyceae) algae (Appendix S1: Table S2).

Relationships with soil biogeochemistry
Nearly all biogeochemical measures (except for DON) were significantly different among sites along the chronosequence (Appendix S1: Table S4) and correlated with beta diversity of all soil groups ( Fig. 2A-D). Unlike soil pH with its highest values (8.5) in the youngest soils, other soil measures were positively correlated with the age of soils (Fig. 2). While bacterial and fungal communities positively correlated with both C and N soil biochemical measures ( Fig. 2A, B), the nematode and other eukaryotic communities were correlated only with N measures (Fig. 2C, D).
Community composition of bacteria, fungi, nematodes, and other eukaryotes was highly predictable using biogeochemical data based on GDM. GDM models explained 83.4% of community dissimilarity for nematodes, 83.7 for fungi, 84.7% for other eukaryotes, and 75.8% for bacteria (Fig. 5). The heights of GDM splines indicate their relative importance in models (Ferrier et al. 2007, Fitzpatrick andKeller 2015), meaning PHOS was the most important predictor of bacteria, and other eukaryotes when all other variables were accounted for in the model. PHOS was fairly important for fungi and nematodes as well. Microbial C was the strongest predictor of nematode and fungal community composition, although microbial C was strongly correlated with microbial N, so this result could relate to either. Successional age was not an important variable in models for nematodes, fungi, or other eukaryotes. This result does not necessarily mean that age is not correlated with eukaryote community composition, but instead that age does not explain any variability in community composition that is not explained by other variables in the model. Age, however, was marginally important for bacterial community composition.
Microbial biomass C and N and activity of microbial enzymes associated with C and N acquisition significantly increased along the chronosequence (Appendix S1: Table S4) whereas the activity of phosphatase did not change along the chronosequence. The mean activity of phosphatase (491 nmolÁh À1 Ág À1 dry soil) was three times higher than that of C-associated enzymes (138 nmolÁh À1 Ág À1 dry soil), and over an order of magnitude higher than that of N-associated enzymes (35 nmolÁh À1 Ág À1 dry soil). This stoichiometric comparison by using the C:N and N:P ratios (Hill et al. 2012) is visualized in Appendix S1: Fig. S4 suggesting strong soil P deficiency and some C deficiency along the entire chronosequence.

DISCUSSION
We investigated communities of major soil groups including bacteria, fungi, nematodes, and other non-fungal and non-nematode eukaryotes (e.g., protists, algae, and microarthropods) in the forefield of the Puca Glacier spanning 9-to 89year-old soils. Previous work at this intensively studied site focused on bacteria and plants in relation to the succession of soil biogeochemical processes (Nemergut et al. 2007, Schmidt et al. 2008, Knelman et al. 2014, Castle et al. 2016. The study reported here builds on these earlier studies to add a more complete picture of the entire soil microbial community and its interactions with soil biogeochemistry along the chronosequence. As with previous studies of this chronosequence (Nemergut et al. 2007, Schmidt et al. 2008, Knelman et al. 2014, this study shows the expected increases in bacterial diversity and biogeochemical functioning with increasing soil age. More importantly, the present study demonstrates for the first time that diversity of other soil microbial groups increased in parallel with the bacterial community (Table 1) and that these groups were similarly and predictably related to soil biogeochemical characteristics (Figs. 2,5). This indicates that even in this extreme, high-elevation (5000 m a.s.l.) environment, all of the various components of the community increased in a more-or-less concerted fashion. While this may not be surprising for more temperate soils from lower elevations, it is somewhat surprising for an extreme environment where we previously hypothesized that bacteria would have an advantage due to daily temperature extremes, high soil pH, and the paucity of nutrients (Nemergut et al. 2007). Recent work analyzing the broader microbial communities in other extreme environments, such as Antarctic cryoconite holes and high-elevation soils of the Atacama region, also point to bacteria not necessarily being the primary components of the microbial communities in extreme environments (Costello et al. 2009, Solon et al. 2018, Sommers et al. 2018. For example, Sommers et al. (2018) found a strikingly high correlation between bacterial and eukaryotic community assembly in cryoconite holes along an extreme nutrient-limitation gradient in the McMurdo Dry Valley of Antarctica.
Because previous studies of the bacterial community of the youngest soils of the Puca Glacier showed a predominantly photoautotrophic bacterial community and an abundance of bacteria in the Comamonadaceae (Nemergut et al. 2007, Schmidt et al. 2008, Knelman et al. 2014, we expected to see mostly bacterial-feeding nematodes in the youngest soils. However, the dominance of fungal-feeding nematodes (rather than bacterial feeders) in the youngest soils suggests Scatter plots show how well GDM fit to the data, with observed community distance plotted against GDM's predicted community distance, and 1:1 lines shown in black (R 2 is for 1:1 line). Beneath each scatter plot are GDM i-splines for variables in the model. These splines are partial regression fits, and the shape of the spline reflects the rate at which community composition changes (y-axis) with increasing dissimilarity in that variable (x-axis). All variables are scaled in this plot such that minimum dissimilarity observed is zero and maximum dissimilarity observed is 1. The maximum heights of i-splines reflect the magnitude of total compositional change associated with that variable, for example, the variable's relative importance to the model. PHOS (blue) was the most important variable for both (A) bacteria and (D) other eukaryotes (B) and was important for (B) fungi and (C) nematodes as well, especially at smaller differences in PHOS activity. that there may be a significant fungal biomass to support nematodes in these young soils. It was also surprising that fungal guild analyses which showed a high relative abundance of putatively lichenized fungi in the youngest soils. If these fungi are actually lichenized, they are probably associated with cyanobacteria since there were few algal sequences detected in the youngest soils. Any fungal-cyanobacterial association at this site still awaits a direct confirmation. However, based on their dominance, ascomycete Verrucariaceae and cyanobacterial Nostoc are the most likely candidates as both have been recognized as lichen partners, particularly in alkaline post-glacial soils (Blundon andDale 1990, Zhang et al. 2016). These findings show why an analysis of the whole microbial community is needed if we are to begin to understand the microbial ecology and biogeochemistry of early successional soils. It should be noted that there were no lichens nor any other macroscopic structures visible in the 9-year-old sites at the time of sampling. Obviously, more in depth work is needed (including an analysis of lichens and explicit quantification of bacterial and fungal biomass) to elucidate the role of fungi in the youngest soils of the Puca and other high-elevation chronosequences.
Despite the parallel patterns of microbial diversity increase along this chronosequence (Fig. 3), the magnitude of increases in diversity was higher for all of the eukaryotic groups compared with bacteria. That is, the rate at which richness increased along the chronosequence was much lower for bacteria (twofold) than the other groups studied (threefold for other eukaryotes, fourfold for fungi, and 10-fold for nematodes). In fact, nematode communities in 9-and even 24year-old soils supported predominantly fungalfeeding species (as discussed above). Likewise, other metazoans (e.g., tardigrades, springtails, and enchytraeids) were only detected in the 89year-old soils suggesting that assembly of soil communities at higher trophic levels is slow at this high-elevation site and that fully functioning higher trophic-level interactions are not occurring during at least the first 24 yr following glacial retreat at this site. This finding emphasizes the extreme nature of the Puca Glacier site compared with other chronosequences at which soil animals have been studied. For example, even at the high-latitude (78.9°N) proglacial foreland of the Midtre Lov enbreen Glacier in Svalbard, multiple detritivore/omnivore microarthropod taxa were already present in 2-year-old soils (Hodkinson et al. 2004).
The present study also adds to our understanding of how nutrient limitation may be controlling patterns of microbial diversity along cold-dry environments such as the Puca chronosequence. Indeed, our measurements of enzyme stoichiometries along this chronosequence indicate extreme P limitation and a moderate level of C limitation across the entire chronosequence (Appendix S1: Fig. S4). These new findings support previous soil microcosm studies and a 6-yr field-fertilization study that also indicated strong P limitation at this site ). In the present study, the activity of soil phosphatases (PHOS) was more than an order of magnitude higher than the combined activities of N-processing enzymes (NAG and LAP in Appendix S1: Table S2) pointing to P limitation as has been observed at other cold-dry sites , Bueno de Mesquita et al. 2019. PHOS also was the strongest predictor of community composition for both bacteria and other eukaryote communities and was a significant predictor of fungal and nematode communities as well (Fig. 5). Thus, it is likely that P limitation, rather than a harsh climate per se is the primary factor slowing the development of multitrophic diversity along the Puca chronosequence.
Despite indications of severe P limitation along the chronosequence, all measures of C and N increased along the chronosequence. In fact, the oldest soils had an order of magnitude higher values of C and N than the youngest soils. Microbial C and N exhibited a parallel pattern, indicating a buildup of microbial biomass as soils aged and supporting the ecological principle that more productive ecosystems support higher biodiversity (Waide et al. 1999, Geyer et al. 2017). This positive relationship in the system studied here appears to partially result from facilitation of microbial processes and taxa through the course of succession. For instance, in the youngest soils with limited C and N resources, the dominant microbes included autotrophic and N-fixing cyanobacteria and lichenized fungi, taxa that not only are able to exist under extreme oligotrophy but also add C and N to the system which allows colonization of the site by higher trophic groups such as nematodes. It is also likely that taxa at higher trophic levels (e.g., microbial grazers) contribute to a buildup of N pools since they require less N than their prey provides (Ferris et al. 1997, Chen and and therefore release N from bacterial and fungal biomass that can be used for algae and plants growth. This is also reflected in the molar C-to-N ratios of the microbial biomass that start out at 10.9 in the 9-year-old soils and decline to 8.4 and 8.9 in the 24-and 89-year-old soils, respectively. The C:N of the biomass in the 24and 89-year-old soils are in line with the global mean of microbial biomass in vegetated soils of 8.6 (Cleveland and Liptzin 2007), whereas the 9year-old soils show a slight excess of C in the biomass compared with N.
In addition to biogeochemistry and diversity, there were some significant shifts in apparent microbial functional groups along the gradient. For instance, cyanobacteria (e.g., Nostoc and Microcoleus) were dominant in the youngest soils but declined as soils aged. In contrast, many potentially N-fixing species within Rhizobiales known to be symbiotic with plants (e.g., Bradyrhizobium, Microvirga, and Rhizobium) increased in abundance with the age of soils. Among fungi, the most dominant taxa in the youngest soils were identified as lichenized fungi (mostly from ascomycete Verrucariaceae that accounted for more than 50% of the fungal community) and saprotrophs (e.g., Rhinocladiella and Tetracladium). In contrast, plant parasites were present only in older (especially 89-year-old) soils and endomycorrhizal fungi only in the oldest (89-yr) soils. Colonization of the youngest soils (9-year-old) was limited to fungal-feeding nematode species, whereas the oldest soils contained other feeding groups (e.g., a bacterialfeeding Rhabdolaimus and an omnivorous/ predatory Nygolaimus). The unexpected early predominance of fungal-feeding nematodes suggests a significant fungal biomass must reside in these soils. It appears that lichenized fungi could be the most likely candidate food source for fungal-feeding nematodes.
It is important to note that functional traits assigned to taxa may be imprecise in some cases. First, relatively short length of sequences can be insufficiently informative to uncover the exact identity of microbial species (Clarke et al. 2017, Ahmed et al. 2019. Second, identified taxa within all soil groups were rarely exact matches to known reference sequences. For instance, among a total of 68 nematode species only two produced 100% hits. Third, many taxa were simply unknown. For instance, a large component of fungal communities, regardless of their successional stage, was unidentifiable even at the class level not only with the use of databases specifically designed for fungi (e.g., UNITE), but also databases that are more inclusive and updated on more frequent basis (e.g., NCBI). And fourth, functional traits of most microbial taxa are largely unknown. For example, while FUNGuild database classifies Exophiala as animal pathogens, a closer analysis shows that our sequences are not related to any animal pathogens. Many Exophiala have been recently shown to be DSE (dark septate endophytes) of alpine plant roots, and the closest BLAST matches to our Exophiala sequences were all from studies of plant roots or from rhizosphere soils. For example, one of our Exophiala ASVs was 100% identical to endophytic and rhizosphere fungi from roots of Paeonia ludlowii in Tibet (MK983935) and to montane pines in Montenegro (Lazarevi c and Menkis 2018) and a 99% match to high-elevation root endophytes from Colorado (Bueno de Mesquita et al. 2018). Likewise, another Exophiala ASV was a 98% match to Colorado root endophytes , and another was a 97% match to phylosphere fungi of European Beech (Cordier et al. 2012). Also, many Exophiala are Black Yeasts and are among the most abundant fungal phylotypes in sequencing studies of the soils of the Dry Valleys of Antarctica (Dreesens et al. 2014), and in desert crusts of the western USA (Bates et al. 2006). Therefore, the Exophiala sequences found in our study are likely root endophytes or perhaps crust associated fungi, but they are unlikely to be related to animal pathogens as they are classified in FUNGuild.
Finally, while a minimal number of samples were used for inference (four replicates at three distances for a total of 12 samples), they were sufficient to illustrate the key patterns reported here. Each of the 12 collected samples represented three homogenized subsamples capturing potential spatial heterogeneity on the scale of v www.esajournals.org less than one meter. Moreover, the four replicate composite samples representative of each soil age category were sampled to minimize autocorrelation and maximize the integration of spatial variation (King et al. 2010). Despite inherent differences in natural variation among investigated biotic groups, the observed diversity patterns for all groups were consistently similar, with statistical results producing a robust and coherent story. And finally, the results from this research were very much in agreement with results from previous research from the Puca Glacier chronosequence. In all, despite the shallow sampling, the analytical approaches used in this study provide for high confidence in the validity of the results.

CONCLUSIONS
We investigated all major soil groups including bacteria, fungi, nematodes, and other eukaryotes along the Puca Glacier chronosequence spanning 9-to 89-year-old soils. Here, we show that contrary to our expectations: (1) All of the microbial soil groups showed similar diversity patterns along this chronosequence with all groups becoming more diverse and more complex as soil became older, (2) the fungal component of the community was likely well established early in succession as indicated by the dominance of fungal-feeding (not bacterialfeeding) nematodes in the youngest soils, and (3) all soil groups were highly correlated with soil biogeochemical characteristics indicating that at these early stages of primary succession all of the various soil groups are really members of the same developing community.

ACKNOWLEDGMENTS
We thank Preston Sowell, Kelsey Reider, and Anton Seimon for their field and logistical support. This work was supported by NSF grant DEB-1258160 to Steve K Schmidt and the late Diana Nemergut. Travel and fieldwork were supported by a grant from the National Geographic Society Committee for Research and Exploration. SKS and JLD conceived the study and collected the samples, DLP carried out sample processing for all measures, DLP, PS, JLD, SKS, and WH analyzed the data, and DLP, SKS, and WH wrote the manuscript. All authors provided critical feedback and contributed to the final version of the manuscript. All authors have no conflict of interest.