High throughput sequencing combined with null model tests reveals specific plant‐fungi associations linked to seedling establishment and survival

Plant‐fungal interactions are important for plant community assembly, but quantifying these relationships remains challenging. High throughput sequencing of fungal communities allows us to identify plant‐fungal associations at a high level of resolution, but often fails to provide information on taxonomic and functional assignment of fungi. We transplanted seeds of Pinus cembra across an elevational gradient (1850–2250 m a.s.l.) and identified environmental factors and known fungal associates important for seedling establishment and survival. We then applied null model tests to identify taxonomically unassigned fungi associated with pine recruitment. Early seedling establishment was determined by abiotic environmental factors, while seedling survival was predominantly affected by biotic environmental factors (i.e., the abundance of a fungal pathogen known from literature and the distance to adult trees). Null model tests identified known mycorrhizal partners and a large number of unknown operational taxonomic units (OTUs) associated with seedling survival, including saprotrophic and pathogenic species. These results highlight that unknown fungal OTUs, which are usually discarded from analyses, could play a crucial role for plant survival. Synthesis. We conclude that high throughput metabarcoding paired with null model tests, is a valuable approach for identifying hidden plant‐fungal associations within large and complex DNA metabarcoding datasets. Such an approach can be an important tool in illuminating the black box of plant‐microbe interactions, and thus understanding ecosystem dynamics.

Virulence of the pathogens may depend on environmental factors (Bever et al., 2015;Packer & Clay, 2000), including soil nutrient availability and temperature (Alexander, 2010;Bever et al., 2015). Such effects can be further modified by mycorrhizal fungi, which can provide protection against disease, greater access to soil nutrients and resistance to climate extremes (Jung, Martinez-Medina, Lopez-Raez, & Pozo, 2012). It is also known that different aspects of plant-soil interactions affect different life stages of plants (van der Putten et al., 2013). In the early stages of recruitment, e.g. during germination, plants are especially vulnerable to abiotic factors (e.g. water availability), whilst seedlings may be more affected by the abundance of fungal pathogens and mutualists, as well as other biotic factors (such as competition with neighbouring plants or herbivory; Hersh et al., 2012;Packer & Clay, 2000;van der Heijden & Horton, 2009).
To date, most knowledge of plant-soil interactions has been obtained from lab and pot experiments, which have revealed the relative importance of positive and negative interactions between entire microbial communities associated with a particular plant species or plant community (Kardol, Deyn, Laliberté, Mariotte, & Hawkes, 2013). These experiments often only exhibit the overall net balance of positive and negative interactions (Klironomos, 2002;van der Putten et al., 2013), and the role of individual soil microbial species is not identified. Accordingly, conclusions from these studies remain broad and lack mechanistic detail. Further, experimental settings often over-simplify the complexity of real ecosystems, although it is known that the environmental context affects the outcome of plantsoil interactions (Kardol, Martijn Bezemer, & Putten, 2006;Manning, Morrison, Bonkowski, & Bardgett, 2008;van der Putten et al., 2013van der Putten et al., , 2016. The quantification of plant-soil interactions during different plant life-history stages across 'real world' environmental gradients has been challenging (Kardol et al., 2013). This is mainly due to methodological constraints, including difficulties in characterizing complex microbial communities, and difficulties in understanding the ecological roles of recorded species, given that only about 7% of an estimated 1.5 million fungal species are taxonomically and functionally assigned (Blackwell, 2011;Tedersoo et al., 2014).
High throughput metabarcoding has advanced our understanding of fungal diversity (Tedersoo et al., 2014), distribution , community composition (Mucha et al., 2017), ecological roles of members of the community (de Vries et al., 2018) and species interactions (Bogar et al., 2018). This method allows microbial taxa to be grouped into ecologically meaningful molecular operational taxonomic units (OTUs; Taberlet, Coissac, Hajibabaei, & Rieseberg, 2012), and to assign putative ecological functions such as saprotrophs, pathogens or mutualists (Nilsson et al., 2018). However, the identification of fungal diversity relies on curated databases (e.g. UNITE; Koljalg et al., 2005) and specific pipelines (e.g. FUNGuild; Nguyen et al., 2016), and the sparsity of data in these databases severely limits the taxonomic and functional assignment of taxa detected by metabarcoding. As a result, a large fraction of OTU-level diversity is often excluded from downstream ecological analyses (Blackwell, 2011;Tedersoo et al., 2014). For example, in a previous study on the spatial occurrences of plant-associated fungi (Merges, Bálint, Schmitt, Böhning-Gaese, & Neuschulz, 2018) we obtained a dataset of over 1,000 distinct fungal OTUs, indicating a highly diverse fungal community. However, BLAST searches of the acquired fungal OTUs against reference sequence databases (i.e. UNITE and GenBank; Koljalg et al., 2005) allowed us to assign only 58% of the OTUs to genus level. Only these taxa could be included in a literature research to identify candidate OTUs associated with our focal plant species, and the remaining 42% of unassgined OTUs were omitted from further analysis . This phenomenon has occurred in many similar studies where only a small fraction of the acquired OTUs were included in the final analyses (e.g. Glynou, Nam, Thines, & Maciá-Vicente, 2017;Mundra et al., 2015;Schmidt et al., 2017). It is unlikely that the problem of OTU identification due to incomplete databases will disappear soon given the sheer diversity of fungi (Blackwell, 2011;Tedersoo et al., 2014) and our limited biological knowledge of many taxa, especially those which are only known by internal transcribed spacer (ITS) sequences from environmental samples. Nevertheless, some unassigned OTUs could be major players in ecological dynamics (Toju et al., 2017;Toju, Yamamoto, Tanabe, Hayakawa, & Ishii, 2016), and analysing all community members, regardless of taxonomic or functional assignment status, could greatly facilitate our understanding of ecosystem dynamics.
Detecting influential species in a community, so-called keystone taxa, is important for understanding ecosystem function (Banerjee, Schlaeppi, & van der Heijden, 2018) and ecosystem responses to disturbance (Stinson et al., 2006). Different approaches have been proposed to pinpoint species with putative key ecological functions in natural communities including microbial network analysis (Banerjee et al., 2018;Berry & Widder, 2014). Network analysis has been applied to visualize co-occurrences between microbial species and statistically determine 'keystone' taxa (i.e. highly connected taxa, with a strong effect on the network structure (Banerjee et al., 2018;Berry & Widder, 2014). While network approaches identify possible interactions within soil fungal communities, we suggest a complementary null model approach, (Gotelli, Ulrich, & Maestre, 2011), which can detect whether the presence, absence or abundance of a 'species' is associated with a target variable, such as an ecosystem function or performance variable (Gotelli et al., 2011;Soliveres et al., 2016;Ulrich, 2010;Ulrich, Piwczyński, Maestre, & Gotelli, 2012). Null model randomization tests provide the opportunity to link a functional response to the abundance of species, functional groups or OTUs (Gotelli et al., 2011;Soliveres et al., 2016;Ulrich et al., 2012).
Further advantages are that (a), randomization tests are distribution free and well-suited to the non-normal distributions of OTU data; (b) potential species interactions are taken into account by maintaining the community structure during modelling (Gotelli et al., 2011;Ulrich et al., 2012); (c) null model randomization tests have a clear cause-and-effect hypothesis, where species affect the functional response. This underlying hypothesis is advantageous where a directional hypothesis is reasonable, e.g. for investigating effects of species on ecosystem functions (Gotelli et al., 2011;Soliveres et al., 2016). Null model approaches may allow the identification of species that influence other species or ecosystem functions (Gotelli et al., 2011;Soliveres et al., 2016;Ulrich, 2010), but their application to OTU-based microbial community data remains to be tested.
In this study, we used DNA metabarcoding to assess which fungal species were linked to the performance of Swiss stone pine (Pinus cembra) at two life stages, establishment and survival, across the pine's natural distribution. The elevational distribution of the pine, ranging from 1,850 to 2,250 m a.s.l., exhibits pronounced changes in both abiotic and biotic conditions, such as microbial community composition Neuschulz, Merges, Bollmann, Gugerli, & Böhning-Gaese, 2018;Roll-Hansen, 1989). Seedling establishment of P. cembra is affected by abiotic factors, such as light and water availability (Neuschulz, Mueller, Bollmann, Gugerli, & Böhning-Gaese, 2015), but also by biotic factors, such as competition or plant-soil interactions (Barbeito, Dawes, Rixen, Senn, & Bebi, 2012;Neuschulz et al., 2015). While P. cembra seedlings are thought to need several mycorrhizal species to efficiently acquire nutrients and water, seedlings are also severely affected by needle pests, such as the snow blight fungus Gremmenia infestans (Roll-Hansen, 1989;Smith & Read, 2008).
We hypothesized that seedlings experience different feedback effects from fungal communities during a) seedling establishment (0-4 months) than b) later stage survival (5-15 months). We predicted that initial seedling establishment is primarily vulnerable to abiotic environmental stresses (Kueppers et al., 2016;Walck, Hidayati, Dixon, Thompson, & Poschlod, 2011), while seedling survival will be more strongly affected by the availability of mycorrhizal partners and the presence of known fungal pathogens (Bardgett, Bowman, Kaufmann, & Schmidt, 2005). Furthermore, we predicted that a null model approach could help to identify fungal associates important for seedling establishment and survival that are so far taxonomically unassigned. We tested these hypotheses by conducting seed translocation experiments and DNA metabarcoding of fungal communities along replicated altitudinal gradients that span the entire elevational distribution of P. cembra and reach beyond the tree's upper distribution limit. We build this present study on that of Merges et al. (2018), which explored the occurrences of a P. cembra-associated subset of taxonomically assigned fungi in relation to environmental factors and adult conspecific trees. First, we used the full DNA metabarcoding dataset containing known plant associates  and tested the association of known plant associates with seedling establishment and survival of P. cembra. Second, we applied null model randomization tests to (a) reveal whether this approach can detect the same known fungal associates as used in the first analysis and (b) identify previously taxonomically unassigned fungi that are associated with seedling establishment and survival.

| Study area and design
We conducted this study in the Central Alps within the core of P. cembra's distributional range. We selected two elevational gradients

| Seed translocation experiments
To study the effects of plant-fungal associations on plant establishment and survival in an environmental context, we conducted seed translocation experiments in both valleys in 2014 and 2015. We divided each valley into nine elevational belts with 50 m elevational intervals ranging from 1,850 to 2,250 m a.s.l. Following a random stratified sampling design we installed ten seed bags in 2014 and 20 seed bags in 2015 at each elevational belt ( Figure 1). Each seed bag was made of 1.5 mm wire-mesh (to prevent loss of seeds) and contained five P. cembra seeds. Seed bags were placed 4 cm deep in the soil and fixed by metal pins. To break dormancy the seeds were placed in a wet clay-sand mixture and exposed to temperature shifts between 5 and 25°C for 22 weeks (simulated seasonal variation, G.  Figure   S1). We further monitored the survival of seedlings until the end of their second growing season in the following year ( Figure S1).

| Environmental factors
We tested the effects of environmental factors (i.e. light availability, temperature, soil moisture, vegetation cover and distance to conspecific adults) and fungal occurrences on the establishment and survival of P. cembra seedlings. We focused on these early life stages, since they are most vulnerable, and thereby the bottleneck of plant regeneration (Vitasse et al., 2012). Light availability was measured as canopy openness above each seed bag with a spherical densitometer. Ericaceous vegetation cover was recorded by estimating the percentage cover of dominant ground flora species: Loiseleuria procumbens (L.) Desv., Vaccinium spp. L. and Rhododendron ferrugineum L. within 1 m 2 of each seed bag (Braun-Blanquet, 1964).
Ericaceous plants compete for resources with P. cembra seedlings but also act as nurse plants by ameliorating harsh environmental conditions (Bardgett & Wardle, 2010;Castro, Zamora, Hódar, & Gómez, 2002). Distance to conspecific adults was measured as the distance from each seed bag to the closest adult P. cembra by estimation and using a laser range finder (Nikon 800S) for distances over 10 m. The distance to conspecific adults can be an important factor in plant recruitment, as species-specific pathogens and herbivores accumulate close to conspecific adult plants and negatively affect seedling survival (Janzen-Connell Hypothesis;Connell, 1971;Janzen, 1970).
Temperature was measured for 16 of the 30 seed bags per elevational belt per valley using iButton data loggers (Maxim). Soil surface temperatures were recorded every 4 hr over the duration of the study. We calculated the mean of daily maximum temperatures of the hottest three months (MMaxST), since extreme maximum temperatures can induce rapid drying of soils, and thereby induce desiccation of seedlings (Andrus, Harvey, Rodman, Hart, & Veblen, 2018;Kueppers et al., 2016;Tingstad, Olsen, Klanderud, Vandvik, & Ohlson, 2015). We also calculated the mean daily minimum temperatures of the coldest three months (MMinWT), since extreme frost events can lead to high seedling mortality in high-elevation ecosystems (Kueppers et al., 2016;Lenoir, Gégout, Marquet, Ruffray, & Brisse, 2008). We measured soil moisture under dry weather conditions in September by averaging five tensiometer (Theta-Kit version 3) measurements for each seed bag. Water availability is especially important during early establishment, where root systems are barely developed and desiccation can rapidly occur (Kueppers et al., 2016).

| Soil fungal communities
Pinus cembra is obligately mycorrhizal, i.e. it requires an ectomycorrhizal (ECM) mutualist for survival in field conditions (e.g. to acquire soil nutrients; Smith & Read, 2008). In contrast, the presence of a known pathogen, the snow blight fungus (G. infestans), severely limits P. cembra survival (Barbeito, Brücker, Rixen, & Bebi, 2013;Roll-Hansen, 1989), mostly by infecting needles covered by snow in winter. To assess the impact of these and other fungi on seedling P. cembra occurrence performance we re-analysed a DNA metabarcoding dataset of the soil fungal communities that have been recorded at each elevational belt . Occurrence data for pine-associated fungi obtained from this dataset are published in Merges et al., (2018).
This previous study focussed on occurrence patterns of pine-associated fungi in relation to environmental factors and their host, but the effect of fungal occurrences on seedling establishment and survival was not addressed . Soil fungal communities were recorded by collecting soil samples near the eight seed bags planted per elevational belt following the stratified-random sampling design, focusing on the two most distinct microhabitats (i.e. ericaceous vegetation cover and close to adult P. cembra individuals, Figure 1; Merges et al., 2018). Although soil microbial communities are known to be very stable across years, they often show high within-year seasonality (Lipson & Schmidt, 2004;Rudolph, Schleuning, & Piepenbring, 2018;Schadt, Martin, Lipson, & Schmidt, 2003). Therefore, soil was sampled twice (May and September 2015) around each seed bag, to account for seasonality, resulting in a total of 288 soil samples . ECM roots of P. cembra were collected to establish a reference database for P. cembra-relevant ECM fungi . Roots of saplings and adult trees were sampled at 1,850, 2,050 and 2,200 m a. s. l. in both valleys in May 2015, of which 100 ectomycorrhizal root tips per elevation were collected for DNA extraction. DNA extraction, amplification and sequencing are described in Merges et al. (2018).
In brief, 300 mg of each root and soil sample were used for DNA extraction (Cubero & Crespo, 2002). During polymerase chain reaction (PCR) the ITS2 region was amplified. An Illumina MiSeq was used for paired-end sequencing (2 × 300 bp) at Fasteris SA, Plan-les-Ouates, Switzerland . Default options for the Illumina pipeline developed by Bálint, Schmidt, Sharma, Thines, and Schmitt (2014) were applied. Fungi were identified by blasting the OTU representative sequences against UNITE database and GenBank nucleotide database using a 97% similarity threshold .
Rare OTUs from potential erroneous sequences were excluded, following the recommendation of Bokulich et al. (2013).
Where possible, the OTU reads were classified into pathogens and mutualists . For pathogens, a systematic search of peer-reviewed journal articles was conducted . For ECM, the P. cembra ECM root samples taken at the study site were used to identify a list of candidate ECM OTUs, which resulted in 35 species known to be ECM mutualists (Table S1, Merges et al., 2018). For further details on the classification of pathogens and mutualists of P. cembra see Merges et al. (2018).

| Statistical analyses
We fitted two models describing the determinants of establishment (i.e. seed germination and seedling establishment within the first growing season) and survival (i.e. survival from the end of the first growing season to the end of the second growing season). First, we used generalized linear mixed models with a binomial error distribution in the r package 'lme4' (Bates, Maechler, Bolker, & Walker, 2015). The predictor variables were environmental factors (i.e. light availability, temperature, soil moisture, vegetation cover and distance to conspecific adults), the diversity of P. cembra-associated mycorrhiza (OTU antilogarithm of the Shannon diversity) and the abundance of pathogenic fungal OTUs, as identified by literature research (see Merges et al., 2018). Abundance data were calculated at the plot level based on the number of soil cores in which an OTU was detected (i.e. between 0-8, Figure 1). We assume that OTUs with high plot abundances have a higher likelihood to interact with seedlings growing on the plots. Shannon diversity was chosen to represent the frequency of potential positive interactions since P.
cembra is known to be associated with a large number of mutualists, whereas abundance was chosen for pathogens, since there are only few fungal pathogens reported to be explicitly associated with P. cembra (Barbeito et al., 2012;Rainer et al., 2015). We also included plot ID, region and year as random effects to account for spatial and temporal autocorrelation. Observational level random factors were included in the seedling establishment model to account for overdispersion (Albrecht et al., 2015;Bates et al., 2015). We selected the models by adding explanatory variables in a stepwise manner based on a hypothesized 'hierarchy of controls' using the Akaike information criterion correcting for small sample sizes (AICc; Burnham & Anderson, 2002) and likelihood ratio deletion tests to find the most parsimonious models explaining either establishment or survival ( Figure S2; Table S1; Diaz et al., 2007;Manning et al., 2015). The first step consisted of a model only containing abiotic environmental factors (i.e. maximum summer temperature, minimum winter temperature, mean soil moisture, light availability), which we hypothesized to be the underlying ultimate cause for distribution of all subsequently added explanatory variables in the following steps.
In the second step, we added biotic factors, such as vegetation cover and distance to conspecific adults, since these proximate causes are potentially shaping the fungal community variables added in the following step. Finally, the third set of terms was fungal community data: Antilogarithm of Shannon diversity of ectomycorrhizal communities and the abundance of the two pathogenic OTUs (282,1,198).
We calculated marginal and conditional R 2 to explore the variance explained by the models (Nakagawa & Schielzeth, 2013). Marginal R 2 describes the variance explained by the fixed effects, whereas conditional R 2 is the variance explained by both fixed and random effects (Nakagawa & Schielzeth, 2013).
Second, to detect associations of potential mycorrhizal and pathogenic fungal OTUs with seedling establishment and survival that are not described in the existing literature, we applied null model randomization tests on the full OTU table (Gotelli et al., 2011;Soliveres et al., 2016). These tests are based on a null model approach where one linear regression is performed between the response and each given species (Soliveres et al., 2016;Ulrich et al., 2012). The observed regression slope is then compared to 1,000 random permutations of the species' values and a standardized effect size (SES) is calculated for each species according to: SES = S obs − S sim )/SD; where S obs is the observed regression slope , S sim is the average of the 1,000 simulated regression slopes and SD is the standardized deviation of the slopes of these 1,000 randomizations (Gotelli et al., 2011;Soliveres et al., 2016;Ulrich et al., 2012). Significant relationships between the response variable and the abundance of each species are assumed when SES values are higher or lower than an SES of 2 and −2 respectively (Gotelli et al., 2011;Soliveres et al., 2016;Ulrich et al., 2012). As the association between seedling performance and an OTU may be driven by shared responses to environmental drivers, we also conducted a more conservative test of association that accounted for the effects of environmental factors. To do this, we extracted the unexplained residual variances of models that included environmental factors and conducted a second round of null model randomization tests. Residuals were extracted from GLM models with the identical structure as the GLMM models to allow standardization of residuals with the 'rstudent' command in r to compensate for difference in leverage (i.e. obtaining studentized residuals; Kutner, Nachtsheim, Neter, & Li, 2004). Based on the increased risk of type II errors from multiple testing, we expected 5% of the OTUs to have a significant relationship with each response variable by chance (Gotelli & Ellison, 2013;Soliveres et al., 2016). All OTUs that were significantly associated with seedling establishment and survival were parsed against the fungal community database FUNguild (Nguyen et al., 2016). We only retained matches with over 97% similarity and highly probable confidence rankings to examine which ecological guilds OTUs were assigned to.

| RE SULTS
In the seed translocation experiments, 411 (15%) of the 2,700 planted seeds established as seedlings within the first growing season (i.e. the period between May and September in 2014 and 2015; Figure S1). Of these seedling cohorts 68 (17%) survived to the end of the following growing season (i.e. the period from September to September of the following year; Figure S1).
We identified a total of 1,074 OTUs across the whole study . Of these, we found two pathogenic OTUs (assigned to the snow blight fungus G. infestans) in 83 (29%) soil samples (Table S2). In P. cembra root samples, we found 35 ECM OTUs and 184 (64%) of the soil samples contained at least one of these (Table S2).
F I G U R E 2 Generalized linear mixed effects model showing the effect of mean daily maximum temperature of the hottest three months, soil moisture, light availability, fungal pathogen abundance operational taxonomic unit (OTU 282) abundance and the distance to an adult conspecific as well as their interacting effects on the probability of (a) seedling establishment and (b) survival of Swiss stone pine (Pinus cembra). Models were fitted assuming a binomial error distribution. Black points show the model fit (p < .01) with standard error added as black lines. Asterisks indicate the level of significance obtained from likelihood ratio deletion tests  Seedling establishment was negatively associated with maximum summer temperature and positively associated with light availability and soil moisture (Figure 2a, Table 1). However soil moisture was negatively associated with the establishment of seedlings at sites with high light availability (significant soil moisture x light availability interaction; Figure 2a, Table 1). In contrast, ectomycorrhizal diversity (i.e. antilogarithm of Shannon diversity) and pathogen abundance did not affect the establishment of seedlings, i.e. these variables were not selected for the most parsimonious model (i.e. Table 1). Seedling survival was significantly positively associated with light availability and light availability at long distances to conspecific adults (significant light availability x distance to conspecific adult interaction; Figure 2b, Table 2). In contrast to seedling establishment, survival of seedlings was negatively associated with the abundance of the pathogenic snow blight fungus, G. infestans, and the distance to conspecific adults (Figures 2b and 3, Table 2).
Survival of seedlings was estimated to be 16% lower when G. infestans (OTU 282) was present in one soil core per plot and 69% lower when present in six soil cores per plot (i.e. highest abundance) relative to seedling survival at plots without G. infestans (Figure 3).
For seedling establishment, the null model randomization test identified 29 (3%) of the 1,074 OTUs as significantly positively associated with established seedlings and 32 (3%) as negatively associated ( Figure 4, Table S2). Here, 82% of the OTUs were unassigned, 11% could be assigned to ectomycorrhizal fungi, 2% to (predominately) animal pathogen-saprotroph taxa, 3% could be assigned as lichenized and 2% as root endophyte. For seedling survival, the null model randomization test showed 296 (28%) of the 1,074 OTUs as significantly positively associated with survival of seedlings and 217 (20%) as negatively associated (Figure 4, Table S2). Here, 87% of the OTUs were unassigned, 10% were ectomycorrhizal fungi, 1% belonged to (predominately) animal pathogen-saprotroph taxa, 1% was assigned as lichenized and 1% as root endophyte. Twenty-two of 35 ECM OTUs (63%) present in our ECM reference database were identified to be significantly associated with the survival of P. cembra seedlings. For seedling establishment, the null model randomization test in which environmental factors had been accounted for revealed that, of the 1,074 OTUs, two (0.2%) were positively associated and eight (0.7%) were negatively associated with established seedlings. Here, 80% were unassigned, 10% were assigned as lichenized and 10% as root endophytes. (Figure 4, Table   S2). For seedling survival, the environment-controlled null model randomization test yielded 247 of the 1,074 OTUs (23%) as significantly positively associated with survival of seedlings and 185 as negatively associated (17%) (Figure 4, Table S2). Here, 91% of the OTUs were unassigned, 6% were ectomycorrhizal fungi, less than 1% belonged to (predominately) animal pathogen-saprotroph taxa, 2% were lichenized and less than 1% were root endophytes. Eight of 35 ECM OTUs (23%), which were present in our ECM reference database, were significantly associated with the survival of P. cembra seedlings.
The number of OTUs that were significantly associated with seedling establishment in the randomization tests was close to or lower than that expected by chance. The null expectation was 54

| D ISCUSS I ON
We found clear evidence that the very early life stages of juvenile trees (i.e. establishment of seedlings during the first four months after translocation) were determined by abiotic environmental factors, whereas the later stages of recruitment (i.e. survival of seedlings until the age of 15 months) were predominantly affected by biotic environmental factors (i.e., the abundance of a known fungal pathogen and the distance to adult trees). Using null model randomization tests, we revealed patterns of association between unassigned fungi and the establishment and survival of seedlings. This demonstrates the general potential of this method for identifying microbial species involved in the plant-microbe interactions, which drive plant community assembly.
Our results concur with previous studies demonstrating that young seedlings are especially vulnerable to abiotic factors, but that the later stages of recruitment depend more on the biotic environment (van der Heijden & Horton, 2009;Hersh et al., 2012;Packer & Clay, 2000). The early establishment of P. cembra was fostered by high light availability and soil moisture and limited by high maximum temperature. These findings are in accordance with a previous study on coniferous subalpine tree species in North American that showed similar positive responses to light availability (Kroiss, Hillerislambers, & D'Amato, 2015). Several studies have shown that warm and dry conditions during seedling establishment limit seedling survival rates (Andrus et al., 2018;Kueppers et al., 2016;Tingstad et al., 2015). Extreme temperature events are often related to summer drought stress, as the soil surface dries out rapidly and establishing seedlings with poorly developed root systems could fail to assimilate sufficient water for physiological processes (Andrus et al., 2018;Brodersen, Germinot, & Smith, 2018). The survival of P. cembra seedlings was reduced at sites with a high abundance of the pathogenic snow blight fungus G. infestans, and by increasing distance to conspecific adults. Snow blight fungi, such as G. infestans, are known to infect pine needles mostly under winter snow cover, where infected needles are killed (Burdon, Wennstrom, Ericson, Muller, & Morton, 1992;Roll-Hansen, 1989). Negative distance dependency of seedling survival is a well-known mechanism in the Janzen-Connell framework (Connell, 1971;Janzen, 1970) in which pathogens and herbivores accumulate close to adult plants, thereby creating negative above-and below-ground feedbacks for conspecific seedlings (Bell, Freckleton, & Lewis, 2006;Liang et al., 2016;Merges et al., 2018;Packer & Clay, 2000). However, in this study there was a decrease of seedling survival with increasing distances to conspecific adults. This could be explained by more favourable conditions for snow blight fungus G. infestans above the tree line (i.e. higher density snowpack), harsher abiotic conditions away from adult trees and missing positive interactions with the mycelia of suitable ectomycorrhizal partners.
We could not detect any effects of ground vegetation cover on the establishment and survival of seedlings, although several previous studies found significant relationships between these factors (e.g. Andrus et al., 2018;Bardgett & Wardle, 2010;Castro et al., 2002;Kueppers et al., 2016;Tingstad et al., 2015). Ectomycorrhizal diversity was also not associated with the establishment and survival of seedlings. This could be explained by the fact that mycorrhizal diversity was generally high, with a minimum of 16 ECM species being present at each elevational belt. This high ECM diversity could have provided ample inoculum, as well as a suitable mutualistic partner, F I G U R E 4 Operational taxonomic units (OTUs) identified by null model randomization tests as being significantly associated with (a) seedling establishment, (b) seedling establishment after accounting for environmental factors, (c) seedling survival and (d) seedling survival after accounting for environmental factors. Percentage of the total of OTUs identified presented at the top of each bar. Functional assignment was done with FUNguild (Nguyen et al., 2016). Taxonomic assignments and effect sizes are presented in Table S2   Apart from the ECM OTUs present in our reference database that were yielded in the null model tests as significantly associated with seedling survival, the majority of fungi identified were not present in the UNITE database (Koljalg et al., 2005) or the fungal community dataset FUNGuild (Nguyen et al., 2016). However, those few fungal species that were assignable belonged mostly to the ecological guild of mycorrhiza, thereby representing a group of fungi known to interact with plants (

| CON CLUS IONS
Identifying microbial species involved in plant-soil interactions is a major challenge in ecology, particularly in field-based studies.
Here we present one of the first field studies to identify a subset of functionally distinct members of the soil fungal community that are likely to affect plant establishment and survival. Our results show that the proportion of taxonomically and functionally unknown fungi that are associated with seedling performance is particularly high during early plant establishment, which may hint at previously unknown roles of fungi in seed germination and survival. The results presented here outline future research directions in above-and below-ground species interactions, for instance testing whether individual soil microbial species can accelerate or retard plant community succession (e.g., Fukami & Wardle, 2005;Wardle et al., 2004) or promote species coexistence (e.g., Connell, 1971;Hersh et al., 2012;Janzen, 1970;Packer & Clay, 2000). We demonstrate that DNA metabarcoding coupled with an ecologically relevant classification of ITS sequences and experimental linkage to plant life stages is a promising approach to unravel potential plant-fungal associations at a previously unattainable resolution.
Furthermore, this approach may also be applied to a wide range of other poorly understood soil taxa, which determine plant performance (e.g. oomycetes and bacteria). In the past, incomplete reference databases and limitations in taxonomically assigning fungal OTUs have led researchers to disregard a substantial part of the observed biodiversity. Thus, only a small fraction of plant-fungal interactions are known, and ecologically important relationships remain hidden. Our study reveals that the combination of community barcoding and a null model approach has the potential to overcome some of these constraints. For example, applying our approach in other systems could help to identify 'keystone' OTUs with an important role in structuring plant communities over a wide range of habitats and ecosystem types.

ACK N OWLED G EM ENTS
We are grateful for permission from landowners for working on their