Experimental evidence that root‐associated fungi improve plant growth at high altitude

Unravelling how species communities change along environmental gradients requires a dual understanding: the direct responses of the species to their abiotic surroundings and the indirect variation of these responses through biotic interactions. Here, we focus on the interactive relationships between plants and their symbiotic root‐associated fungi (RAF) along stressful abiotic gradients. We investigate whether variations in RAF community composition along altitudinal gradients influence plant growth at high altitudes, where both plants and fungi face harsher abiotic conditions. We established a translocation experiment between pairs of Bistorta vivipara populations across altitudinal gradients. To separate the impact of shifting fungal communities from the overall influence of changing abiotic conditions, we used a root barrier to prevent new colonization by RAF following translocation. To characterize the RAF communities, we applied DNA barcoding to the root samples. Through the utilization of joint species distribution modelling, we assessed the relationship between changes in plant functional traits resulting from experimental treatments and the corresponding changes in the RAF communities. Our findings indicate that RAF communities influence plant responses to stressful abiotic conditions. Plants translocated from low to high altitudes grew more when they were able to associate with the resident high‐altitude RAF compared to those plants that were not allowed to associate with the resident RAF. We conclude that interactions with RAF impact how plants respond to stressful abiotic conditions. Our results provide experimental support that interactions with RAF improve plant stress tolerance to altitudinal stressors such as colder temperatures and less nutrient availability.


| INTRODUC TI ON
Abiotic and biotic filtering are major processes determining the dynamics and structure of ecological communities.While a large body of literature has focused on disentangling the roles of these two processes in shaping ecological communities (Cadotte et al., 2015;Carmona et al., 2019;Kichenin et al., 2013), they do not act in isolation of each other.The environmental context may influence the outcome of biotic interactions (Maron et al., 2014;Pellissier et al., 2018).Along stressful abiotic gradients for example, the role of positive interactions in shaping ecological communities is accentuated (Callaway et al., 2002;He et al., 2013).In return, biotic interactions may influence how communities respond to environmental variation and ultimately drive local adaptation to abiotic conditions (Johnson et al., 2010;Runquist et al., 2020).Thus, to predict how species communities change along abiotic environmental gradients, we need to understand not only the direct responses of the species to their abiotic environment but also how these responses are indirectly modulated through biotic interactions.
Plants interact with a large diversity of microorganisms, including bacteria, fungi, protists, nematodes, and viruses, collectively known as the plant microbiota (Trivedi et al., 2020).One prime example of a highly interactive and diverse community involves plants and rootassociated fungi (RAF).Plant-RAF associations are mostly mutualistic and are of vital importance for the establishment and growth of both plants and fungi (Smith & Read, 2010).Furthermore, according to the Driver hypothesis (Hart et al., 2001), changes in root-associated fungal communities are primary drivers of plant community composition (Tedersoo et al., 2020).In symbiotic plant-RAF interaction networks, plants depend on fungi for resource allocation, defence against pathogens, and abiotic stress tolerance, whereas fungi gain direct access to nutrients and habitat (Rodriguez et al., 2009;Smith & Read, 2010).Symbiotic RAF include both endophytes and mycorrhizae.Endophytes are a guild that encompasses fungi inhabiting living plants without causing disease (Bacon & White, 2000), and mycorrhizal fungi exchange nutrients with plants through nutrient uptake from the soil (Smith & Read, 2010).Yet, plant-RAF interactions are not always positive, as mycorrhizal and endophytic fungi may also act as parasites (Rodriguez et al., 2009).
Plants from the same species may associate with different communities of fungal symbionts depending on the surrounding abiotic environmental conditions (Abrego, Huotari, et al., 2020;Cobian et al., 2019;Jarvis et al., 2015).Along altitudinal gradients, the composition of RAF communities changes within plant species (Abrego, Huotari, et al., 2020;Jarvis et al., 2015).Most studies have considered this variation to be due to the direct effects of the abiotic environment on the fungi, which follows the Habitat Hypothesis (Zobel & Öpik, 2014).For instance, higher altitudes are characterized by thinner soil depth, limited nutrient availability, colder temperatures, and stronger winds (Bahram et al., 2012;Jarvis et al., 2015;Matsuoka et al., 2016;Miyamoto et al., 2014).However, a less studied explanation is that, rather than the occurrences of fungi being directly influenced by the abiotic conditions, they could be shaped by the interactions between plants and fungi.In fact, laboratory experiments have shown that abiotic conditions can affect which fungal species plants associate with (Hoeksema et al., 2010), but there have been no field studies to confirm this pattern in harsh high-altitude regions.
Abiotic conditions at high altitudes result in smaller, more compact plants with lower reproductive fitness (De Villemereuil et al., 2018;Freschet et al., 2018).Under conditions of resource stress, plants tend to adjust their functional traits to alleviate stress levels by enhancing their capacity for resource uptake (Freschet et al., 2018).Furthermore, higher abiotic stress at high altitudes enhances positive interactions among plants, which enable the occurrences of stress-intolerant species by expanding their realized niche (Stress-Gradient Hypothesis; Brooker et al., 2008;Callaway et al., 2002;He et al., 2013).Likewise, it is increasingly recognized that symbionts alter species' niches (Afkhami et al., 2014;Brown & Vellend, 2014).In particular, a few observational studies have suggested that symbiotic fungal communities enable the establishment of plants in high-stress habitats such as high altitudes (Lynn et al., 2019;Pellissier et al., 2013;Rodriguez et al., 2008).
RAF communities are key components in ecosystem functioning, influencing plant growth by modifying the quality and flow of nutrients and water from soil to plants (Wardle et al., 2004), and contributing to soil carbon storage (Clemmensen et al., 2013).The roles of RAF in nitrogen cycling and plant nutrition are particularly important in arctic ecosystems, where nitrogen limitation is a common feature of soils (Hobbie & Hobbie, 2006;Robinson et al., 2020).Because arctic ecosystems are among the most climate-sensitive ecosystems on Earth, they are likely to experience substantial landscapelevel changes, resulting in the release of carbon stores from the soil, thus creating a positive feedback to climate change (Jansson & Hofmockel, 2020).The influence of plant-fungal symbiotic networks in controlling soil carbon and nitrogen acquisition emphasizes the need for comprehensive research in this area to inform climate change mitigation efforts.
While observational evidence suggests that associations with RAF communities are a key component of plant stress tolerance in high-stress situations such as higher altitudes (Lynn et al., 2019;Pellissier et al., 2013;Rodriguez et al., 2008;but see Wutkowska et al., 2021), experimental evidence is limited.Here, we experimentally investigated whether the changes in RAF community composition along arctic altitudinal gradients influence the growth of the associated plants at higher altitudes.To test this, we established a translocation experiment between pairs of higher-and lower-altitude populations of Bistorta vivipara (Alpine bistort) along altitudinal gradients in the Arctic.The RAF communities were characterized by applying DNA barcoding based on the ITS2 region to root samples.
To disentangle the influence of the changing fungal community from the general influence of changing abiotic conditions, we used a root barrier to exclude new colonization by RAF after translocation.While we expected abiotic variation to have a greater effect than fungal community composition on plant functional traits (i.e., we expected all plants to grow more at lower altitudes compared to higher altitudes, regardless of whether new fungal colonization was allowed or not), we expected RAF communities from higher altitudes to enhance growth of translocated plants compared to the growth of those plants not allowed to establish new fungal associations.Our main hypothesis was that the growth of the translocated plants is facilitated by allowing the plants to be colonized by those fungi that occur at higher altitudes.Among the fungal guilds, we expected those taxa known to establish ectomycorrhizal associations to be the main drivers of increased plant functional growth at higher altitudes.

| Study area and study design
The study area is located in the arctic-alpine tundra of Kilpisjärvi, northwestern Finland (69.0629°N, 20.8145° E).This area is characterized by monthly mean temperatures averaging −11.3°C in winter to 9.5°C in summer, an annual precipitation of 546 mm, and is largely represented by low tundra plant species.
We selected Bistorta vivipara as the focal study plant species.
B. vivipara is a common, long-living perennial plant in alpine and arctic environments of the Northern hemisphere.Due to its compact and relatively small root system, and the fact that it establishes ectomycorrhizal associations, B. vivipara is a model system to study root-associated microbial communities (Blaalid et al., 2014;Gardes & Dahlberg, 1996;Yao et al., 2013).Importantly to test our main hypothesis, the RAF communities of B. vivipara have been shown to distinctly change along altitude (Abrego, Huotari, et al., 2020).
In Kilpisjärvi, we selected as study sites three altitudinal gradients: one located in the northern side of the Saana fell, one located in the western side of the Jehkas fell, and one located in the southern side of the Jehkas fell (Figure 1b).Along each of the three altitudinal gradients, we determined the lowest and highest altitudes at which B. vivipara occurred, and then selected five sampling plots at least 100 m apart at both altitude extremes (30 locations in total, see Figure 1b; Table S1).The altitude of the sampling plots ranged from 557 to 632 m.a.s.l.(meters above sea level) for the lower-altitude plots, and 711 to 814 m.a.s.l. for the higher-altitude plots.We used satellitederived normalized difference vegetation index (NDVI), snow cover duration, air-and soil temperature, and moisture information to characterize and confirm the abiotic differences between the lower-and higher-altitude plots (Appendix S1: Figure S1).The lower altitude plots had on average higher NDVI (mean = 0.725 ± 0.015) than the higheraltitude plots (mean = 0.624 ± 0.017), experienced earlier snowmelt (June 3rd ± 1 day) than the higher-altitude plots (June 9th ± 1 day), and had higher average annual near-surface (+15 cm) and soil (−6 cm) F I G U R E 1 Experimental setup (a) and map of the study plots (b).In panel (a), A shows plants collected in the first year to account for species turnover, B plants collected in the second year with no treatment, C plants translocated from high to low or low to high altitude and excluded from resident fungi colonization, D plants translocated from high to high or low to low altitude and excluded from resident fungi colonization, E plants translocated from high to low or low to high altitude and allowed colonization from resident fungi, and F plants translocated from high to high or low to low altitude and allowed colonization from resident fungi.In panel (b), the map shows the exact locations of the plots, with black circles showing each of the three gradients.
The altitudinal translocation experiment was set up in June 2018.
Six B. vivipara individuals were randomly selected within each sampling plot (in total 180 individuals).From the six selected plant individuals in each plot, one was uprooted and its fine roots were collected for RAF fungal identification already in June 2018.The other five plant individuals were collected in June 2019, after they had been exposed to the following experimental set-up.One plant individual was left untreated as a general control, and the remaining four were translocated: two were excluded from colonization by resident RAF fungi and either translocated within their plot or between pairs of higher-and loweraltitude sampling locations within the same altitudinal gradient, and two were allowed to be colonized by resident RAF and either translocated within their plot or between pairs of higher-and lower-altitude sampling locations within the same altitudinal gradient (Figure 1a).For the translocation, 20 cm × 20 cm × 20 cm soil cubes were excavated around each plant individual.For the exclusion of the resident fungal community, we used a 1.5 mm thick polypropylene fabric with a random orientation of fibres and maximum apparent opening size (AOS) of 21 μm, which prevented the penetration of roots and most fungi but was permeable to water and nutrients.To avoid destruction due to reindeer grazing, the plant individuals were covered by a 30-40 cm tall wire net of 3 cm grid.
For all selected plant individuals, we measured the following functional traits in June 2018 before the experiment, as well 1 year after in June 2019 (except for those plants already collected in 2018): number of leaves, number of inflorescences, and maximum leaf length and width.Using size-related functional traits of plants, which serve as reliable indicators of plant performance (Freschet et al., 2018), is an effective method of inferring changes in function without the need to destructively measure the plant.Additionally, soil carbon and nitrogen (both measured with a LECO CHN-1000 elemental analyser, 1997) and pH (determined using a combination electrode in a 1:5 (v/v) suspension of soil in water) were measured at the plot level.For the latter, we pooled soil samples collected near (ca.30 cm) the focal B. vivipara individuals within each plot.Additionally, the fine roots of all specimens were immediately collected for RAF fungal identification.
The fine root samples were hand-cleaned for soil particles, first in the field and then more thoroughly in the laboratory.Within 2 days after the first root cleaning in the field, the roots were further cleaned in the laboratory once or twice, where the absence of soil particles was verified with the use of a magnifying lens.The root samples were stored by wrapping them in tissue paper and placing them in plastic bags containing moisture-indicating silica gel.Molecular analyses for samples from both 2018 and 2019 were initiated in early autumn 2019.

| DNA extraction and sequencing
DNA was extracted from root samples using plant DNA extraction protocol (Ivanova et al., 2008) with modifications that are described in the Supporting Information (Appendix S2).
PCR conditions in 12.5 μL reactions followed CCDB Platinum Taq protocol as described in Ovaskainen et al. (2020) with ITS3 and ITS4 primers (White et al., 1990) modified with 6 N heterogeneity spacers and Illumina adapters as described in Abrego, Roslin, et al. (2020) with synthetic control (consisting of 9 plasmids derived from Palmer et al. ( 2018) with modifications prepared and diluted as described in Ovaskainen et al. (2020)).Sequencing on Illumina MiSeq followed methodology described in Ovaskainen et al. (2020) and Abrego, Roslin, et al. (2020).

| Bioinformatic analyses
Bioinformatic processing was performed using a pre-release development version of the OptimOTU pipeline, which is implemented using the {targets} pipeline management package 0.11.0 (Landau, 2021) in R version 4.0.5 (R Core Team, 2022).The exact version of the pipeline used in this paper is available at https:// github.com/ brend anf/ bisto rta_ vivip ara_ trans location, and is also described in more detail in the Supporting Information (Appendix S2).The steps of processing included raw read filtering and trimming, denoising to form amplicon sequence variants (ASVs), ASV filtering and trimming, clustering to form operational taxonomic units (OTUs), removal of non-fungi, and taxonomic identification.
The first round of filtering and trimming was applied to raw demultiplexed R1 and R2 read pairs using Cutadapt 4.0 (Martin, 2011).
After trimming, read pairs where either read was less than 100 bp were also removed, as were pairs where either read contained ambiguous bases ('N's).A second quality filtering step was then performed using DADA2 version 1.18 (Callahan et al., 2016) to remove read pairs where R1 had more than 3 expected errors or R2 had more than 5 expected errors, and to remove matches to the PhiX genome.Following quality filtering, R1 and R2 reads were separately dereplicated, denoised, and merged according to the ITS variant of the DADA2 pipeline tutorial (Callahan et al., 2020), separately for each sequencing run.ASVs were then taxonomically identified using ProtaxFungi, resulting in taxonomic assignments of each ASV at ranks from phylum to species, along with a calibrated probability that each assignment was correct (Abarenkov et al., 2018).We chose a cutoff of 90% probability for reliable identification and discarded other identifications.
ASVs were then clustered using a 3-stage approximately singlelinkage clustering process at each rank from phylum to species.In the first stage, ASVs which were reliably identified as belonging to the same taxon were merged to form cluster cores.In the second stage, unidentified ASVs were matched to the nearest identified ASV within a rank-and taxon-specific threshold using the usearch_global command in VSEARCH, and joined to the relevant cluster.Finally, in the third stage, remaining unidentified ASVs were de novo singlelinkage clustered using BLASTCLUST version 2.2.26 (Dondoshansky & Wolf, 2000) based on pairwise distance matrices calculated by the calc_distmx command in USEARCH version 11.0.667(Edgar, 2010) which were processed into BLASTCLUST format using custom R code, and passed to BLASTCLUST using its -r option.The entire clustering process was performed at each rank successively, with clustering at lower ranks constrained by the results at higher ranks.Optimal rank-and taxon-specific clustering thresholds were determined using the same general method as described for DNAbarcoder (Vu et al., 2022), but using our hybrid USEARCH+BLASTCLUST single linkage clustering method rather than the DNAbarcoder software, and using identified sequences from the Global Spore Sampling Project (Ovaskainen et al., 2020), based on the same ITS3-ITS4 amplicon as in this study, as referenced.The end result of the clustering process was a hierarchical classification of all ASV sequences, where some taxa had reliable names based on Protax results, and other taxa were the result of de novo clustering.The latter taxa were assigned placeholder names of the form 'pseudo{rank}_NNNNN', for example 'pseudogenus_00012'.We chose taxa (whether identified or pseudo) at the species rank as the primary unit of analysis, and refer to them in the rest of the manuscript as OTUs.

| Statistical analyses
To evaluate how the changes in plant functional traits due to the experimental treatments were connected to the changes in the RAF communities, we analysed the data with the joint species distribution modelling (JSDM) framework of Hierarchical Modelling of Species Communities (HMSC) (Ovaskainen & Abrego, 2020;Ovaskainen et al., 2017).This modelling framework allowed us to simultaneously examine how the RAF communities and plant functional traits responded to measured environmental variables, and to the experimental treatments in particular.
We constructed three models.The first model (which we henceforth call the 'fungal model') focused on evaluating how the RAF communities responded to the measured environmental variables and experimental treatments.The second model (which we henceforth call the 'plant model') measured the responses of plant functional traits to the environmental variables and experimental treatments.
The third model (which we henceforth call the 'fungal-plant model') focused on measuring how RAF communities were linked to plant functional traits along the environmental gradient.For each of the three models, the response matrices varied as follows.

| Fungal model
We constructed a response data matrix where the matrix elements describe for each sample (i.e., plant individual) the number of sequencing reads assigned for each OTU-level fungal species.The full fungal data matrix consisted of 3010 fungal species-level OTUs across the 180 plant individuals.To obtain sufficient statistical power in the JSDM model, we included only common RAF species, here interpreted as those that occurred on at least 50 plants (across all sites and treatment types), resulting in 108 species (Table S3).To assess whether the 2902 rare species that were excluded showed different responses than the common species that were included, we jointly modelled the species richness of the common (i.e., those occurring in ≥50 plants) and of the rare species (i.e., those occurring in <50 plants) in a bivariate Poisson regression model, henceforth called the species richness model (Appendix S4).The results of these analyses showed that, at the level of species richness, the responses of the common species were consistent with those of the rare species (Table S4).Additionally, we removed one sample from the main fungal model which contained less than 10,000 fungal sequence read counts, resulting in 179 sampling units containing between 24,124 to 237,907 fungal read counts with a median of 116,454 read counts.Due to the zero-inflated nature of the data, we fitted to the fungal data a hurdle model consisting of two parts: presence-absence (modelled with probit regression) and abundance (measured as the number of reads) conditional on presence (modelled with linear regression, with absences declared as missing data, and counts log-transformed and then normalized to zero mean and unit variance).We fitted these two parts separately and included for both parts a phylogenetic correlation term to quantify to which extent the responses were similar for related species.The phylogenetic tree was approximated using the taxonomic classification of the OTUs, including pseudotaxa, with equal branch lengths for each Linnaean rank.

| Plant model
For the plant model, we constructed a response matrix which included B. vivipara's leaf width and length, presence or absence of flowers, and number of leaves.We also measured approximate leaf area (length × width), but it was highly correlated with leaf length, thus not included in the analysis.The leaf width and length were modelled with normal distribution, the presence-absence of flowers with probit regression, and the number of leaves with Poisson distribution.These variables were modelled jointly, as HMSC allows the application of different link functions for the columns in the response matrix.

| Fungal-plant model
For the fungal-plant model, we constructed a mega-response matrix which included all the elements from the previous two models simultaneously: presences-absences and abundances of fungi, as well as leaf width and length, the presence or absence of flowers and number of leaves of B. vivipara plants.
The core explanatory part of the models was the same for all three models.We controlled for the study design by including as random effects the altitudinal gradient (three levels), sampling plot (30 levels), and plot of origin (30 levels).These random effects aimed to capture spatial covariates that influenced the plant functional traits and RAF communities but were not explicitly accounted for in the models, such as microclimatic variation in snow melt affecting plant-fungal interaction phenology (Mundra et al., 2015).In the fungal and fungal-plant models, we further included individual sample as a random effect (179 levels) to estimate fungal-fungal and plantfungal associations.In the plant model (n = 330 sampling units), we further included the plant individual as a random effect (180 levels) to account for the fact that individuals that were sampled in the second year resulted in two data points.We parameterized the experimental treatment effects measuring how altitude and colonization by the resident RAF community affected plant growth traits and RAF communities through 12 variables (Table 1).We additionally controlled for the effects of the following covariates that were measured in the field: total carbon in the soil (we also measured nitrogen, but it correlated highly with carbon so was not included) and pH (Figure S5).For translocated individuals, these covariates were measured in the plot of origin.We further controlled for log-transformed sequencing depth.
We fitted all models with the R-package 'Hmsc' (Tikhonov et al., 2020) assuming the default prior distributions (see Ovaskainen & Abrego, 2020).We sampled the posterior distribution with four Markov chain Monte Carlo (MCMC) chains, each of which was run for 37,500 iterations, of which the first 12,500 were removed as burn-in.The chains were thinned by 100 to yield 250 posterior samples per chain and so 1000 posterior samples in total.We examined MCMC convergence by the potential scale reduction factors (Gelman & Rubin, 1992) of the model parameters (Appendix S7: Table S7).
The explanatory power of the joint species distribution model was assessed through AUC (Pearce & Ferrier, 2000) and Tjur's R 2 (Tjur, 2009) values with the presence-absence part of the model, and through R 2 with the abundance part of the model, and for the species richness model through a pseudo-R 2 .We applied variance partitioning to investigate the proportion of the explained variance in species occurrences, abundances, and species richness attributed to each fixed and random effect included in the models.We also evaluated the β-parameters describing species-level responses to the predictors.To examine if related species showed similar responses to the environmental predictors, we examined the level of phylogenetic signal in the species responses in the fungal model by estimating the parameter ρ.In particular, we examined which specific taxonomic groups showed systematic responses to the treatments.
To examine the statistical associations among fungi co-occurring in the same plant individuals and the statistical associations between the fungal species and plant growth, we estimated the Ω-parameters in the fungal and plant-fungal models, respectively.
We conducted all statistical analyses using R version 4.2.0 (R Core Team, 2022).

| RE SULTS
In total, 3010 fungal species-level OTUs were detected.The average number of fungal OTUs per plant individual was 150 ± 35.
Ascomycota was the most abundant and species-rich phylum found in the root samples, with 62.7% of all sequence reads and 57.6% of all species belonging to this group (Table S8).Helotiales was the most abundant Ascomycota order, accounting for 45.4% of all sequence reads.Basidiomycota was the second most abundant phylum found (36.1% of reads), with Russulales being the most abundant order in this group at 18.5% of all sequence reads (Table S8).

| How much variation does translocation explain compared to natural variation in RAF communities and plant functional traits?
Compared to the natural variation in RAF community composition and functional plant traits along the altitudinal gradient, the experimental treatments involving translocation explained little variation (Table 2).In all models, the variables measuring natural variation (covariate group N in Table 2) explained two to three times the variation explained by the treatments (covariate group T in Table 2).The variables measuring natural variation and the experimental treatments explained, respectively, 76.3% and 19.5% in the 'Fungal presenceabsence model', 65.9% and 30% in the 'Fungal abundance model', 74.2% and 25.8% in the 'Plant model', and 76% and 20.5% in the 'Fungal-plant model'.Out of the experimental treatments, the largest amount of variation was found when measuring to what extent the effect of within-plot translocation was different for plots located in the lower and higher altitudes.
Among the variables measuring natural variation, the identity of plant individual explained most of the variation in RAF communities (explaining 35.1% and 15.2% of the variation in species occurrences and abundances, respectively), followed by altitude (explaining 7.1% and 5.4% of the variation in species occurrences and abundances, respectively), and soil pH (explaining 7.1% and 8.9%, respectively).On the other hand, for plant functional traits, the identity of the gradient explained most natural variation (explaining 26.2% of the variation in plant functional traits), followed by altitude (9%), while pH explained little variation (2%).The fact that most of the variation was natural rather than treatmentrelated was also reflected by the fact that the plot of origin explained more variation than the plot of sampling for both RAF communities (12.3% and 13.1% vs. 6.8% and 10.5%, for species occurrences and abundances, respectively) and plant functional traits (14.9% vs. 8.6%).

| How do RAF communities respond to translocation across altitudes?
In line with the results from the variance partitioning, the results from the analyses evaluating the species-level responses to the predictors revealed that the experimental treatments involving translocations had generally mild effects on the occurrence and abundance of most fungal species (Figure 2).However, the occurrences and abundances of several species were statistically affected by the experimental treatments.Furthermore, these responses were phylogenetically structured, as confirmed by the high value of the ρ parameter (ρ = 1.00 with posterior probability ≥ 95%) measuring phylogenetic (taxonomic) signal (Table S7).Namely, experimental treatments consistently influenced the occurrences and abundances of species within higher taxonomic groups.
Within-plot translocation generally had a negative effect on RAF occurrences and abundances, except in the case of the Ascomycota genera Pezicula and Mollisia, which increased in abundance.The within plot translocations had more negative effects when being carried out at the higher altitudes than the lower altitudes.Fungal OTUs showed both negative and positive responses to exclusion of colonization by the resident community, both in their occurrences and their abundances.OTUs assigned to the genus Mortierella, for example, increased in abundance when translocations with exclusion were carried out at higher altitudes.Generally, species responded negatively to translocations from lower to higher altitudes, decreasing in abundance.In the case of translocation with exclusion, RAF communities responded negatively in occurrence but positively in abundance, especially when translocated from higher to lower altitudes.This was especially the case for OTUs assigned to the orders Helotiales, Chaetothyriales, Dothideales, and Pleosporales, including OTUs that could be assigned to the genera Phialocephala, Mollisia, Pezicula, Cladophialophora, Sporormiella, Phaeosphaeria, Herpotrichia, and Leptosphaeria, and OTUs that could be assigned to the families Leotiaceae, Didymellaceae, and Dothioraceae.Among the other predictors, both altitude and pH influenced RAF species occurrences and abundances the most, with many OTU's occurrences and abundances decreasing with altitude, and with higher levels of soil pH.In particular, OTUs assigned to the genera Mortierella, Mucor, Umbelopsis, Mollisia, and Pezicula responded negatively to increased altitude.

TA B L E 1
Variables used to parameterize the experimental treatments in the models, and their ecological interpretations.The technical implementation of the experimental treatment variables in the models is explained in Appendix S6.

Intercept
The mean value of the response variables in an average year and at an average altitude for plant individuals that did not undergo any manipulative treatment Whether the effect of exclusion is different for within-plot translocation in lower-and higher-altitude plots Low to high Categorical variable with two levels corresponding to whether the plant individual was translocated from lower to higher altitude or not Changes for a plant individual that was translocated from lower-to higher-altitudes (without exclusion), compared to a plant individual that was translocated within a plot in lower altitudes (without exclusion) High to low Categorical variable with two levels corresponding to whether the plant individual was translocated from higher to lower altitude or not Changes for a plant individual that was translocated from higher-to lower-altitudes (without exclusion), compared to a plant individual that was translocated within a plot in higher altitudes (without exclusion) Excluded × Low to high Interaction between Excluded and Low to high Whether the effect of exclusion was different for a plant individual that was translocated from lower-to higheraltitudes, compared to a plant individual that was translocated within a plot in lower altitudes

Excluded × High to low
Interaction between Excluded and High to low Whether the effect of exclusion was different for a plant individual that was translocated from higher-to loweraltitudes, compared to a plant individual that was translocated within a plot in higher altitudes

| How do plant functional traits respond to translocation across altitudes?
Bistorta vivipara leaves were generally smaller in width and length at higher altitudes.Additionally, among those plants that were translocated from lower to higher altitudes, leaf width and length decreased especially when the colonization of resident fungi was excluded (≥95% posterior probability, see Figure 3; Table S9).In other words, when plants were translocated from lower to higher altitudes, the plants that were able to associate with the local higheraltitude RAF grew more than those plants that were not allowed to associate with the local RAF.Likewise, after translocation, plant leaf TA B L E 2 Variance partitioning quantifying the proportion of variation that can be attributed to the measured fixed and random effects (rows) in each model (columns).The predictors have been grouped into those modelling natural variation (N), variation due to the experimental treatments (T), and technical variation due to the sequencing (S).S3 for a list of the fungal species.width was larger when fungal colonization was not excluded compared to when it was excluded both when plants were translocated from low to high and from high to high (Figure 3; Table S9).

| How are changes in RAF communities linked to changes in plant functional traits?
The 'fungal-plant model' revealed that specific taxonomic fungal groups influenced plant growth.The presence and abundance of certain taxonomic groups of Ascomycota (specifically, various OTUs identified to the genera Pezicula and Mollisia, the orders Capnodiales and Pleosporales, and the classes Sordariomycetes and Pezizomycetes) were associated with increased leaf length and flower presence (Figure 4).Similarly, the presence of Leucosporidium (Basidiomycota) was associated with increased leaf length and flower presence (Figure 4).Additionally, the increased abundance of some genera, including OTUs assigned to Phialocephala, Cladophialophora, and Mortierella, was associated with decreased leaf width.S9.In agreement with our main hypothesis, the results suggested that plants at higher altitudes exhibited greater growth when grown in the presence of RAF from higher altitude environments.Earlier studies based on observational data have likewise suggested that plant fitness responses to drought, heat, and salt as environmental stressors are influenced by microbial communities that were adapted to that stressful environment (Lau & Lennon, 2012;Rodriguez et al., 2008).Our results provide experimental support to these previous findings by showing that RAF improve stress tolerance to altitudinal stressors, such as colder temperatures and lower nutrient availability at higher altitudes.

| DISCUSS ION
Our results suggest that plants may not need to adapt to stress as much as we previously thought, as their RAF can alleviate that stress adaptation.The RAF may be mediating plant stress responses, which is known to affect leaf size.For instance, ectomycorrhizal fungi have a mantle, which provides a physical barrier between the plant and the environment, which can alleviate abiotic stress and negative biotic stress (Branco et al., 2022).Additionally, mycorrhizal fungi can manipulate plant signalling pathways such as the jasmonic acid and salicylic acid pathways, both of which trigger plant immunity and stress responses (Vishwanathan et al., 2020).Mycorrhizal fungi and rhizosphere-associated organisms, such as saprotrophs, can facilitate decomposition, which enables mycorrhizal fungi to transfer nitrogen and phosphorus to plants, thereby enhancing plant health and conferring stress tolerance (Hestrin et al., 2019).
Another possibility is that RAF are modulating root exudates (Frew et al., 2020), which influences recruitment of other rootassociated microbes (Emmett et al., 2021).Furthermore, RAF at high altitudes likely improve plant growth because the benefits plants derive from these associations are more important in stressful environments (Pellissier et al., 2013).For instance, nitrogen is known to enhance plant growth and fitness, but soil nitrogen tends to decrease at increased altitudes.Therefore, nitrogen acquisition from fungal symbionts is considered more beneficial at high altitudes (Pellissier et al., 2013;Wagg et al., 2011).
Alternatively, it has been previously suggested that plants may be actively filtering for beneficial fungi at high altitude where abiotic conditions are harsher (Mundra et al., 2016).In this case, plants could be selecting for certain fungal species based on functional diversity.Plants may be able to influence the nearby soil fungal communities (Dassen et al., 2017;Schmid et al., 2021) through physical or chemical defences against different types of fungi (Högberg & Högberg, 2002).However, such hypothesis cannot be disentangled from the present study.
Our results showed that plant growth traits were improved not just when specific fungal species were present, but when entire clades were present or abundant.Increased leaf length was positively associated with Pezicula and Mollisia, both of which are part of the Helotiales order, which are often found at high latitudes (Newsham et al., 2009).Pezicula species are functionally listed as a harmless endophyte, a dark septate endophyte (an endophyte often found in cold-stressed habitats; (Read & Haselwandter, 1981)), or a weak plant pathogen (Chen et al., 2016;Tedersoo et al., 2014).
Endophytes have been shown to provide habitat-specific stress tolerance to plants (Rodriguez et al., 2008).Endophytes often increase host shoot or root biomass, possibly by inducing or biosynthesizing plant hormones (Tudzynski & Sharon, 2002) or through protection against other fungal pathogens (Rodriguez et al., 2009).
Saprotrophs are a guild that receive nutrients by breaking down dead organic plant material (Rodriguez et al., 2009).Root endophytes can occur along the saprotrophic-symbiotic spectrum (Koide et al., 2008;Rodriguez et al., 2009), which may explain why the species characterized as endophytic are found both to positively and negatively affect plant growth.
The contrasting functional roles of Helotiales (Wang et al., 2006) are further validated in our study, as Helotiales species had idiosyncratic effects.Our results revealed that the Helotiales genera listed above improved plant functional growth, suggesting they may facilitate plant performance through enhanced biomass, phosphorus concentration, and nitrogen uptake (Jumpponen et al., 1998;Mullen et al., 1998;Read & Haselwandter, 1981;Schadt et al., 2001).
However, other species -particularly those in the Phialocephala genus, which are characterized as a dark septate endophyte (Newsham, 2011) -decreased plant functional growth when present in our root samples, agreeing with previous descriptions (Jumpponen & Trappe, 1998).
Several saprotrophs were found in plants with increased leaf width and length, including Trichocladium opacum, Herpotrichia pinetorum, and Phaeosphaeria.One explanation as to why some saprotrophic fungi were found to increase plant functional trait growth is that they may have been rhizosphere-associated but not necessarily endophytic in the plant samples.In this case, they could be helpful to the plant by breaking down root exudates or by mineralizing nutrients which are then absorbed by plant roots (Baldrian et al., 2011;López-Mondéjar et al., 2018).Alternatively, healthier plants may have more root exudates, so saprotrophs occur near them.
While the RAF communities showed systematic differences between higher and lower altitudes, a major part of their variation was captured by the plant individual-level random effect.
While our study design does not allow disentangling whether such unpredictable variation was generated by biotic interactions, responses to small-scale environmental variation, or some other factors, our results suggest that the variation in the RAF communities influences plant growth, pointing out specific fungal groups that facilitate plant growth.We however note that this result is of correlative nature, and hence that providing further causal evidence would require alternative experimental designs, such as inoculation experiments with specific fungi growing in higher altitudes.Furthermore, our results show that the RAF communities acquired by B. vivipara plant individuals are relatively stable, as the experimental treatments had generally only rather mild effects.
Such results could however also be because B. vivipara plants were translocated with their surrounding soil, which may have hampered the colonization by the local resident communities within a one-year period.
We conclude that interactions with RAF influence how plants respond to stressful abiotic conditions.However, it remains to be experimentally tested whether these specific species, their functional traits, or if other interactions within the plant microbiome lead to improved plant growth.Categorizing fungal species into guilds further characterizes their functional role in relation to their environment, and is often associated with their responses to warming, nitrogen, CO 2 , and plant species richness and identity (Alzarhani et al., 2019).
While knowledge of fungal functional traits is expanding (Abarenkov et al., 2010;Nilsson et al., 2019), there is still a large knowledge gap associated with arctic fungi (Abrego, Huotari, et al., 2020) Species communities are known to change based on their abiotic environment, but how these responses are related to their biotic interactions has remained largely experimentally untested.In this F I G U R E 3 Bistorta vivipara leaf width and length for each treatment before and after translocation, with statistical support (***) based on ≥95% posterior probability.Numerical values of the posterior probabilities are given in Table investigated whether biotic interactions, specifically between RAF communities and plants, affect plant growth along an abiotic stress gradient.
due largely to incomplete fungal reference databases.The potential for RAF to impact plant functional trait responses to stressors suggests that these communities could aid in preserving plant fitness amidst rapid global change.AUTH O R CO NTR I B UTI O N S Nerea Abrego conceived the idea and implemented the field experiments, Natalia Ivanova and Arusyak Abrahamyan carried out the laboratory work, Brendan Furneaux applied the bioinformatics, Otso Ovaskainen assisted in the field experiments and led the statistical analysis, and Skylar Burg contributed to the statistical analysis.Nerea Abrego and Skylar Burg wrote the first draft of the manuscript with specific contributions from Natalia Ivanova, Arusyak Abrahamyan, Panu Somervuo, Pekka Niittynen, Brendan Furneaux, and Otso Ovaskainen.All authors contributed critically to the drafts and gave final approval for publication.

Variable group Covariate Fungal model for presence-absence (%) Fungal model for abundance (%) Plant model (%) Fungal-plant model (%)
F I G U R E 2Fungal species-level responses to each predictor, both for presence-absence and abundance conditional on presence fungal models.The darker red blocks indicate a positive response of fungal species with ≥95% posterior probability, and the darker blue blocks indicate a negative response with ≥95% posterior probability.The lighter blue and red blocks indicate a response with ≥75% posterior probability.See Table