Can root-associated fungi mediate the impact of abiotic conditions on the growth of a High Arctic herb?

Arctic plants are affected by many stressors. Root-associated fungi are thought to influence plant performance in stressful environmental conditions. However, the relationships are not transparent; do the number of fungal partners, their ecological functions and community composition mediate the impact of environmental conditions and/or influence host plant performance? To address these questions, we used a common arctic plant as a model system: Bistorta vivipara. Whole plants (including root system) were collected from nine locations in Spitsbergen (n=214). Morphometric features were measured as a proxy for performance and combined with metabarcoding datasets of their root-associated fungi (amplicon sequence variants, ASVs), edaphic and meteorological variables. Seven biological hypotheses regarding fungal influence on plant measures were tested using structural equation modelling. The best-fitting model revealed that local temperature affected plants both directly (negatively aboveground and positively below-ground) and indirectly - mediated by fungal richness and the ratio of symbio- and saprotrophic ASVs. Fungal community composition did not impact plant measurements and plant reproductive investment did not depend on any fungal parameters. The lack of impact of fungal community composition on plant performance suggests that the functional importance of fungi is more important than their identity. The influence of temperature on host plants is therefore complex and should be examined further.


Introduction
Arctic plants are facing many environmental constraints for growth, such as short vegetation season, consistent cold, limitation of nutrients or cyclic physical disturbances, i.e. cryoturbation (Billings, 1987). These plants have evolved a range of adaptations to cope with the prevailing conditions, including being perennial and allocating most of their biomass below-ground (Iversen et al., 2015;Poorter et al., 2012;Qi et al., 2019). Being perennial provides a resource-saving advantage in nutrient-poor habitats with low temperatures that slow down biochemical reactions and therefore also growth; whereas the benefits of biomass allocation to below ground parts include increased area of nutrient absorption. The interface between plant and soil is thought to mediate plants' stress tolerance (Hou et al., 2021). Because of common nutrient scarcity and a large portion of plant biomass located underground this interface seems to be of relatively greater importance in the Arctic than in other biomes (Poorter et al., 2012). A significant part of the soil-plant interface is inhabited by microbes, including root-associated fungi (RAF). Arctic plant RAF consist mostly of symbiotrophic fungi, especially ectomycorrhizal (Bjorbaekmo et al., 2010;Blaalid et al., 2012;Mundra et al., 2015 a). These fungi efficiently increase the volume of soil that can be penetrated in search for resources, such as nutrients from seasonally or newly thawed permafrost (Hewitt et al., 2020). The most severe limitations for growth observed in arctic plants are low temperatures and resource limitation (Billings, 1987;Chapin and Shaver, 1985), suggesting that the relationship with RAF might play a crucial role in plant survival and growth.
Multiple characteristics of species communities play an essential role in the functioning of ecosystems, such as richness, abundance or community structure (Hooper et al., 2005;Loreau, 2001;Maestre et al., 2012;Tilman et al., 2012). Based on previous findings, we expect that the more diverse the RAF community, the better for a host plant (Baxter and Dighton, 2001). However, it is not clear how these characteristics of RAF communities impact their host plants, especially in cold biomes. Symbiotic fungi provide resources and probably additional benefits such as mitigating harmful effects of environmental stressors, thus enhancing plant growth and productivity (Wardle et al., 2004). However, releasing root exudates of primary metabolites that can be absorbed by members of its microbiota does come with a cost for a plant (Brzostek et al., 2014;Shi et al., 2015). In nitrogen-limited tundra in Alaska, 61-88% of plant nitrogen was supplied from mycorrhizal fungi; in exchange, the plant delivered 8-17% of carbon produced photosynthetically to the fungi (Hobbie and Hobbie, 2006). A plant could perhaps increase the amount of released nutritious root exudates to attract more species of symbiotrophic fungi that in turn could increase the amount of nitrogen delivered. However, higher fungal richness would increase competition for limited space in the rhizosphere and possibly for resources, although the mechanism is not yet fully understood. Therefore, plants 'living on the edge' in the High Arctic may benefit from selectively choosing their RAF community members, favouring the most beneficial fungal partners for plant growth or stress mediation (Mundra et al., 2016). In this scenario, species richness in RAF communities would be irrelevant for plant performance. On the other hand, the presence of specific functional traits rather than the identity of particular species could be more important (Louca et al., 2017). The vast array of interconnected biotic and abiotic factors occurring in natural systems complicate uncovering if and how plants prefer specific root-associated fungi from the pool of species present in the soil (Jones et al., 2019).
One approach to disentangle these often-confounded factors are controlled experiments. Most of the experiments assessing the impact of RAF diversity on host plant performance have focused on arbuscular mycorrhiza in crops (Begum et al., 2019); whereas similar studies on ectomycorrhizal (EcM) plant species come mostly from the pre-high throughput sequencing era and have focussed on trees (e.g. Baxter and Dighton, 2001). Several experiments under controlled settings have shown that EcM host plants benefit from increased fungal richness, however, the tested level of richness was often incomparable with natural environments, such as an increase from 1 to 4 species of EcM fungi (Baxter and Dighton, 2001). Some studies, however, did not find enhanced plant performance mediated by EcM fungi or concluded that the influence of EcM species richness on plant productivity is context dependent (Jonsson et al., 2001). RAF diversity was shown to be particularly sensitive to experimental conditions compared to fungi that inhabit space further from the roots in the rhizosphere or bulk soil (Almario et al., 2017). Moreover, morphology and physiology of lab-grown plants differ from those in the natural system, e.g. by increased growth rate and higher concentrations of nutrients in tissues (Poorter et al., 2016). All these differences could affect and alter plant-associated organisms, such as RAF. Experimental procedures cannot consider all the complexity of natural systems and their effects do not always reflect those observed in the wild. Thus, observational studies can provide crucial complementary knowledge, in particular for extreme environments like the High Arctic.
Species response to environmental shifts, including ongoing climate changes, is one of the crucial questions in natural sciences. It is a particularly outstanding issue in the Arctic where temperatures are increasing at the fastest pace in the world, and are predicted to continue rising rapidly (Bintanja and Selten, 2014;Post et al., 2019). These changes impact mechanisms that alter biogeochemical cycles and determine critical ecosystem-climate feedback processes, such as the release of organic carbon pool of which nearly half of the global stock is stored in the Arctic soils (Schuur et al., 2015;Tarnocai et al., 2009). Such ecosystem feedbacks, which are essential bricks in the understanding of global change, depend on complex relationships between abiotic and biotic factors in arctic soils (Wookey et al., 2009). However, the biology of these arctic soils remains at present an understudied 'black box' (Metcalfe et al., 2018;Virkkala et al., 2019).
To shed some light onto these soil processes, we used a plant-centric approach to study the impact of the root-associated fungal community on the growth and reproductive investment of a wide-spread arctic plant, Bistorta vivipara. We took into account the most important abiotic factors, both edaphic and climatic, affecting the host plant and its RAF community. We used structural equation modelling (SEM) to assess whether the fungal community mediates the effect of abiotic conditions on plant performance and to disentangle direct from indirect effects. We tested the following hypotheses: (i) Plant morphological measurements (considered as a proxy for plant performance) depend both on abiotic conditions and on the fungal community, and (ii) only richness and functional traits, but not the specific species composition of the RAF community affects plant morphology. Moreover, we tested which measurements of plant parts involved in different processes such as energy storage, energy acquisition and reproduction depend among others on the RAF community.

Datasets
We combined and reanalysed datasets from nine different locations collected in 2012 and 2013 in Spitsbergen, the largest island of the High Arctic Archipelago Svalbard, Norway (Table 1; Fig. 1; precise description of location and environmental conditions of each sample is available at https://github.com/magdawutkowska/bistorta/blob/master/en v_bistorta_temp.txt). Each dataset consisted of host morphology and molecular description of the RAF community, together with associated Table 1 Overview of the data included in this study. Each dataset was generated to investigate specific topics regarding Bistorta vivipara root-associated fungi (RAF). References are given for previously published data.  (Table 1). Each of the studies established a randomized sampling scheme in the locality of choice, also assuring that sampled plants were of different age. B. vivipara grows, sets leaves and develops inflorescence early, which means that even plants sampled in June were already fully grown. Whole plants with an intact root system were excavated. To explore the associations between plant performance, allocation patterns and its environment we measured three morphological features of the B. vivipara individuals hosting the analysed RAF communities (Supplementary 1). The rhizome is an underground storage organ that accumulates assimilated biomass as non-structural carbohydrates, and was therefore used as a proxy for overall plant performance (Hartmann and Trumbore, 2016). Rhizome dimensions were measured and used to calculate an approximate volume (RV) by multiplying its length, height and width. Length of the longest stem leaf (LL) was used as a proxy for photosynthetic capability of the plantthe longer the leaf, the bigger the photosynthetic area. In the upper part of the stem, B. vivipara produces flowers and bulbils for sexual and asexual reproduction, respectively. We used the ratio of the length of the stem covered by flowers and bulbils (inflorescence) to the total stem length (I/S), as a proxy for the plant's investment in reproduction.

Meteorological and edaphic variables
Meteorological data were obtained for each sampling point from the high-resolution 1 km-gridded dataset Sval_Imp_v1 (Schuler and Østby, 2020). We extracted the sum of average monthly precipitation (p) and average July air temperature (t), both from the year of sampling. Svalbard is characterized by steep temperature gradients and the main climatic contrasts in our study were between sites. However, to take into account the two different years of sampling at the Isdammen site, we used weather data of the sampling year, as interannual variability in temperature may impact the RAF community. Soil samples were collected from the same sampling spot as plants. The following edaphic parameters, representing critical properties of the abiotic environment, were measured in all datasets: pH, soil nitrogen concentration (N) and carbon to nitrogen ratio (C/N; used as an indicator for soil nitrogen availability or soil fertility). Edaphic variables were obtained in the same way for all datasets (described in detail in Mundra

Fungal data
Bistorta vivipara roots were cleaned within a day after sampling and fixed in a 2% CTAB extraction buffer until DNA extraction (details described in each of the publications; Table 1). All datasets targeted the same fragment of internal transcribed spacer 2 amplified with fITS7a forward primer (Ihrmark et al., 2012) and reverse primer ITS4 (White et al., 1990). The amplified fragments were sequenced with Illumina MiSeq (300bp paired-end reads).
Each dataset was a mixture of sequences located in 'forward' and 'reverse' direction. Thus, first, a mapping file with variable length barcodes and primer sequences was used to identify sequences in each location using sabre (https://github.com/najoshi/sabre) and generating separate R1 and R2 files for each read direction. Next, primers were clipped, and sequences with ambiguous bases (Ns) were removed using cutadapt v. 2.5 (Martin, 2011). Python script FastqCombinePairedEnd. py (https://github.com/enormandeau/Scripts) was used to assure that each sequence had its pair and were in the matching order for further analyses. We used an amplicon sequence variants (ASVs) approach implemented in DADA2 v. 1.11.1 (Callahan et al., 2016) and executed in R v. 3.5.2 (R Core Team, 2018; for details see Supplement 2 and scripts generated for this study). The datasets were analysed using DADA2 ITS workflow (https://benjjneb.github.io/dada2/ITS_workflow.html). Fungal data were produced independently for each study; therefore, they were initially analysed separately due to different error rates for each sequencing run. Separate ASVs tables were then merged. A consensus method was used to remove chimeras (3759 out of 11243 input sequences). Sequences shorter than 200bp and six samples with a very low number of reads were removed. Due to profound differences in depth of sequencing the ASV table was randomly subsampled (21639 reads per sample; the number of detected ASVs before and after subsampling was highly correlated; Kendall's τ = 0.95). Taxonomy was assigned using the RDP naive Bayesian classifier implemented in DADA2 using the full UNITE + INSD reference dataset for fungi (UNITE Community, 2019; sh_general_release_dynamic_02.02.2019). All the ASVs were functionally annotated using the FUNGuild database (Nguyen Fig. 1. Bistorta vivipara plants from the four concatenated datasets were collected from nine localities on Spitsbergen. The localities represent different habitats: glacier forefronts (Renardbreen and Hørbyebreen), soil around hot springs (Trollkjeldene), arctic steppe (Ringhorndalen), hydrocarbonrich site (Kvalvågen), nutrient-rich mine-contaminated site (Bjørndalen, tailings from a mine), nutrient-rich tundra located under bird cliffs (Vestpynten) and two natural tundra patches both located in Adventdalen (snow fence experimental setup and Isdammen). et al., 2016).
Differences in community composition were summarized through non-metric multidimensional scaling (GNMDS; vegan package (Oksanen et al., 2019). The first axis that most clearly separated localities and habitats was used as a proxy for composition in further analyses. We used both presence-absence based metrics and parameters based on read abundance to describe RAF communities: ASV diversity (D), a ratio of symbio-to saprotrophs (Sy/Sa) and GNMDS values for 1st axis as a proxy for community composition (CC; Table 2).

Statistical analyses and model selection
The statistical analyses were performed in R v. 3.5.2 (R Core Team, 2018). Based on available literature of soil and weather influence on fungi and plant interactions in the Arctic (Table 3), we built seven hypothetical causal path models relating abiotic variables to the three metrics characterizing the fungal community and plant morphological measurements (solid lines in Fig. 2). The unbranched rhizome of B. vivipara elongates with age, providing space for new roots to stem from its distal end (Diggle, 1997) and therefore increasing the richness of recruited RAF . Randomised sampling schemes in each of the studies included in our study excluded the potential influence of plant age on the results. For the full model, we assumed that all three fungal parameters influence all three plant measures, in addition to abiotic factors impacting both fungal and plant variables.
There is a growing body of evidence suggesting that RAF might not always provide benefits for its hosts (Ågren et al., 2019;Corrêa et al., 2006;Noë and Kiers, 2018). Additionally, resources provided by RAF can be used to address many current needs within a host, which perhaps could not be manifested in a larger investment in growing different body parts. Thus, we tested models hypothesizing that changes in fungal richness and community composition along the sample gradient are not essential for specific plant measurements. Therefore, in the three subsequent models, we preserved all the relationships omitting only the fungal variables' link to a specific plant response (I/S, RV or LL was not influenced by fungi). In the next model, we, hypothesized that CC is not an important parameter for any of the plant measurements. Additionally, we combined this last model with the best model obtained from simplifying the relationships between fungi and plants responses.
Finally, in our null model, to evaluate whether fungal parameters have any impact on plant measurements, we removed all connections between fungal parameters and plant measurements. In all the models, we treated edaphic and meteorological variables as independent. We are aware that they can affect each other, but this was not the focus of the study. The strongest correlation among them was between N and C/N (r = − 0.64). We also did not hypothesize any causal links between the fungal parameters. Concerning the plant variables, we assumed a causal link between rhizome volume and leaf length, because leaf growth in the start of the season depends on stored resources. Locality was used as a random effect in all the models because: (i) fungal community composition usually shows a high spatial variation (e.g. Blaalid et al., 2014), (ii) preliminary ordinations showed that in our dataset fungal communities differed between localities and (iii) plants may have local adaptations influencing growth. Additionally, for the site we accounted for in-between year variation by using year-specific climatic variables.
We applied structural equation modelling (SEM) to carry out an exploratory path analysis of these models, using the psem function in the piecewiseSEM package (Lefcheck, 2016). The SEM was composed of linear mixed-effects models (LMMs) for each fungi parameter and plant measurement, which was fitted using the lme function in nlme package (Pinheiro et al., 2020). The fit of the separate LMMs were assessed graphically for normality of the residuals. Residuals clearly deviating from the expected distribution on a quantile-quantile plot with standardised residuals > |3| were considered as outliers and these data points were therefore excluded.
The analysis was performed using both presence-absence and read abundance metrics for the fungal community. Because some of the fungal parameters were correlated, we included non-directed correlations among them in the SEM to make it possible to estimate the paths in Table 2 Metrics used to describe the fungal community used in this study for presence absence data and number of reads, respectively. All the parameters were calculated using a rarefied table containing amplicon sequence variants (ASVs).

Fungal parameter
Presence-absence   our exploratory model. This was the case for CC and Sy/Sa based on presence-absence and for Sy/Sa and D based on read abundance. The distributions of all variables were assessed graphically, and some were log-or logit-transformed to assure roughly normal distributions. All variables were scaled to 0 mean and a standard deviation of 1 to make effect sizes comparable. A prerequisite for a SEM model to be considered as fitting was Fisher's C p-value > 0.05 (Shipley, 2013). The best models among the candidate sets described above were chosen based on the lowest AIC values. Both of these values were calculated within the psem function. We used statistically significant estimates from the best fitting presence-absence model to calculate indirect effects of abiotic factors on plant measures.
The combined dataset consisted of 214 B. vivipara plant measurements with associated edaphic data and corresponding RAF data. For the SEM, we excluded all observations with missing values resulting in a final dataset with 188 plants (after excluding outliers the presenceabsence dataset had 187 data points and the abundance dataset had 185).

Models based on presence-absence fungal parameters
The best-fitting presence-absence path model (Table 4) supported the hypothesis that fungal CC does not impact plant measurements within the range of CC variation in the sampling gradient in our study. Simultaneously no fungal parameters affect the I/S. The second bestfitting model with a relative difference ΔAIC <1, supported the related hypothesis that I/S does not depend on any of the fungal parameters included in this study, but included the effect of CC on other plant parameters.
In the best-fitting and most parsimonious model, fungal community richness and the ratio of symbiotrophic to saprotrophic species were related to plant measurements as follows (Fig. 3): fungal richness was positively related to RV (path coefficient ± standard error = 0.26 ± 0.07, p < 0.001); full list of all the effect sizes in Supplement 4a) and Sy/ Sa was negatively related to LL (− 0.20 ± 0.07, p = 0.004). Except for the fungal metrics, the RV also showed positive correlations with p (0.29 ± 0.11, p = 0.01). LL was negatively impacted by N content (PC ± SE = − 0.20 ± 0.08, p = 0.02) and t (− 0.34 ± 0.08, p < 0.001). The highest estimate in our model suggested a positive correlation between RV and LL (0.53 ± 0.06, p < 0.001).
Meteorological data had a clear effect on fungal parameters: p had a positive effect on Sy/Sa (0.44 ± 0.21, p < 0.04), and t on fungal CC (0.27 ± 0.09, p = 0.003), but t had a negative effect on D (− 0.45 ± 0.13, p < 0.001). Based on the best fitting presence-absence model, edaphic variables did not seem to impact any fungal parameters and plant measurements except the already mentioned N content impact on LL. On the other hand, t correlated with multiple fungal and plant variables.
Among abiotic factors impacting plant measurements, t affected LL over three pathways: one direct (negative, path coefficient = − 0.34) and two indirect: positive through RV (path coefficient = 0.29 * 0.53 = 0.154) and negative through fungal D (path coefficient = − 0.45 * 0.26 = − 0.117). The direct effect was therefore the strongest and the two indirect effects were of comparable magnitude, but opposite directions.

Abundance model
The best-fitting path model based on read abundance supported the hypothesis that fungal parameters do not impact any plant measurements within the range of fungal parameters' variation in the sampling gradient in our study (Table 4). Another model that differed by ΔAIC = 0.25 supported the same hypothesis as the best fitting presence-absence model: fungal CC does not impact plant measurements, and I/S is not affected by other fungal parameters.
Although the role of fungi in the best model based on abundance differed substantially from the best model based on the presenceabsence ASV table, some of the statistically significant relationships between environmental variables and plant measurements were preserved (Fig. 3, full list of all effect sizes from both types of models in Supplement 4). This included the negative correlations between N content and LL (− 0.23 ± 0.08, p = 0.005), t and LL (− 0.37 ± 0.08, p < 0.001), as well as the effect of t on CC (0.31 ± 0.09, p < 0.001). Also, the relationship between the two plant variables RV and LL, showed the same magnitude as in the best fitting presence-absence model (0.54 ± 0.06, p < 0.001). The best fitting abundance model did not support any indirect effects of abiotic factors mediated by fungal parameters.
The abundance-based model revealed links between edaphic and fungal parameters that were not statistically significant in the presence- Fig. 2. Schematic illustration of a conceptual plant-centric model representing relationships between variables suggested by the literature and tested in this study. Solid lines are associations that were researched by studies from the Arctic; dashed lines were described by fewer studies, mainly from other regions. The full model includes all possible links between each abiotic, fungal and plant variable. Abbreviations and symbols: N -soil nitrogen content; C/N -the ratio of soil nitrogen to soil carbon content; p-precipitation; t -temperature; D -diversity; Sy/Sa -the ratio of symbio-to saprotrophs; CC -fungal community composition; I/S -the ratio of inflorescence to stem length; RV -rhizome volume; LL -leaf length of the longest leaf.

Table 4
Summary of the models and statistics used to select the best fitting model. Each model reflects a separate hypothesis. The full model includes all possible links between each fungal variable and each plant variable. Subsequent models exclude some of the links, as indicated in the name of each model. Abbreviations: I/S -ratio of inflorescence to stem length; RV -rhizome volume; LL -leaf length; CC -root-associated fungal community composition. The best model for each approach is highlighted in bold. Models which don't fit based on the test of directed separation are in italics.

Variance in fungal and plant response variables
In both best fitting models, the variance in plant measurements was on average better explained by fixed factors than the variance in fungal parameters (marginal R 2 = 0.02-0.44 vs 0.07-0.26, Table 5). However, overall, the variance explained by fixed factors was rather low. On the other hand, locality included as a random factor explained on average more variation in fungi than in plants (conditional R 2 -marginal R 2 = 0.03-0.58 and 0.01-0.33, respectively). The high proportion of variance explained for fungal response variables was especially pronounced in presence-absence compared to the abundance model (conditional R 2marginal R 2 = 0.40-0.58 and 0.03-0.48, respectively).

Discussion
Our exploratory study revealed that measurements of below-and aboveground plant organs responded in opposite ways to average July temperature, the effects of which were both direct and mediated by parameters of the RAF community. Regarding fungal parameters, both species richness and functional diversity were important for plant performance measurements, but not the specific community composition.
Our study revealed that among the abiotic factors temperature was the most important element, reflecting its immense significance in physical constraints for arctic biota (Billings, 1987) and the general tendency of modifying interactions between organisms (Bideault et al., 2019). However, our results also suggest that the impact of temperature on an arctic plant is far more complex than previously thought (Chapin, 1983;De Long et al., 2015;Rinnan et al., 2020) and perhaps unpredictable (De Long et al., 2019). Our results showed different influences of temperature on below-and aboveground plant measurements (Fig. 3), which question current methods of monitoring changes in arctic vegetation, such as the use of normalized difference vegetation index (NDVI) as a proxy for plant biomass. This technology advanced the understanding of vegetation biomass dynamics simultaneously over vast and otherwise under-sampled areas of the Arctic (e.g. Myers-Smith et al., 2020;Myneni et al., 1997;Phoenix and Bjerke, 2016). However, it is based on remote measurements of Earth's surface reflectance, and therefore takes into consideration only aboveground changes in foliage. In these methods plants' below-ground productivity and biomass are disregarded, probably resulting in underestimation of the overall impact of increased temperatures on plants, such as B. vivipara, which is an ubiquitous species in the Arctic and essential food source for ptarmigans (Steen and Unander, 1985), geese (Anderson et al., 2012) and reindeer (Bjørkvoll et al., 2009). Temperature had a direct opposite effect of similar magnitude on LL and RV (− 0.34 vs 0.29, respectively), additionally complicated by indirect fungal effects. This suggests that NDVI can easily underestimate the impact of warming on overall plant biomass and confuse understanding of carbon stocks dynamics. Presently, there are no tools that can be used to scan below-ground plant biomass at scales similar to NDVI. However, more laborious in situ methods, e.g. minirhizotrons, are used to measure below-ground biomass (Wilson, 2014). Their use significantly enhances our understanding of the dynamics in belowground biomass allocation. Nevertheless, temperature affecting a host plant through multiple direct and indirect pathways generates major difficulties in projections of the future response of ecosystems to warming.
The temperature also affected two measures of RAF: diversity and community composition. Parameters associated with RAF communities, such as richness and the ratio between symbiotrophic and saprotrophic ASVs, impacted the plant both positively and negatively, thus balancing themselves out. The mechanism behind fungal mediation of temperature Fig. 3. Path diagram showing resulting connections between predictor and response variables in the best fitting models. Statistically significant (p < 0.05) links are depicted by arrow colours (positive or negative nature of the relationship) and thickness (relationship magnitude); the numbers are estimates from the models. Abbreviations and symbols: N -soil nitrogen content; C/Nthe ratio of soil nitrogen to soil carbon content; p -precipitation; t -temperature; D -diversity; Sy/Sa -the ratio of symbio-to saprotrophs; CC -fungal community composition; I/S -the ratio of inflorescence to stem length; RVrhizome volume; LL -leaf length of the longest leaf. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Table 5 Proportion of variance explained without (marginal R 2 ) and with random factors (conditional R 2 ). Locality was used as a random factor in all of the models. Abbreviations: D -diversity; Sy/Sa -the ratio of symbio-to saprotrophs; CCfungal community composition; I/S -the ratio of inflorescence to stem length; RV -rhizome volume; LL -leaf length of the longest stem leaf. is not clear. However, there are other molecular and physiological characteristics that could explain the influence of fungi on plant performance mechanistically. For instance, secretion of fungal signalling molecules, such as volatile organic compounds (Schenkel et al., 2018) or plant-like hormones (Chanclud and Morel, 2016;Han and Kahmann, 2019) that can be translocated to host plant cells and there elicit a physiological response. Release of these molecules could be temperature-dependent. Similarly, plant-based responses to these signals could also be at least partly temperature-dependent, e.g. release of root exudates (Canarini et al., 2019). The finding of negative impact of nitrogen on leaf length was unexpected in light of previous findings ( Fig. 3; Wookey et al., 1994). Bistorta vivipara is regarded as a pioneer plant (Dormann et al., 2002), able to cope with severe conditions and resource limitations (Davey et al., 2015;Mühlmann et al., 2008). In a High Arctic nitrogen-rich habitat, such as bird cliffs, where the competition between organisms is high, it is most likely outcompeted by other plants. Additionally, these nutrient-rich habitats are characterised by an increased number of plant interactions with herbivores, such as reindeers or geese, that can eliminate foliage and may change the ratio of below to above ground plant parts.
Almost all symbiotrophic RAF of B. vivipara in Svalbard are ectomycorrhizal (Davey et al., 2015;Mundra et al., 2015 a). Since these fungi exchange nitrogen with plants in return for versatile carbon metabolites (Hobbie and Hobbie, 2006), we hypothesized that in a resource-limiting environment this fungal trophic mode could promote bigger plants (Tedersoo et al., 2020) with bigger leaves. This way, fungi could potentially influence the number and amount of metabolites that the plant could produce in return and share in its rhizosphere. However, our results showed the opposite scenario, where Sy/Sa had a negative effect on leaf length, suggesting that more fungal partners enhance competition over scarce resources (Kennedy, 2010). When accounted for separately, the richness of symbio-and saprotrophs did not show any associations with plant measurements (data not shown). However, the ratio of their richness did, perhaps reflecting the characteristics of soil conditions in different localities. Small ratio of Sy/Sa was found in localities with little organic matter (Supplementary 3), suggesting that this parameter mirrors fertility properties of soil (soil organic matter content could not be included in the SEM analysis, because it was lacking for 5% percent of the samples). When soil organic matter content is low, colonizing plant roots ensure fungal access to an easily accessible pool of carbon from root exudates (Baldrian and Kohout, 2017). Although the B. vivipara root system is relatively compact and flexible, growing in mineral soils could promote longer roots to assure access to quickly drained soil water, e.g. at early stages of soil development in glacier forefronts (Frenot et al., 1998). Intense disturbance caused by periglacial processes in these habitats may contribute to physical breaks in fine roots or associated fungal mycelium, perhaps leading to an increase in the number of saprotrophic fungal species. Alternatively, saprotrophic fungi could be among the first organisms in primary community assembly, based on their ability to use organic carbon from heterotrophic communities of invertebrates that feed on allochthonous organic matter. This is now recognized as a crucial step in primary succession before establishment of autotrophs (Hodkinson et al., 2001(Hodkinson et al., , 2002Jumpponen, 2003).
Our finding that fungal community composition does not affect plant measurements could perhaps originate from strong environmental filtering of root-associated fungal communities (Blaalid et al., 2014). High physicochemical heterogeneity of arctic soils corresponds with distinct RAF community composition observed at different scales (Bjorbaekmo et al., 2010;Blaalid et al., 2012;Botnen et al., 2019;Mundra et al., 2015 a). A set of physicochemical conditions that translate into ecological niches selects fungal species that can withstand and thrive in these locality-specific combinations of factors. Principally, abiotic factors have previously been shown to affect fungal parameters (Botnen et al., 2019;Canini et al., 2019;Ni et al., 2018;Siciliano et al., 2014). Relationships between variables collected from literature searches (Table 3) were, in general, poorly reflected in the results of our models (Fig. 3). In most cases, we saw no effect of these abiotic drivers on neither plants nor fungi. It was especially pronounced in RAF community composition, suggesting other sources of the differences that are specifically connected to locality (Mundra et al., 2016). These could be other edaphic factors not included in this study (e.g. phosphorus (Darcy et al., 2018) or heavy metal concentrations (Hanaka et al., 2019), competition (Bell et al., 2013;Kennedy, 2010) or other factors that historically impacted the community assembly (Nemergut et al., 2014)). Nevertheless, the fact that arctic ectomycorrhizal RAF display little or no affinity to host species (Botnen et al., 2014) suggests that the fungal contribution to plants reflects mitigation of locality-specific conditions, rather than individual species needs. Similar conclusions were made in edge soil habitats beyond the Arctic. For instance, RAF communities in soil characterised by combined effects of poor nutritional and water status (Marasco et al., 2018) or high contamination levels (Gil-Martínez et al., 2018) seem to also be host-independent and highly variable between sites. In this study we did not look into the percentage of root colonisation in the different habitats, taxonomically specific analyses, the fraction of fungi that are metabolically active or their specific activities, but we think that these could be next steps in future analyses.
To explain discrepancies in results between presence-absence and read abundance models, it is necessary to identify possible sources of variation in read abundances in fungal metabarcoding studies. Fungal species vary in the copy number of ribosomal DNA (rDNA; 14-1442), and this number is independent of genome size or ecological roles, such as guild or trophic mode (Lofgren et al., 2019). Strains of the same fungal species, especially yeast, can exhibit high variation of rDNA copy number (Kwan et al., 2016;Liti et al., 2009). Relative abundances of reads are sometimes used as a proxy for the relative biomass contributions of some species (Deagle et al., 2019). However, a quantitative meta-analysis found only a weak relationship between the two (Lamb et al., 2019). Read abundance can be profoundly affected by methodological biases at several steps during metabarcoding procedures, starting from the choice of primers through wet-lab methods, including sequencing, to bioinformatic pipelines (Lindahl et al., 2013;Nguyen et al., 2015;Song et al., 2015;Taylor et al., 2016). However, in our study, the main pathways affecting plants directly and not through fungal parameters remained present in both best-fitting models. This supports prevalence of a biological signal over methodological biases from abundance data. On the other hand, the abundance-based model in this study showed clear links between fungal parameters and soil fertility (N and C/N) mirroring the stoichiometric state of the environment (Elser et al., 2000) and temperature that controls the rate of biochemical reactions.
Here we demonstrated that fungal parameters, such as richness and functional diversity, could mediate the influence of abiotic factors on host plants, but the underlying mechanisms remain unknown. It is not clear how different fungal species contribute to plants' biometrics, how many resources are being exchanged with plants and how that changes with RAF variation in time and space. Not only molecular identification, but also establishing biomass estimations for both fungi and bacteria could help to understand below-ground dynamics. Low proportion of variance explained by fixed factors showed that there is a strong need to obtain and include more abiotic and biotic variables that were not considered in this study but are of high importance for fungi and plants. Controlled experiments could potentially help to address these uncertainties. Additionally, morphological characterization of multiple plant species, biomass and nutrient concentration measurements in separate plant parts would ensure precise comparisons between plant life strategies in variable habitats and distant locations. Another critical aspect in making these links is to include the host plant genotype to accurately tie its phenotype to the influence of the environment (de Villemereuil et al., 2016). A comprehensive interdisciplinary study employing various methods could help to develop a mechanistic understanding of links between above-and below-ground biota, including other taxonomic groups. Elucidating these functional relationships between the components of Arctic ecosystems is a critical step towards understanding the dynamics of soil carbon pools and decrease associated uncertainties (Wieder et al., 2019;Zha and Zhuang, 2018).

Author contributions statement
MW, DE, SM, PB, AV -wrote and edited the manuscript; MW, AVanalysed sequencing data; MW, DE -did statistical modelling; SMcollected and processed soil samples, measured abiotic parameters, generated sequencing data, PB -developed the idea and provided funding.

Additional information
The already published datasets are available online at Dryad repository (Botnen et al., 2020;Mundra et al. 2015Mundra et al. , 2016. The other dataset containing B.vivipara root-associated fungal sequences used in this study are available at 10.5281/zenodo.4301997. Environmental data, scripts generated for bioinformatic and statistical analysis are available at https://github.com/magdawutkowska/bistorta.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.