Accounting for environmental variation in co‐occurrence modelling reveals the importance of positive interactions in root‐associated fungal communities

Understanding the role of interspecific interactions in shaping ecological communities is one of the central goals in community ecology. In fungal communities, measuring interspecific interactions directly is challenging because these communities are composed of large numbers of species, many of which are unculturable. An indirect way of assessing the role of interspecific interactions in determining community structure is to identify the species co‐occurrences that are not constrained by environmental conditions. In this study, we investigated co‐occurrences among root‐associated fungi, asking whether fungi co‐occur more or less strongly than expected based on the environmental conditions and the host plant species examined. We generated molecular data on root‐associated fungi of five plant species evenly sampled along an elevational gradient at a high arctic site. We analysed the data using a joint species distribution modelling approach that allowed us to identify those co‐occurrences that could be explained by the environmental conditions and the host plant species, as well as those co‐occurrences that remained unexplained and thus more probably reflect interactive associations. Our results indicate that not only negative but also positive interactions play an important role in shaping microbial communities in arctic plant roots. In particular, we found that mycorrhizal fungi are especially prone to positively co‐occur with other fungal species. Our results bring new understanding to the structure of arctic interaction networks by suggesting that interactions among root‐associated fungi are predominantly positive.


| INTRODUC TI ON
Understanding the role of interspecific interactions in shaping ecological communities is one of the central goals in community ecology (Bairey, Kelsic, & Kishony, 2016;Werner & Peacor, 2003;Wootton, 1994). Species can interact in two ways: directly, when individuals of one species affect the fitness of a second species, or indirectly, when individuals of one species affect the fitness of a second species through their direct interaction with a third species. In interaction networks, such as those formed by plants and microbes in the rhizosphere, direct and indirect interactions play important roles in determining the community composition and fitness of both plants and microbes (Artursson, Finlay, & Jansson, 2006;Frey-Klett, Garbaye, & Tarkka, 2007;Gould et al., 2018;Johansson, Paul, & Finlay, 2004). Furthermore, interactions in the rhizosphere influence ecosystem processes, such as nutrient cycling (e.g., Lindahl, de Boer, & Finlay, 2010;Toro, Azcon, & Barea, 1997). Therefore, unravelling the structure of rhizosphere interactions and in particular disentangling direct and indirect interactions has important implications for understanding both plant diversity and ecosystem functioning.
In this paper, we focus on root-associated fungi, which form highly interactive communities with important effects on host plant fitness (Rodriguez, White, Arnold, & Redman, 2009;Smith & Read, 2008).
Root-associated fungi directly and indirectly interact with each other, at the same time as they interact with the host plants. Mycorrhizal fungi facilitate nutrient and water uptake as well as protect roots from pathogens and toxic compounds. In return, they gain direct access to carbohydrates from plant roots (Smith & Read, 2008). A less studied group of root-associated fungi are endophytic fungi, which use living plant tissues as their habitat and may sometimes improve plant fitness (Newsham, 2011;Rodriguez et al., 2009). Yet, root-associated fungal-plant interactions are not always positive, and mycorrhizal and endophytic fungi may also act as parasites (Johnson, Graham, & Smith, 1997;Rodriguez et al., 2009).
Within host plants, root-associated fungi may interact directly or indirectly with each other in a number of ways. Direct interactions among root-associated fungi involve trophic relationships in which some species constitute the resource for others (Lindahl et al., 2010;Tedersoo et al., 2009). Another type of interaction is direct competition during colonization of roots and for root resources (Moeller & Peay, 2016). Root-associated fungi may also interact indirectly, by some species modifying exudate production by plant roots (Ingham & Molina, 1991). However, the outcomes of interactions between fungal species are not always consistent and may vary with the timing of root colonization (i.e., priority effects, see Kennedy & Bruns, 2005;Kennedy, Peay, & Bruns, 2009) or with the abiotic environment (Erland & Finlay, 1992;Mahmood, Finlay, Fransson, & Wallander, 2003).
There is increasing evidence that interactions among co-occurring microbes within host plant individuals affect plant fitness.
Fungi and bacteria can modify the nutritional influence of each other on the associated plants (Frey-Klett et al., 2011;Johansson et al., 2004). The so-called "Mycorrhiza Helper Bacteria" stimulate the symbiotic association between plants and mycorrhizal fungi, ultimately improving plant growth (Founoune et al., 2002;Frey-Klett et al., 2007). Among root-associated fungi, mycorrhizal fungi and endophytes are usually studied separately, and thus their interactions and the effects of their interaction outcomes on plant fitness are largely unknown. Some authors have found positive co-occurrence patterns between endophytic and mycorrhizal fungi and proposed that the positive co-occurrences may mostly be due to them occupying different niches within plant roots (Wagg, Pautler, Massicotte, & Peterson, 2008), or alternatively due to their facilitative interactions (Yamamoto et al., 2014). In contrast, Wearn, Sutton, Morley, and Gange (2012) found that mycorrhizal and endophytic fungi mostly avoided each other within herbaceous plants, and suggested that mycorrhizal fungi and endophytic fungi interact antagonistically by competing for plant roots.
Evidence of interactions among root-associated fungi is best achieved through direct observations in the laboratory. However, root-associated fungal communities are highly species-rich and laboratory experiments are feasible for only a small proportion of the whole root-associated fungal diversity (see Singh, Millard, Whiteley, & Murrell, 2004). An indirect way of assessing the role of interspecific interactions in determining community structure from field observations is to evaluate patterns of co-occurrence (Gotelli, 2000). Studies of co-occurrence in root-associated fungal communities have found co-occurrences to range from predominantly positive (Gorzelak, Hambleton, & Massicotte, 2012;Kennedy, Nguyen, Cohen, & Peay, 2014) to predominantly negative (Koide, Xu, Sharda, Lekberg, & Ostiguy, 2005;Pickles, Genney, Anderson, & Alexander, 2012;Pickles et al., 2010;Wubet et al., 2012). However, and as noted in the studies cited, co-occurrence patterns do not solely reflect interactive associations, but also environmental filtering. If the effect of the environment is not taken into account, positive co-occurrences may reflect that species tend to occur under similar environmental conditions, whereas negative co-occurrence may reflect that species sort according to different environmental niches. To formulate more reliable hypotheses on species interactions based on co-occurrence patterns, it is thus crucial to determine which co-occurrences cannot be explained by the environmental variation (Ovaskainen, Hottola, & Siitonen, 2010).
In this study, we investigated co-occurrences among root-associated fungi, asking whether mycorrhizal and endophytic fungi co-occur with each other more or less strongly than expected based on their responses to environmental conditions and the host plant species. For this, we generated DNA-based community data on root-associated fungi of five plant species evenly sampled along an elevational gradient at a high arctic site (Zackenberg Valley, Northeast Greenland). We hypothesized that co-occurrence patterns among arctic root-associated fungi would be predominantly positive, reflecting facilitative interactions. In particular, we ex-

| Study design
Endophytic fungi were expected to be found in all plant species (Strobel, 2018).
In July 2015, the plant individuals were sampled at 18 locations randomly distributed along the elevational gradient. The lowest sampling location was located at 33 m a. s. l. and the highest sampling location at 479 m a. s. l. which was the upper limit for the distributions of the focal plants. The sampling locations were separated by an average distance of 2.7 km, with the shortest distance between two sampling locations being 373 m and the longest 6.4 km. Within each sampling location, we uprooted five plant individuals from each of the focal species along a transect of at most 50 m. The plant individuals were selected by moving from a starting point towards the southeast along an elevational isocline. Sampled individuals of different species were sampled at any distance, but individuals of the same species were separated by at least 1 m from each other. The whole root system of each plant individual was uprooted, and the fine roots (<2 mm) were collected as samples for the forthcoming molecular analyses, which were initiated in September 2015. The fine root samples were hand-cleaned for soil particles, first in the field and then more thoroughly in the laboratory. Within the next 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. The following environmental variables, expected to shape niches of plants and/or fungi (e.g., Erlandson, Savage, Cavender-Bares, & Peay, 2015), were directly measured at each sampling location: soil pH (in soil-water suspension using a Direct Soil Measurement pH Portable Meter, Hanna Instruments), soil water content (using a HydroSense Handheld Soil Moisture Sensor, Campbell Scientific), the depth of the active soil layer (measuring the distance until the frozen horizon with a metal bar) and vegetation cover (visually estimated cover percentage of all vascular plants in an 1 × 1-m area). These measurements were averaged over three points for each sampling location.

| DNA extraction and sequencing
The root samples were first weighed and then ground into fine powder using a ball mill (Retsch Mixer Mill MM400). Ten milligrams of each sample was used for DNA extraction (all available material was used if less than 10 mg). DNA was extracted using the NucleoSpin Plant II kit (Macherey-Nagel).
As fungal DNA marker, we used the internal transcribed spacer (ITS). The ITS region shows marked variation among fungal species but less variation within species and it is thus considered the universal genetic barcode for fungi (Schoch et al., 2012). We amplified the internal transcribed spacer region 2 (ITS2) using the forward primer fITS7 (GTGARTCATCGAATCTTTG) (Ihrmark et al., 2012) and the reverse primer ITS4 (TCCTCCGCTTATTGATATGC) (White, Bruns, Lee, & Taylor, 1990). We note, however, that the ITS2 is not an optimal marker for discriminating among arbuscular mycorrhizal taxa (Krüger, Krüger, Walker, Stockinger, & Schüßler, 2012), and thus this group of root-associated fungi is likely to be underrepresented in our data.
Before polymerase chain reaction (PCR), all the DNA extracts were diluted to 0.5 ng/µl with ddH 2 O based on NanoDrop Lite (Thermo Scientific) measurements. PCR amplifications were performed using primers fITS7 and ITS4 tagged with 104 unique identification tags following Clemmensen, Ihrmark, Durling, and Lindahl (2016). PCRs were run for 22-35 cycles in a total volume of 50 µl.
The number of cycles was adjusted on a sample-by-sample basis to yield bands of approximately the same strength for all samples on the agarose gel. PCR products were cleaned using the AMPure kit (Beckman Coulter), and DNA concentrations were determined using the Qubit dsDNA HS Assay Kit (Life Technologies). All 450 samples were successfully amplified and pooled into six composite samples, which were further purified using the Cycle-Pure Kit (Omega), and verified for quality on a Bioanalyzer (Agilent Tech). Pooled amplicon mixes were sequenced on a PacBio RS II system (Pacific Biosciences) at SciLifeLab (Uppsala, Sweden) using six SMRT cells. The sequence data have been published in the Dryad data repository (Abarenkov, Somervuo, Nilsson, et al., 2018a).

| Bioinformatics analyses
We clustered the sequences to operational taxonomic units (OTUs) with a 98.5% clustering threshold using the SCATA pipeline (Sequence Clustering and Analysis of Tagged Amplicons, http://scata. mykop at.slu.se). Sequences were screened for tags and primer sequences, requiring a 90% match to the primer sequence. Sequences with a mean quality score lower than 20 or containing bases with a score lower than 3 were discarded. The remaining sequences were aligned pairwise, using usearch, and clustered into OTUs by single linkage clustering with 1.5% maximum distance allowed for sequences to enter clusters, homopolymers collapsed to 3 bp, miss match penalty 1, gap open penalty 0 and gap extension 1. Remaining singletons were removed. Species identification of the OTUs was conducted using the probabilistic taxonomic placement method of Protax-Fungi (Abarenkov, Somervuo, Nilsson, Kirk, et al., 2018b).
This method quantifies the probabilities of all possible taxonomic placements for each query sequence. Importantly, the identification probabilities given by Protax-Fungi account for the possibility that some of the reference sequences are mislabelled, or that the species behind the environmental sequences are missing from the taxonomy or that there are no reference sequences for them. Such sources of uncertainty are particularly common among data on the fungal kingdom (Somervuo et al., 2017). As parameterization data, Protax-Fungi used a fungal taxonomy derived from Index Fungorum (Index Fungorum, 2016), and a fungal reference database derived from UNITE (Kõljalg et al., 2013).
Identified fungal taxa were classified as mycorrhizal or endophytic based on the FUNguild database (Nguyen et al., 2016) as well as the expertise of the authors (N.A. and B.L.). Mycorrhizal species were mostly ectomycorrhizal, but also some ericoid and arbuscular mycorrhizal fungal species were included. Fungal species with uncertain interaction ecology, for which different data sources provided contrasting information or that could not be taxonomically assigned (i.e., no hit), were grouped as "unclassified." The molecular data included 2,874 OTUs, out of which 695 were classified as mycorrhizal and 659 as endophytic fungi. Because our interest here was to estimate species associations, OTUs occurring in less than 5% of the sampling units were excluded from the statistical analyses, leaving 231 OTUs (78 mycorrhizal and 63 endophytic, see Table S1).

| Statistical analyses
To estimate the pairwise co-occurrences among the root-associated fungi, we used the joint species distribution modelling approach (Warton et al., 2015) of Hierarchical Modelling of Species Communities (HMSC; Ovaskainen et al., 2017). HMSC is a hierarchical generalized linear mixed model and is thus structured by fixed effects and random effects. This approach allowed us to jointly model the occurrences and abundances of the fungal species identified in each plant individual as a function of the environmental variables, while accounting for the structure of the study design.
Importantly for reaching the aim of this paper, this modelling approach allowed us to estimate to what extent co-occurrence could be explained by the environmental conditions, and to identify species pairs for which statistically supported co-occurrence remained also after accounting for the environmental conditions and the host plant species. This is achieved by comparing the co-occurrences estimated from a model which does not include any explanatory variables (i.e., raw co-occurrences) with the co-occurrences estimated from a model including explanatory variables (i.e., residual co-occurrences). Residual co-occurrences are more likely to reflect interactive associations than raw co-occurrences (Ovaskainen, Abrego, Halme, & Dunson, 2016;Ovaskainen et al., 2010;Zurell, Pollock, & Thuiller, 2018). The structure of the two fitted models is as follows: 1. "Raw co-occurrences model"-In this model, we accounted only for the variable of sequencing depth (continuous variable, log-transformed), which represents observation effort rather than variation in environmental conditions. While in both models the random effect of the plant individual models the species pair that co-occurs more or less often than expected by random, the two models differ in what "by random" exactly means, as this depends on which fixed effects are controlled for in the models. Thus, in the raw co-occurrences model, "by random" refers to the null expectation that root-associated species would be distributed among the plant individuals fully independently of each other. In the residual co-occurrences model, "by random" refers to the environmentally constrained null expectation that, after accounting for the species affinities to the environmental covariates and host plant species, the root-associated species would be distributed among the plant individuals independently of each other.
To account for the zero-inflated nature of the data, we fitted hurdle-type models. For this, we first modelled presence-absences using a probit-link function, and then abundances (sequence counts) conditional on presence with a log-normal model. We fitted the models using the R-package hmsc 3.0 (Tikhonov et al., 2020) assuming the default prior distributions. The settings applied for posterior sampling through a Markov chain Monte Carlo (MCMC) method and the results of MCMC convergence are provided in the Appendix S1.
We considered that a co-occurrence between a pair of fungal species was well supported if the posterior probability of the association being positive or negative was at least 0.9. To assess how different types of fungi varied in the proportion of positive and negative co-occurrences that they established, we computed the proportions of species pairs that showed well-supported positive or negative co-occurrences within and among the three groups of root-associated fungi.
To evaluate the phylogenetic structure in the co-occurrence patterns, we focused on the taxonomic composition of the groups of positively co-occurring species revealed by the residual co-occurrences model (Figure 1b; Table S1). We evaluated the differences in phylogenetic composition in two ways. We first tested whether higher taxonomic levels (families in the case of mycorrhizal species and orders in the case of endophytic fungi because for the latter the family-level taxonomy is unresolved for many species) were underor over-represented among OTUs that show statistically supported co-occurrences. To make the comparison possible, we conducted this analysis for those taxonomic levels that were represented by at least five species. For each family, we counted how many species were included in the four groups, and compared this number with the null expectation based on 10,000 permutations. In our second phylogenetic analysis, we asked whether the groups of positively co-occurring species differed from each other in their phylogenetic composition. For this, we counted the number of species in each family or order, and computed Bray-Curtis dissimilarity to evaluate how different the groups of positively co-occurring species were in terms of phylogenetic composition. We compared this dissimilarity to the null expectation based on 10,000 permutations.
To evaluate the environmental drivers that cause the difference between the raw and the residual co-occurrences, we quantified how much of the community variation explained by the residual co-occurrence model was attributed to the plant species, elevation, soil pH, soil water content, depth of the active soil layer and vegetation cover, using the variance partitioning approach of Ovaskainen et al. (2017). We evaluated the explanatory power of the presenceabsence model by the area under the curve (AUC) statistic and that of the abundance model by the R 2 value, both averaged over the species (Norberg et al., 2019).

| RE SULTS
The number of statistically supported positive raw co-occurrences was higher than the number of statistically supported negative raw co-occurrences (Figures 1 and 2). Taking the effects of the environmental covariates into account, the majority of the positive co-occurrences remained, whereas the majority of the negative co-occurrences disappeared (Figures 1 and 2). Thus, in line with our original hypothesis, we found the root-associated fungal communities of our high-arctic site to be predominantly structured by positive co-occurrences. The co-occurrences that disappeared between the "raw co-occurrences model" and the "residual co-occurrences model" were mostly attributed to plant species and elevation, the other environmental variables explaining a small part of the community variation ( Table 1). The "residual co-occurrences model" had a high explanatory power, with an AUC value of .86 for the presenceabsence model and R 2 of .30 for the abundance model.

F I G U R E 1
Raw and residual co-occurrences among root-associated fungal OTU pairs estimated from the presence-absence model. The co-occurrences are shown only for cases for which the co-occurrence was either positive (red) or negative (blue) with at least 90% posterior probability (the remaining cases are indicated in white). The OTUs have been ordered emphasizing the network structure, so that OTUs within group 1 (including the mychorrizal group M1, the endophyte group E1 and the unclassified group U1) are found especially often together, the OTUs in group 2 (including M2, E2 and U2) are found especially often together, and the OTUs of group 1 are found especially seldom with the OTUs of group 2. The identity of the OTUs displayed in the panels is given in Table S1, based on identification through Protax-Fungi [Colour figure can be viewed at wileyonlinelibrary.com] Interestingly, most of the positive residual co-occurrences were found among mycorrhizal fungi (30%) or between mycorrhizal and endophytic fungi (25%, Figure 2). Thus, especially mycorrhizal fungi co-occurred positively with either other mycorrhizal or endophytic fungi. The first permutation test showed that no family or order was over-or under-represented among the groups of positively co-occurring species (p > .05 for all families or orders, Table 2; Table S2).
The second permutation test, which asked whether the groups of positively co-occurring species were phylogenetically similar, re- Note: The OTUs have been taxonomically assigned using the probabilistic method Protax-Fungi, provided in Table S1. In cases where the family-level taxonomic assignment was unknown, the order-level taxonomic assignment is given. sampling units where both species are present, a high abundance of one of the species did not imply high abundance of the other species.

| D ISCUSS I ON
In both macroscopic and microbial communities, negative interactions have, thus far, gained more attention than positive interactions.
However, community ecologists now increasingly acknowledge the importance of positive interspecific interactions in shaping communities (e.g., Bruno, Stachowicz, & Bertness, 2003;Tirado, Pugnaire, & Eriksson, 2005), also among microorganisms (Elias & Banin, 2012;Freilich et al., 2011;Nadell, Xavier, & Foster, 2008;but see Foster & Bell, 2012). Our results from arctic plant root systems suggests that positive interactions play an important role in shaping microbial communities. By showing that co-occurrences among fungal species are predominantly positive, our results have important consequences for understanding the structure of arctic interaction networks, where associations with symbiotic fungi play a pivotal role in plant nutrient acquisition (Hobbie & Hobbie, 2006).
In line with our main hypothesis, root-associated fungi formed predominantly positive co-occurrences that were particularly prevalent between mycorrhizal fungi and other root-associated fungi.
Some previous studies have that found predominantly positive co-occurrences among root-associated fungi have suggested that such patterns reflect facilitative interactions (Pan & May, 2009;Yamamoto et al., 2014). This possibility seems plausible also in the present study, especially given that mycorrhizal fungi are known to establish facilitative interactions with bacteria (Frey-Klett et al., 2011;Johansson et al., 2004). Some bacteria enhance fungal growth through improved nutrition and reduction of the impact of adversities (Frey-Klett et al., 2007). We propose that in the same way as mycorrhizal fungi exhibit facilitative interactions with bacteria, they may also form facilitative interactions with other mycorrhizal or endophytic fungi. Just as similar associations between the mycorrhizal helper bacteria and mycorrhizal fungi have been found not to be restricted to particular bacterial or fungal taxa (Frey-Klett et al., 2007), we found no phylogenetic structure among fungi-fungi co-occurrence patterns. As one example of facilitative interactions, some fungi may change the host plant's physiology in a way stimulating or benefitting colonization by other fungi (Kennedy et al., 2014).
Another possibility is that rather than facilitative interactions, the positive co-occurrences detected in our study reflect either mycoparasitic or trophic interactions among pairs of root-associated fungi (Lindahl et al., 2010;Tedersoo et al., 2009). Finally, it might be that positive co-occurrences among fungi arise due to their shared helper bacteria, so that the presence of some fungi and their associated helper bacteria may make it easier for other fungi to colonize the plant.
Most studies analysing observational data have concluded that root-associated fungi mainly compete during colonization of plant roots (e.g., Koide et al., 2005;Pickles et al., 2012;Pickles et al., 2010;Wubet et al., 2012). However, studies that have reported such negative co-occurrences have typically been based on samples that do not represent the roots of single plant individuals, and that have applied co-occurrence analyses without accounting for the effects of environmental variation or spatial structure in the data.
Compared to these previous studies, our sampling design and joint species distribution modelling approach allowed us to unveil patterns of co-occurrences that were not explained by the environment. These co-occurrences were scored at the level of the roots of a plant individual (i.e., at the scale at which fungi-fungi interactions are expected to primarily take place). We found that negative co-occurrences were mostly explained by environmental variation, whereas most positive co-occurrences remained unexplained by environmental variation. However, we note that because our approach of estimating species interactions relies on co-occurrence patterns rather than experimental validation, it is highly susceptible to the choice of environmental covariates that are controlled for in the analysis (Dormann et al., 2018). For instance, our models do not include microscale environmental data (e.g., soil conditions changing at the scale of 1 m or less), which may have explained part of the variation in species occurrences that remained unexplained in our data. While a statistical model will always miss at least some of the factors behind the data, given the high explanatory power achieved by our models, we argue that the general patterns captured by our analyses would not drastically change even if we had accounted also for those "unknown" environmental constraints.
We note that one reason why we captured more positive than negative associations may be that co-occurrence data are typically more informative about positive associations than about negative associations. Namely, if two species are rare in the sense that they inhabit only a small fraction of sampling units, it is unlikely that just by chance they would co-occur in the very same sampling units, and thus the data may contain a strong signal on positive co-occurrence of such rare species. While this may partially explain why we recorded a higher fraction of positive than negative co-occurrences in the raw co-occurrence model, we do not consider that this is the case for the residual co-occurrence model. This is because most of the positive co-occurrences remained statistically supported after environmental variation was taken into account, even if the environmental variables explained a large part of the variation in species occurrences. Hence, we consider the dominance of residual positive co-occurrences to reflect ecological signal rather than a statistical artefact.
Previous studies have shown that root-associated fungal communities show marked variation along elevational gradients, due to the changes in vegetation structure and abiotic conditions (Bahram, Põlme, Kõljalg, Zarre, & Tedersoo, 2012;Coince et al., 2014;Jarvis, Woodward, & Taylor, 2015;Matsuoka, Mori, Kawaguchi, Hobara, & Osono, 2016). In line with these results, we found that plant species and elevation were major factors explaining fungal occurrences. The effect of elevation is likely to reflect the effect of temperature, with high-elevation sites being colder than low-elevation sites (Bahram et al., 2012;Jarvis et al., 2015). The effect of plant species reflects the specialization of root-associated fungal communities to the host plant species (Põlme et al., 2018). One environmental aspect that was not taken into account in our study is differences in microhabitat between mycorrhizal and endophytic fungi (i.e., that the two groups colonize different tissues in the roots). Such differences have been previously suggested to explain their positive co-occurrences within root systems (Wagg et al., 2008). Yet, even if mycorrhizal and endophytic fungi colonize different parts of the root, the chemical changes caused by a given fungus are likely to cause systemic rather than local effects, affecting the physiology of the whole plant individual (Bonfante & Genre, 2010).
Given our intriguing finding of widespread positive co-occurrences between root-associated fungi from a high-arctic site, we propose that future studies should focus on assessing the generality of our results for arctic ecosystems and the mechanisms driving fungal co-occurrences and their significance for plant growth. What remains unresolved is whether the patterns found in our study extend to other fungal groups inhabiting other plant growth forms than those considered here. In our study, the focal plant species establish mostly ectomycorrhizal associations, and thus plants establishing arbuscular associations such as herbaceous plants were not studied. Also, effectively all fungal species classified as mycorrhizal were ectomycorrhizal, because to efficiently detect arbuscular mycorrhizal fungi, a different set of primers is needed (Lekberg et al., 2018).
Studies focusing on the mechanisms driving co-occurrence patterns will be particularly relevant for arctic ecosystems, where plant growth is hampered by the short growing season and the scarcity of nutrients. At a more applied level, added proof for the ecological significance of fungus-fungus co-occurrences might have important implications in biotechnology and soil restoration. Some studies have demonstrated that co-inoculation of some specific pairs of mycorrhizal and endophytic fungal species is more beneficial for plant growth than single inoculations (Reininger & Sieber, 2012.
Our study provides a starting point for investigating fungal species pairs with beneficial effects on plant fitness.

ACK N OWLED G EM ENTS
We thank the three anonymous reviewers for the constructive

DATA AVA I L A B I L I T Y S TAT E M E N T
The species and environmental data matrices used in the data analysis as well as R-scripts for reproducing the results are published in the Dryad data repository, https://doi.org/10.5061/dryad.9kd51 c5dp (Abrego, 2020