The study was conducted in the Zackenberg Valley in northeast Greenland (74°28'N, 21°33'W). The Valley is located in the high Arctic, with a mean annual air temperature of -8.6 °C (measured between 1996 and the year of data collection 2013; Hansen et al. (2017)). Precipitation falls mainly as snow and snow depths vary considerably from year to year (Pedersen et al. 2016). The area is ca. 25 km² in size and is covered by a mosaic of tundra vegetation types (Bay 1998).

Ecological research in this region is supported by two important resources. First, the long-term Greenland Ecosystem Monitoring programme documents inter-annual variation as well as long-term trends in plant and animal biodiversity (Schmidt et al. 2019). Second, floristic work in the region is greatly facilitated by a comprehensive DNA barcode library encompassing more than 80% of the local diversity of vascular plant taxa. Established by Wirta et al. (2015), this library comprises the rbcLa and ITS2 gene regions of all plant species for which samples could be acquired. Thus, the Zackenberg Valley offers major advantages for exploring the representation of plant DNA in soil samples, for comparison with plant communities observed aboveground and for establishing plant-microbe associations, based on both types of data.

To assess the linkages between plant community composition and soil microbes and to compare observational and eDNA based methods, we used an existing dataset from 2013 collected by Stewart et al. (2018). In 2013, Stewart et al. (2018) recorded the vegetation composition in 200 plots across a large-scale environmental gradient. The sampling area extends from sea level (7 m a.s.l.) to the upper slopes of the Aucella Mountain at 770 m a.s.l. The sites were selected using a random stratified design (Hirzel and Guisan 2002), with all sites marked by aluminium pegs for future re‐analyses.

Plant communities were examined within a circle of 1 m² at the centre of each plot. Within each circle, a list of all vascular plant species was compiled, using the taxonomy and nomenclature of the Annotated Checklist of the Panarctic Flora (PAF; Elven et al. (2011)). This material is henceforth referred to as ‘observations’ (of the vegetation). For each plot, each plant species encountered was first scored on a binary presence/absence scale and then on a semi-quantitative scale, based on species cover with the following categories: rare (less than two individuals), < 1%, 1–5%, 6–13%, 14–25%, 26–50%, 51–75% and 76–100% (for more information, see Stewart et al. (2018)).

In 2017, we revisited the 200 plots examined by Stewart et al. (2018) to characterise the plant community by molecular means and to characterise the local community of associated microbes. We note that vegetation turnover in the area is very slow (Schmidt et al. 2012) and that the aboveground plant community is thereby likely to have changed only very little between the two field campaigns. Moreover, low temperature promotes DNA preservation in the soil, which provides optimal conditions to assess plant diversity through eDNA metabarcoding of soil samples (Edwards et al. 2018). During this second visit, a soil core containing 5–10 g of soil was sampled within each of the 200 plots. For this purpose, we used a 5-cc sterile syringe which was twisted into the soil by hand to the depth of ca. 5 cm. For each one-m² plot, only one soil core was sampled to minimize the disturbance and placed in a ziplock bag. Samples were stored on ice in the field and frozen at -20 °C at the end of each field day. Samples were then transported frozen in coolers from the field to storage at -20 °C until further isolation and for characterisation of their DNA.

To characterise the abiotic environment, Stewart et al. (2018) focused on four characteristics found relevant for plant growth and vegetation dynamics (Chapin 1983; Ehrenfeld et al. 2005; Aalto et al. 2013): temperature, moisture and soil pH, as well as soil type. For the current study, we used the existing abiotic data measured in 2013, collected by Stewart et al. (2018). Soil temperature at a depth of 10 cm was measured during the monitoring of each plot, using a General digital soil thermometer model 6300 (Secaucus, NJ, USA). Soil moisture was recorded by Fieldscout TDR300 with 10-cm probes (Chicago, IL, USA) and pH by a Fieldscout Soil stick at a depth of 5–10 cm (Chicago, IL, USA). Soil type was scored, based on visual examination of soil texture using the following categories: Clay; Silt; Fine sand; Coarse sand; Gravel; Humus; Peat. Even though snow cover represents an important variable for vegetation in the Arctic (Niittynen et al. 2020), we were unable to include such data at the resolution targeted in this study (1 m²). For detailed methods on the environmental measurements, see Stewart et al. (2018).

For DNA extraction, 250 mg of each homogenised soil sample was extracted and then purified using the Qiagen DNeasy PowerSoil Pro Kit (QIAGEN, Germany). Samples of DNA-free water were included as blank controls of the extraction protocol. To avoid contamination, all laboratory steps were performed in a laminar flow hood, which was wiped with 70% ethanol and cleansed of potential contaminating DNA with one-hour exposure to mid-range UV light each night. Molecular grade (DNA/RNA-free) tubes, pipette tips, PCR plates and water were used in all protocols.

Our PCR amplification protocols followed those of Wirta et al. (2021), with an initial PCR to amplify the targeted gene region and a second PCR to attach unique indexes with Illumina-specific adapters to the targeted regions. The initial amplifications for each sample were done with a total volume of 10 μl, each containing 5 μl MyTaq Red Mix (Bioline, UK), 1.3 μl DNA-free water (Labnet, Finland), 0.3 μl of each primer (10 μM) and 3 μl of DNA extract.

In order to identify plant taxa from the samples, we amplified the marker ITS2 by using the pair of primers tagF_ITS2-F (5’-ATGCGATACTTGGTGTGAAT-3’) and tagR_ITS2-R (5’-TCCTCCGCTTATTGATATGC-3’) (White et al. 1990; Chen et al. 2010). To amplify the gene region rbcLa, we used tagF_rbcLa-F (5’-ATGTCACCACAAACAGAGACTAAAGC-3’) and tagR_rbcLa-R (5’-CGGTCCAYACAGYBGTCCAKGTACC-3’) (Levin et al. 2003; Ivanova et al. 2016).

For fungi, we assessed part of ITS2 marker by using the pair of primers tagF_ITS3_KYO2 (5’-AHCGATGAAGAACRYAG-3’) and tagR_ITS4_KYO3 (5’-CTBTTVCCKCTTCACTCG-3’) (Toju et al. 2012). For bacteria, we amplified a part of 16S (i.e. V4 and V4–5) using tagF_16S_515FB (5’-GTGYCAGCMGCCGCGGTAA-3’) and tagR_16S_806RB (5’-GGACTACNVGGGTWTCTAAT-3’) (Walters et al. 2016).

The PCR cycling conditions were as follows for the first PCR: the initial denaturation was for 3 min at 95 °C, followed by cycles of 30 s at 95 °C (denaturation), 30 s at 47–55 °C (annealing) and 30 s at 72 °C (extension), ending with final extension for 7 min at 72 °C. For each primer pair, we used a different annealing temperature following Wirta et al. (2021). For plants we, used 55 °C, 47 °C for fungi and 50 °C for bacteria. To amplify plant DNA, we used 35 cycles and 28 cycles were used for fungi and bacteria. To minimise initial amplification bias, each reaction was carried out in two replicates. All the amplicons were checked on a 1% agarose gel and imaged with a BioRad imager to check that the reaction had worked and that the DNA and PCR controls were clean. The PCR replicates were then combined for the second PCR, using 1.3 μl of each PCR product replicate. In the second PCR Illumina‐specific adapters and unique dual‐index combinations for each sample was added to each sample (Vesterinen et al. 2018). The second PCR had a total volume of 10 μl, each containing 5 μl MyTaq Red Mix (Bioline, UK), 0.3 μM of reverse primer, 0.3 μM of forward primer and 2.6 μl of the locus-specific combined PCR product from the first PCR. PCR cycling conditions were the same for all gene regions for the second PCR, starting with 4 min at 95 °C (denaturation), followed by 15 cycles of 20 s at 98 °C (denaturation), 15 s at 60 °C (annealing) and 30 s at 72 °C (extension) and ending with final extension of 3 min at 72 °C. The products of the second PCR were pooled per gene region and per 96 sample plate before concentration using an SPRI bead protocol. The concentrated pooled samples were loaded on a 1% agarose gel and run with 90 V for 120 minutes. The target bands were cut under UV light and the pooled sample was cleaned from gel with the PCR and Gel CleanUp Kit (Macherey-Nagel) and diluted in 2 × 20 μl of the elution buffer. The DNA concentration of the cleaned pools were measured with Qubit 2.0 (dsHS DNA Kit, ThermoFisher Scientific). Based on the compatible lengths of the targeted gene regions, pools of 96 samples were combined in equimolar ratios, with addition of 25% PhiX and sequenced in two MiSeq v.3 2 × 300 runs at the Biomedicum Functional Genomics Unit (FuGU) of the University of Helsinki.

The raw sequences for the plant gene regions ITS2 and rbcLa and the fungal gene region ITS were processed by merging R1 and R2 reads using “pear” (version 0.9.6-bin-64; Zhang et al. 2014), with default values. The resulting merged read pairs were then trimmed using “cutadapt” (version 2.9; Martin 2011) with quality threshold 20 (Phred quality score) and a minimum length of 100 bp.

The taxonomic assignment of plant gene regions ITS2 and rbcLa was done using the local reference databases for the vascular plants of Zackenberg, generated by Wirta et al. (2015). The sequences were assigned to taxa using PROTAX following Somervuo et al. (2016). PROTAX was trained for five taxonomic levels (class, order, family, genus, species) for the relevant plant gene regions (Somervuo et al. 2016). For both gene regions, ITS2 and rbcLa, a separate PROTAX model was built following Roslin et al. (2021). For each taxonomic level, we constructed two community matrices at each taxonomic level, using two probability thresholds, 90% (reliable) and 50% (plausible), for establishing reliable identification.

For the rbcLa gene region, the number of sequences reliably assigned to species proved substantially lower than for ITS2 (Table 1). With a lower probability threshold for species-level taxonomic assignment (Pr > 0.5; plausible), 71.4% of the total of ITS2 reads were assigned to plant species (Suppl. material 1: table S1) compared to 46.4% with the higher threshold of Pr > 0.9 (reliable; Table 1). rbcLa showed a significantly lower success than ITS2 with only 0.1% of the total reads being assigned to species even with this lower threshold (Table 1; for more information on how the assignments varied in relation to the two thresholds, see Suppl. material 1: table S1).

The poor taxonomic assignment success observed for rbcLa resulted in data so sparse that we decided to combine it with data from ITS2, thereby using evidence from both loci to establish species presence (henceforth referred as “eDNA”). For the final community matrix and statistical analysis, we used the community matrix, based on the 90% probability threshold with species-level taxonomy. Here, the read count reflected the number of reads that was reliably assigned to a specific taxon with a high confidence (Pr > 0.9) for ITS2, rbcLa or both.

For both bacteria (16S) and fungi (ITS), we used an alternative approach for bioinformatic analysis and taxonomic assignment. Generally, trimming and quality control of the sequences was conducted according to Vesterinen et al. (2018) and Koskinen et al. (2022). Paired-end reads were trimmed and merged using USEARCH with the “fastq_maxee_rate” algorithm with threshold 1 (Edgar 2010). Primers were removed using software “cutadapt” (version 2.9; Martin (2011)), allowing 20% mismatches and with the minimum length set to 100 bp. Primer-trimmed reads were dereplicated using the USEARCH “fastx_uniques” algorithm with option “minuniquesize 10”, after unique reads were denoised using VSEARCH (Rognes et al. 2016) and uchime3_denovo algorithm used to remove chimeric reads. Then, the denoised unique reads were clustered into operational taxonomic units (OTUs) using the VSEARCH “cluster_fast” algorithm with 97% threshold (Rognes et al. 2016). Finally, reads were mapped back to the original trimmed reads to establish the total number of reads in each sample using the USEARCH “usearch_global” algorithm (~ 93% successfully mapped). OTUs longer than 430 bp with at least 97.0% similarity to the reference database were kept.

For all gene regions, a small number of reads was found in the DNA extraction and PCR controls. Hence, we subtracted the maximum number of reads for a negative sample from all the samples for each OTUs. All samples with fewer than 50 reads in total were subsequently removed. For each sample, OTUs with less than 20 reads per were omitted. Finally, from each sample, we removed OTUs representing less than 0.05% of the total number of reads in the respective sample.

Microbial OTUs were taxonomically assigned to genera using the UNITE usearch/utax database for fungi (Abarenkov et al. 2020) and the Silva database for bacteria (version.123; Quast et al. (2013)). To then assign a functional group to each fungal OTU, we used the microeco package in R, matching taxonomic assignment against FUNGuild (Nguyen et al. 2016). For bacteria, we matched taxonomic assignment against the FAPROTAX database (Louca et al. 2016) following the authors’ instructions (Liu et al. 2023), thereby obtaining the functional role of each bacterial OTU.

The DNA sequencing produced a total of 11.7 M raw sequences, of which 10.8 M passed the quality filters and were assigned to the targeted taxonomic groups (plants, fungi, bacteria; for exact numbers and taxonomic assignment success, see Table 1). Rarefaction curves of fungi and bacteria showed that the sequencing effort was largely sufficient (with all samples recovering full microbial communities; Suppl. material 1: fig. S1).

Overall, out of 9904 microbial OTUs encountered, 2503 were assigned to a specific functional group, with 2099 OTUs of 409 bacterial genera and 404 OTUs of 153 fungal genera successfully assigned a functional role. We grouped the functions assigned into presumptive mutualistic (positive) associations and likely antagonistic (negative) biotic associations, as based on the classification and description given by the FunGUILD database and using the literature for bacteria. Genera assigned multiple functions yielding conflicting assignment to antagonistic versus mutualistic groups were classified into a separate group, “neutral” (or mixed).

To characterise species responses to environmental conditions and to each other, we used Hierarchical Modelling of Species Communities (HMSC, Ovaskainen et al. (2017); Ovaskainen and Abrego (2020)). In this multivariate framework, a matrix of taxon-by-sample observations (the Y community-matrix, with entries y_ij for taxon j at plot i) is modelled as a function of a matrix of plots by environmental covariates (the X matrix, with entries x_ik for covariate k at plot i).

We furthermore separated the different types of response data by treating the observation vs. eDNA methods used to identify each taxon as a taxon-and-method-specific trait (summarised in the T matrix, with entries “OBS” for plant species described by Observation and “eDNA” for plants detected, based on DNA. For plants, we drew on the combined evidence from two loci: ITS2 and rbcLa. For fungal and bacterial OTUs, we drew on the loci ITS vs. 16S, respectively. In other words, any plant species i could potentially occur two times in the Y community-matrix, if recorded as present by ITS2 and/or rbcLa (thus, by eDNA; see Bioinformatics) and by human observation, respectively. Any plant species would then be associated with two different trait states. Our key interest was in testing whether the sampling methodology (i.e. the trait state) affected our estimates of species-specific responses to the environment and species-specific associations with other taxa.

We modelled the taxon presence/absence matrix Y with a generalised linear HMSC model with a probit link (Ovaskainen and Abrego 2020). The occurrences of taxa were modelled as a function of environmental conditions (soil pH, temperature and soil moisture). To control for variation in observation effort resulting from variable sequencing coverage, we included log(sequencing depth) as a taxon-specific co-variate. For plant species described by observation, the value of this co-variate was set to 0.

All co-variates were scaled to a mean of 0 and a variance of 1. Soil type (with six levels, see Environmental co-variates) was included as a random effect. To account for spatial autocorrelation, a spatially structured plot-level random effect was also included and this was modelled, based on the Gaussian predictive process for big spatial data (Tikhonov et al. 2020). However, the spatial structure was negligible, as evidenced by no detectable influence on the model’s explanatory power and no detectable impact on the fixed effects of environmental responses. Spatial structure was, thus, discarded from the downstream analysis. Associations between taxa were examined, based on the residual variance-covariance matrix inspected at the plot level. In the model results, the responses of taxa to fixed effects representing abiotic conditions informed us about individual estimates of environmental responses, whereas the residual variance-covariance matrix informed us about their biotic associations.

Since we explicitly wanted to compare our estimates of taxon-to-taxon associations (i.e. any entries in the residual variance-covariance matrix) to a priori knowledge on the type of association to be expected (positive, i.e. mutualistic or negative, i.e. antagonistic), we only included microbial genera assigned to a specific functional group (see Bioinformatics). Moreover, since taxa with a very low or a very high incidence will contain little information on factors affecting their occurrence, we removed species and OTUs that were present at or absent from less than 5% (n = 10) of the plots (n = 200). Thus, a final set of 44 observed plant species OBS, 37 plants detected by eDNA (out of which 19 species were also observed aboveground), 222 bacterial OTUs and 29 fungal OTUs were included in the model.

The HMSC models were fitted with the R-package Hmsc (Ovaskainen and Abrego 2020; Tikhonov et al. 2020). The models were fitted with four chains and 1,000,000 samples each, which we thinned by 4,000 to yield 250 samples per chain and hence 1,000 samples in total. Since our Y matrix was large and included several different taxonomic groups, we prevented overfitting by modifying the Multiplicative Gamma Process Shrinking Prior parameters a1 and a2 (i.e. the parameters setting control on the level of shrinkage for the species association matrix Omega) to 100 and 100, respectively (for a full justification, see Suppl. material 1: text S1).

MCMC convergence was assessed by examining the potential scale reduction factors of the model parameters. The discriminatory power of the probit model was measured by calculating two different metrics, i.e. species-specific “areas under the curve”, abbreviated as AUC (Pearce and Ferrier 2000) and Tjur’s coefficient of discrimination, henceforth Tjur’s R² (Tjur 2012). These two metrics provide alternative measures of discrimination in the context of presence-absence data, i.e. measures of how well a model discriminates between sampling units where a focal species is present and sampling units from which it is absent (Tjur 2009; Jiménez-Valverde 2012; Norberg et al. 2019).

AUC has become the most commonly applied measure for evaluating model fit of presence-absence species distribution models (Pearce and Ferrier 2000; Elith et al. 2006; Liu et al. 2011). In brief, the AUC equals the proportion of cases for which the occurrence probability for a randomly chosen occupied sampling unit is higher than the occurrence probability for a randomly chosen empty sampling unit or (equivalently) the integral of the Receiver Operating Characteristic (ROC) curve (Hanley and McNeil 1982). An advantage of AUC is that it is not based on any single probability threshold for converting the model predictions into presences or absences, but integrates over all possible such thresholds (Fielding and Bell 1997; Manel et al. 2001). Moreover, it is relatively insensitive to the prevalence of the focal species, i.e. to the fraction of sampling units occupied by the species (Manel et al. 2001; McPherson et al. 2004; Franklin et al. 2009).

A more recent measure of discrimination for presence-absence models is Tjur’s R² (Tjur 2009). This metric was developed as an alternative to other coefficients of determination for logistic regression models (Tjur 2009) and introduced in the SDM literature by Ovaskainen et al. (2016). Tjur’s R² is defined as the difference in the average occupancy probabilities between occupied and empty sampling units (Tjur 2009). One advantage of Tjur’s R² is its resemblance to the “traditional” R² of the linear model, allowing its interpretation as the proportion of variance explained by the model. This proportion can further be partitioned into the portions explained by each of the co-variates included in the model (Ovaskainen et al. 2017). As a challenge, values of Tjur’s R² will typically be (much) lower than a traditional R², urging caution in any comparison to traditional R² values. While Tjur’s R² is still less commonly used than AUC, it has recently gained increasing popularity as a measure of model fit in SDM (e.g. Mang et al. (2018); Zhang et al. (2018); Kotta et al. (2019)).

In attempts to provide clear guidelines for how to assess model fit, AUC values have been split into various “performance classes”, with AUC > 0.9 corresponding to “excellent”, 0.8 < AUC < 0.9 to “good”, 0.7 < AUC < 0.8 to “fair”, 0.6 < AUC < 0.7 to “poor” and 0.5 < AUC < 0.6 to “fail” (Araújo et al. 2005). Such rules-of-thumb for evaluating model fit have been adopted as a common standard (e.g. Thuiller et al. (2006); Marmion et al. (2009); Smolik et al. (2010); Gogol-Prokurat (2011); Wang et al. (2020)). By comparison, values of Tjur’s R² have remained less familiar to ecologists. Since this metric generally achieves lower values than either AUC values of presence-absence models or R² values of linear models, it has often been met with scepticism by ecologists seeking a “good” model fit without appreciating the actual information value of the metric.

Both measures of discrimination can be used for two mutually complementary purposes: on the one hand, we can evaluate the model’s explanatory performance, on the other hand, its predictive performance. In the first case, we ask how well the model is able to predict the same data that it was originally fitted to (explanatory power), in the second we ask how well it predicts validation data independent from the training data used for model fitting (predictive power). To compute explanatory power, model predictions were based on models fitted to all data. To compute predictive power, a four-fold cross-validation was performed, in which the sampling units were assigned randomly to three folds and predictions for each fold were based on a model fitted to data on the remaining four folds. Due to long computational times, cross-validation was run with 25% of the number of iterations used in model fitting. To quantify the relative importance of the various drivers structuring the communities, explained variation was partitioned amongst the fixed and random effects included in the model. To examine each taxon’s response to the model covariates, we assessed the beta parameters (i.e. species- or OTU-specific estimates of environmental responses, equivalent to regression coefficients) in terms of their posterior support and direction (positive vs. negative).

Altogether, we detected 57 plant genera across 24 families. Of these, 68% were identified by both eDNA and direct observation. Fifty-one plant genera were detected by observations and 45 genera by eDNA. Almost a quarter of the genera (21%) were only detected by observation, as compared to six plant genera detected by eDNA alone (Suppl. material 1: fig. S2). Genus-level distribution patterns across the mountain slope were highly similar when established by observation and eDNA (Suppl. material 1: fig. S3).

In total, 107 plant species were detected by observations and 82 by eDNA. Of these, 60 species (46% of the total) were detected by both methods and 47 species (36% of the total) were detected by observation only. Of these 47 species, 15 were not present in the reference barcoding database and, therefore, not detectable by eDNA. In contrast, 27 plant species were detected by eDNA alone. Most plant species were, hence, either efficiently detected by both eDNA and observation or not detected at all by eDNA. The higher the coverage of a plant in a plot, the higher the chance that the same plant species was also detected in soil eDNA (Suppl. material 1: fig. S4). Correspondingly, we found that the frequency at which a species was detected by eDNA was lower for plots in which the plant was scored as absent by direct observation (Suppl. material 1: fig. S5).

Multivariate ordination (NMDS) highlighted partial overlap, but also important differences between observed and eDNA-based plant community composition (Fig. 1A). Plant species composition appeared more homogeneous belowground than aboveground (Fig. 1A) and the belowground samples captured only a subset of the variety in plant species composition that was seen aboveground. A Mantel test of two distance matrices – that of pairwise similarity amongst plots in terms of plants described by observation and that of pairwise similarity amongst plots in terms of plants detected by eDNA – revealed a significant correlation (r = 0.08; p-value = 0.03; number of permutations = 999). Species that were detected by both methods represented, on average, between half and three quarters of all species detected within a plot (Fig. 1B).

Figure 1.

Consistency in community composition observed by eDNA vs. observation. Panel A shows a multivariate ordination (NMDS) of plant community composition across the 200 sample plots, based on observations (blue) vs. eDNA (red). Ellipses summarise the within-community variability (75%) of each method. The NMDS was performed on species presence-absence data using the Jaccard distance. The stress value of the NMDS is 0.185; number of permutations = 999, 2 dimensions, R² = 0.87. In Panel B, the pie charts represent the spatial locations of each plot, and the proportions of plant species detected within each plot by observations only (in blue), by eDNA only (in red) and by both methods (in grey).

In total, 62 plant species (Fig. 2) were found in 5–95% of the plots (and were thereby included in subsequent analyses, see Statistical analysis). Of these, 25 species were only detected by observation, 18 species were exclusively detected by eDNA and 19 species were detected by both methods (Fig. 2). For seven genera, both observations and eDNA indicated a similar incidence amongst plots (Fig. 2; for species-specific distribution patterns, see Suppl. material 1: fig. S3).

Figure 2.

Plant species-specific incidence as scored by the two different methods of identification (observation or eDNA). Bar pairs represent the number of plots where a plant species was detected, with separate colours for detection by direct observation (blue) and/or eDNA (red). Species are sorted according to family (see legend on the left).

Overall, the scoring of a species as locally common or rare (i.e. as having a high or low incidence, equalling occurrence in 10 plots) was roughly consistent between methods (Fig. 2). However, the average number of species detected per plot was significantly higher by eDNA than by observation (11.6 species detected by eDNA vs. 9.4 species by observation, on average; Wilcoxon test; p < 0.05; Suppl. material 1: fig. S4). While species “distributions” (i.e. where a species was recorded in the landscape) differed between methods, there were no signs of DNA “seepage” downslope from where species were observed (Suppl. material 1: fig. S3). For spatial patterns in species distribution amongst methods of identification and how taxon-specific detectability varied across methods, see Suppl. material 1: fig. S3.

The HMSC model was successfully fitted to the data. MCMC convergence was satisfactory and the potential scale reduction factors were close to the theoretical optimum of one (see Suppl. material 1: fig. S6A). The model achieved good discriminatory performance with a mean Tjur R² of 0.30 and a mean AUC value of 0.89. However, explanatory power (Tjur R²) was highly variable amongst taxa (Suppl. material 1: fig. S6B) and the observational data revealed a stronger signal-to-noise ratio than did the eDNA data (OBS mean Tjur R² = 0.20, AUC = 0.84; eDNA mean Tjur R² = 0.08, AUC = 0.75; Suppl. material 1: fig. S6B). The signal was stronger for bacteria than fungi (Bacteria: mean Tjur R² = 0.38, AUC = 0.92; Fungi: mean Tjur R² = 0.17, AUC = 0.85; Suppl. material 1: fig. S6B). Predictive performance assessed by cross-validation was lower than explanatory performance (mean cvTjur R² = 0.08, mean cvAUC = 0.69). On average, it was highest for plants scored by direct observation (mean cvTjur R² = 0.1, cvAUC = 0.72) followed by bacteria (mean cvTjur R² = 0.09, cvAUC = 0.71), fungi (mean cvTjur R² = 0.04, cvAUC = 0.64) and lowest for plants scored by eDNA (mean cvTjur R² = 0.03, cvAUC = 0.56; Suppl. material 1: fig. S6B).

The best predictor of the presence/absence of plant species as described by eDNA was the random effect of the site (accounting for 3% of the total variance; Suppl. material 1: fig. S6B) followed by read count, soil moisture and, equally, soil temperature and the random effect of the soil type (accounting for an average of 2.3%, 1.2% and 0.6%, respectively; Suppl. material 1: fig. S6B). Similarly, for plants described by observation, the random effect of the site accounted for the highest proportion of variation, followed by soil moisture (explaining 4.6% of the total variance), random effect of the soil type (explaining 3.8% of the total variance) and soil pH (explaining 1.9%; Suppl. material 1: fig. S6B). Beyond the random effect, we found that sequencing depth (i.e. read count) accounted for the largest proportion of the total variation amongst bacteria (11.8%) and for a significant part of the variance amongst fungi (3.3%; Suppl. material 1: fig. S6B). Amongst environmental properties, soil pH and soil moisture accounted for 4.9% and 3.8% of the total variance for bacteria and fungi, respectively (Suppl. material 1: fig. S6B). Soil temperature and the random effect of the soil type moisture accounted for 2.7% and 1.7% of the total variance for fungi and bacteria, respectively (Suppl. material 1: fig. S6B).

At the level of individual plant taxa, we found posterior support for a significant directional effect of sampling method on taxon-specific responses to environmental variables (Suppl. material 1: fig. S7A), with stronger responses for data scored by direct observation (Suppl. material 1: fig. S8). When scored by direct observation, a scant majority of plant species showed a negative response to increasing soil pH, whereas when scored by eDNA, no species showed any statistically detectable response to pH (Suppl. material 1: fig. S8). In contrast, plants scored by eDNA exhibited stronger and more positive responses to soil temperature than plants scored by observation. Amongst fungi and bacteria, no individual species responded to soil temperature, whereas both fungal and bacterial species showed a negative response to increasing soil temperature, with stronger responses detected amongst fungi than bacteria (Suppl. material 1: fig. S8).

Relative consistency in species-level responses within organism groups was reflected in general gradients in predicted species richness along environmental gradients. (Suppl. material 1: fig. S9). Overall, plant species richness was predicted to increase across gradients of soil temperature and soil moisture, with generally consistent patterns across data scored by direct observations and eDNA (Fig. 3). For soil pH, both scoring methods suggested no response or a slight decrease (Fig. 3). Across organism groups, predicted fungal richness showed a pattern deviating from that in other organism groups, with a pronounced decrease in richness with increasing soil temperature (Fig. 3).

Figure 3.

Predicted numbers of plant species and fungal and bacterial orders in relation to each environmental co-variate studied a soil temperature b soil pH and c soil moisture. The lines show the predicted relationship, the shaded areas the 95% credible intervals of the predicted relationship. Solid lines represent trends for which there is strong posterior probability (P > 0.95 or P < 0.05). “OBS” stands for plants detected by direct observation, “eDNA” for plants detected by metabarcoding of loci ITS2 and rbcLa, “ITSF” for fungi detected by metabarcoding of locus ITS in fungi and “16S” for bacteria detected by metabarcoding of locus 16S.

To evaluate biotic associations, we first examined residual associations amongst plant taxa across the sample plots, i.e. patterns of co-occurrence unexplained by local soil properties. Plant species-to-species associations were more frequently detected by data generated by direct observations than through eDNA (20% vs. 6% of all possible pairwise associations, respectively). From the subset of plant-to-plant species described by both methods, we found that 75% of the statistically supported interspecific associations detected when plants were recorded by direct observation were not supported when the same species were recorded by eDNA (Fig. 4). The decrease in supported associations primarily concerned positive associations, with a decrease of 82.6% and a decrease of 25% for negative associations detected (Fig. 4). The signs of residual associations detected between plant pairs were generally consistent between the two methods of scoring, with no direct reversal of the estimated sign of the association (Fig. 4).

Figure 4.

Estimated pairwise residual associations amongst plant species at the plot level. Species are ordered by taxonomic groups and by the method of detection (“OBS” for direct observation, “eDNA” for metabarcoding, based on ITS2 and rbcLa). Estimates of taxon-to -taxon associations focus on the 19 plant species identified by both methods. Each plant species is shown as a line and a column within its taxonomic group, with taxa sorted in the same order. Filled cells indicate species pairs showing association with at least 90% posterior probability, with blue for negative associations, red for positive associations.

Pairwise associations amongst plant species and microbial taxa were largely consistent when characterised by direct observation and eDNA. Amongst the large number of taxon-pairs, the majority received no statistical support after accounting for environmental impacts. Associations between plant species and fungi were more frequently detected by data generated by direct observations than through eDNA (27% vs. 13% of all possible pairwise associations, respectively; Suppl. material 1: fig. S10A).

Amongst the pairs showing statistical support, the signs of the associations were largely consistent with a priori expectations. Amongst plants and fungi that were anticipated to have mutualistic interactions (i.e. ectomycorrhizal fungi, endophytes and epiphytes symbiotroph which receive nutrients by exchanging resources with host cells), we found three times more positive (17.5%) than negative associations (5.3%) when the plants were scored by direct observation, whereas for eDNA, the number of negative associations (5.3%) equalled that of positive associations (Fig. 5).

Figure 5.

Residual associations detected amongst plants and different functional groups of fungi. Here, each plant species is shown as a column, including only the 19 plant species identified by both direct observation and eDNA and, thus, allowing direct comparisons. Rows correspond to individual microbial genera, as sorted by functional groups. Red fields indicate antagonistic relationships, blue fields mutualistic interactions and grey fields indicate neutral or mixed interactions (i.e. the same genus being associated with several different functions). For visual comparison, each cell is divided into two, with the upper part describing the association estimated when plant occurrence was detected by eDNA and the bottom part describing the association estimated when plant occurrence was detected by Observation. G corresponds to the functional group, F to the Fungal taxon and P to the plants taxon. For the identity of individual taxa, see key in Suppl. material 2.

Amongst plants and fungi forming presumptively neutral or mixed interactions (i.e. lichens, fungi with imprecise functions, insect pathogens or saprotrophs that receive nutrients by breaking down dead host cells), we found a higher proportion of positive associations (16.5%) than negative associations (8.3%) when plants were scored by direct observation (Suppl. material 1: figs S10B, S11A). When plants were detected by eDNA, we found a lower number of statistically supported associations. Yet, only one discrepancy between scoring methods was observed in the sign of the association (that between Silene acaulis and one fungal genus F15 in the Heliotaceae family; Fig. 5).

Amongst bacterial genera and plants assumed to engage in antagonistic associations (i.e. plant pathogens or intracellular parasites), the associations estimated were fully consistent whether characterised, based on observational or eDNA data (Fig. 5). Unexpectedly, within this group, we found a preponderance of positive (19.8% of taxon-pairs) over negative (9.6%) associations (Fig. 5). For presumptively mutualistic interactions (i.e. bacteria involved in nitrification or denitrification, as indispensable for plants nutrition), we estimated a dominance of positive associations (Fig. 5), with only one discrepancy between the two methods of scoring plants (for Silene acaulis vs. bacterium B252 in Sphingobacteriales; Fig. 5). Amongst plants and bacteria forming presumptively neutral associations (i.e. bacteria involved in hydrogen oxidation, fermentation or ureolysis that can affect the carbon cycle), we found no conflict between plants scored by direct observation or eDNA (Suppl. material 1: figs S10B, S11B).

In terms of diversity of plants at the landscape scale, the traditional scoring of plant presence by direct observation detected more plant species than did scoring through eDNA. However, the number of plant species detected per plot was, on average, slightly higher when using eDNA than when using direct observation of plant individuals. Plant prevalence (i.e. the number of times a plant was recorded across the landscape) was also typically higher when assessed using eDNA rather than direct observations. This suggests that, per plot, approaches based on eDNA may be more sensitive in detecting plant species than are aboveground surveys. This difference may be due to plant DNA commonly being present in the soil even in the absence of the organism in its life stage(s) that can be readily identified (Taberlet et al. 2012; Deiner et al. 2017). Such differences in the plant communities perceived by different methods of scoring may be caused by, for example, the presence of a seed bank or by clonal reproductive structures in the soil. The belowground DNA from roots and tissues is also typically the least degraded and, therefore, might detect past signals up from 30 years (Foucher et al. 2020; Ariza et al. 2023). However, at this point, we cannot distinguish the relative contributions of DNA from four potential sources: plants undetected by observation; plants lying dormant in the soil; locally extinct plants conserved in the soil; and/or plant parts (e.g. seeds or pollen) potentially transported from elsewhere. False positives due to species misidentifications are also typical of data derived from DNA-based methods and represent a topical question in this field of research (Ficetola et al. 2015; Buxton et al. 2022). Despite these shortcomings, we were able to rule out one potential source of exogenous DNA in soil samples, as we detected no signs of DNA for specific species downslope from where they were observed (Fig. 1B, Suppl. material 1: fig. S3). Thus, the transport of DNA with seepage water (Pedersen et al. 2015; Barnes and Turner 2016) seems not to muddle plant distributions in our study region.

Our success rates in the eDNA-based taxonomic assignment and the overlap of the taxon lists between eDNA and the plants described by observation were similar to those in other studies in northern ecosystems (Vasar et al. 2023). The match in detection probabilities was particularly good at the genus level. A large proportion of all plant genera detected by observation was also detected via eDNA. This adds credence to eDNA-based tools as sensitive descriptors of the vegetation – which is hope-inspiring, since while expert botanists will have little trouble in characterising an Arctic vegetation plot, such botanists are in short supply. Moreover, the time during which plants can be reliably identified at the high latitude of our target community is only some months per year, making the characterisation of ecologically-informative numbers of plots a true challenge. Quick sampling of the local community by a single syringe-full of soil is then an attractive alternative.

Importantly, the discrepancy between plant data scored by observation and data recovered by eDNA became much more pronounced when examined at the species level. Here, a large proportion of taxa detected by observation were lacking from the species lists generated by eDNA-based tools and these false absences were particularly centred on selected taxa, such as Cyperaceae or Poaceae (Fig. 2). While the absence of 15 plant taxa could be attributed to an incomplete reference library, most of the species detected by observation in our study (93 species) were also represented in the reference database (Wirta et al. 2015). Thus, database bias is likely to be minor in our case. A main reason for the contrast between genus- and species-level consistency is likely to be found in the probabilistic taxonomic assignment method employed, where only assignments with a confidence level of at least 90% were accepted. In taxa such as grasses, the occurrence of multiple related species in the Zackenberg flora might have meant that no single species passed this threshold – even when we were reasonably sure about the genus. For this reason, a much larger number of sequences was assigned to genera than to species (Table 1). This interpretation is supported by the higher taxonomic assignment success achieved at species-level when applying a lower confidence threshold (i.e. 0.5; Suppl. material 1: table S1).

As another important methodological consideration, different markers will provide widely differing levels of taxonomic detection. This was illustrated by remarkable variation in taxonomic assignment success when using locus rbcLa for taxonomic assignment to different levels. Here, assignment success dropped markedly from the genus to the species level (Table 1). This is only to be expected, given the poor taxonomic resolution provided by this gene region for Zackenberg plants (Wirta et al. 2015; Foster et al. 2021). Multiple studies suggest that the gene region amplified by rbcLa universal primers contains little species-level variation, resulting in low discrimination power amongst plant species (Costion et al. 2011; Foster et al. 2021; Trujillo-Argueta et al. 2022). In contrast, at a larger scale and at higher taxonomic levels, rbcLa has been reported to successfully identify all families and nearly all genera, with species identification rates varying significantly amongst plant groups (Dong et al. 2014).

For the current flora of our high-Arctic study site, we suggest that the marker ITS2 alone may be sufficient to produce reliable data on species’ distributions across the landscape.

In terms of species responses to abiotic drivers, our data revealed ample ecological signal in samples from across the landscape. Nonetheless, the method of identification generated some variation in estimates of the intensity of the abiotic imprints on individual species and predictions of overall response. For plants described by observation, soil pH emerged as an important factor, with predominantly negative responses to increasing pH. However, when plants were scored by eDNA, we found no matching pattern. Plant responses to increasing soil moisture were more consistent across methods.

Beyond the detailed descriptions of the various plant species’ ecological niches (above), the same soil samples also sufficed to characterise the distributions of microbial taxa. When examining the abiotic imprints on our different organisms in further detail, our results reveal some interesting taxon-specific patterns. Our results indicated that, in general, bacterial and fungal community diversity patterns were modulated by distinct ecological drivers. To summarise the relevant discrepancies amongst kingdoms, soil moisture appeared to be a factor relevant to all groups of organisms. This is only to be expected, as high-Arctic ecosystems are typically arid or semi-arid, where soil water is a limiting resource for plant growth as well as for microbial activity (Griffiths et al. 2003; Liu et al. 2010; Nabe-Nielsen et al. 2017). However, bacterial responses to moisture varied between being positive and negative and, hence, these patterns did not translate into any overall trend in predicted taxonomic richness along the soil moisture gradient. Previous studies have shown that bacterial communities are sensitive to water addition and soil moisture (Kaisermann et al. 2015; Umair et al. 2020; Jaeger et al. 2023). Importantly, such responses may also translate into knock-on effects on nutrient availability. In this study, we found that a majority of bacteria that showed a positive significant relationship with soil moisture were associated with denitrification and/or nitrification (Fig. 5). This result was well in line with previous literature showing that, in high-Arctic wet sedge meadows and heath tundra soils, denitrification is mainly controlled by soil moisture (Siciliano et al. 2009; Blaud et al. 2015).

Amongst fungi, soil temperature proved a strong predictor of the occurrences of a third of the individual taxa and a driver of overall fungal richness (Fig. 3; Suppl. material 1: fig. S8). Bacteria, in contrast, showed no overall richness trend with soil temperature (Fig. 3). This finding agreed with previous studies showing that, indeed, fungal activity is more sensitive than bacterial activity to soil temperature (Suppl. material 1: fig. S8). Overall, fungi seem to be more adapted to low-temperature conditions than are bacteria (Pietikäinen et al. 2005; Bárcenas-Moreno et al. 2009). As a consequence, bacterial community composition may be dependent on long-term trends in near-ground temperatures and soil moisture regimes (Frindte et al. 2019). However, the effects of fine-scale heterogeneity and potential interactive effects between temperature and soil moisture regimes on microbial taxa remains unclear (Perez-Mon et al. 2020). In soils subjected to frequent freeze-thaw cycles, such as those typical of Greenland, we tend to find microbial communities better adapted to variation in temperature and moisture (Perez-Mon et al. 2020).

Our results suggested the local fungi to show little response to soil pH. This finding contrasts with previous studies, which identified pH as one of the main determinants of elevational diversity patterns amongst many bacterial and fungal communities across Arctic and similar ecosystems (Timling et al. 2014; Shen et al. 2015; Zhang et al. 2016; Canini et al. 2019). In a similar vein, Fierer and Jackson (2006) pinpointed soil pH as the strongest predictor of bacterial richness, diversity and overall community composition. Likewise, Canini et al. (2019) have emphasised how pH affects the richness of ericoid mycorrhizal fungi, root associated fungi, plant pathogens, animal pathogens and mycoparasites, while ectomycorrhizal fungi richness seemed less affected. Matching such findings are reports of indirect effects of phylogenetic structure on microbial responses to pH (Zhang et al. 2016; Liu et al. 2018).

A few factors may account for the lack of detectable response in the current setting. Overall, we detected relatively few plant pathogens, animal pathogens and ectomycorrhizal fungi amongst the taxa assigned to a clear-cut function. Moreover, the pH effect of previous studies might, in some cases, reflect a cofounding effect, where young soils after snow melt contain only low amount of organic matter and clay minerals, resulting in higher pH than older soils covered first by biological crusts and then by shrubs (Kwon et al. 2015). Such variation may significantly impact the composition and structure of plant communities, which then act as a major driver of variation in fungal communities (Krüger et al. 2017; Maciá-Vicente et al. 2023). Separating between cause and effect can then be a complex issue. Nonetheless, the perhaps most likely reason for the mixed and conflicting patterns detected is simply the limited range of pH variation observed within our study area. Most soils showed a relatively high pH and the full range of variation was only pH 5.06–8.57.

In conclusion, the mixed responses of different organism groups to different environmental gradients supported previous reports of decoupled ecological responses to environmental conditions amongst bacteria and fungi (Lange et al. 2014; Liu et al. 2020; Shen et al. 2020). Such variation will result in a mosaic of organism communities across heterogeneous landscapes, with potential knock-on effects in terms of ecological functioning due to the varying contributions of different functional groups.

In terms of plant species’ associations with other plants, the overall signs (positive or negative) of estimated residual plant-plant associations were largely consistent between sampling methods. Despite some differences in the absolute numbers of associations detected by each survey method, the overall association patterns were similar in terms of the proportions of negative and positive associations amongst groups of taxa. Nonetheless, we observed one contradictory association between Salix arctica and Silene acaulis between methods and some discrepancies in terms of which plants were associated with each other. Slight variation in the set of residual associations can be explained by the fact that the species prevalence was higher when described by eDNA than by observation. Ultimately, two plant species will, thus, have had a higher chance of being found within the same plots when assessed by eDNA than by direct observation. Therefore, the model had access to more data on species co-occurrence for the former than the latter. An important alternative interpretation is still that some apparent associations arose from species’ responses to unmeasured environmental gradients, such as soil chemical properties beyond pH or to snow conditions (Niittynen et al. 2020; Rixen et al. 2022), which we were unable to measure in the current setting. This interpretation was supported by the fact that close to 8% of variance in plant species occurrences (based on both observations and eDNA) was captured by a spatial random effect, indicating the presence of spatial structure in plant distributions beyond that captured by the current model covariates. Clearly, the current data on species co-occurrences alone were insufficient to ultimately prove biotic associations. Rather, experimental tests are required to assess if some of these residual associations arise from direct biotic interactions (Blanchet et al. 2020). For now, a prudent interpretation is that the patterns observed point to potential interactions. Positive associations amongst plant species may then reflect competition or facilitation between plants or that the two species prefer the same unmeasured environmental conditions. A negative association could reflect exclusive competitive or that two taxa are found in significant different environments and are, therefore, found less often than by chance.

Overall, our results reveal that a fully DNA-based approach will suffice to generate a wealth of data-driven hypotheses of not only species responses to the abiotic environment, but also to each other. A wholesale characterisation of species communities, including both the responses to abiotic imprints and the co-association within the communities, is then a truly exciting perspective.

Beyond plant-plant associations, we observed a high number of residual associations between plants and microbes. Amongst these, statistically supported associations between plants and bacteria were more frequent than between plants and fungi. Classifying the microbes into functional groups allows us to speculate about ecological function. Consistent with our a priori assignments, we found a majority of the statistically-supported positive associations between plants and fungi to match with presumptively mutualistic associations – but opposite to our expectations, the majority of positive residual associations with bacteria corresponded to relationships that we expected to be antagonistic.

The variable signs of the associations uncovered reflected the challenges involved in estimating processes from patterns. While a mutualistic association may be detectable as mutual attraction – and thereby higher co-occurrence than expected, based on the joint environment – it is well conceivable that antagonistic interaction partners may likewise show positive associations. To understand why this may be the case, consider the case of a classical predator and its prey. The two frequently co-occur due to the simple reason that predators will accumulates in areas of dense prey. One may then observe a positive association, despite the obvious antagonistic effect of the predator on individual prey. In our case, the antagonistic associations were associated with functional groups that could reflect direct negative impact (i.e. intracellular parasites or predatory and exoparasitic taxa). Within the antagonistic functional groups, we found bacteria belonging to Clostridium, Nocardia, Rhodoplanes, Sphingomonas and Streptomyces genera, all of which are known to dominate the plant-associated communities of Alpine, Arctic and Antarctic Regions (see Cripps and Eddington (2005)). Members of these bacterial taxa have been shown to be cold-adapted and tightly associated with plants, supporting their potential importance for plant fitness and survival (Almario et al. 2022).

Amongst mutualistic associations, we included bacteria connected to soil nitrogen fixation, nitrification, denitrification, ammonification and other major nitrogen transformation processes mediated by soil bacteria. Such processes include the oxidation of aerobic ammonia, nitrate reduction, aerobic nitrite oxidation and ureolysis, which can be involved in denitrification (based on FunGuild classification; Nguyen et al. (2016)). While such processes may increase nitrogen availability, microbes can also reduce ecosystem nitrogen availability by transforming nitrogen to more mobile forms such as nitrate (van der Heijden et al. 2008; Kuzyakov and Xu 2013). Hence, nitrifying bacteria can indirectly reduce plant productivity, as they reduce the availability of a nutrient, which limiting plant productivity (van der Heijden et al. 2008). Ecosystems such as Arctic and Alpine tundra are most likely to be strongly nutrient limited and here soil microbes have been shown to compete effectively with plants for nitrogen (Nordin et al. 2004; Almario et al. 2022). Hence, part of the positive associations detected between plants and mutualistic bacteria can potentially result in competition for nutrients in soil solution, with possible negative effects on plant nutrient acquisition and growth. All in all, further characterisations of plant-associated bacterial communities in Arctic soils are, thus, needed to better clarify the key drivers and outcomes of associations between plants and bacteria.

Amongst plants and soil fungi, we observed predominantly positive residual associations. This finding brought important insights for understanding the structure of Arctic interaction networks. Amongst the fungi classified as mutualistic, we included only taxa forming mycorrhiza and taxa known to be epiphytic or endophytic symbiotrophs. The fungal taxa assigned to ectomycorrhiza were Cortinarius sp. (Agaricales) and Cenococcum sp. (Mytilinidales). Amongst associations with the 19 plant species described by direct observation, we found more than 60% of the statistically-supported associations to be positive. Both Cortinarius and Cenococcum were positively associated with the same plant species, namely Minuartia rubella, Trisetum spicatum, Silene acaulis, Saxifraga oppositifolia and Saxifraga cernua. Cenococcum showed a positive association with Stellaria longipes and Salix arctica, while Cortinarius did not. The epiphytic symbiotroph fungus Knufia was found to be positively associated with Minuartia rubella and Saxifraga oppositifolia. Interestingly, previous studies have reported Minuartia rubella and Silene acaulis to be associated with either arbuscular mycorrhiza or to be devoid of any mycorrhizal associations (see Cripps and Eddington (2005)). Similarly, Trisetum spicatum, Saxifraga oppositifolia are generally associated with arbuscular mycorrhiza (Cripps and Eddington 2005). However, Salix arctica has been recorded as being associated with both ectomycorrhiza and arbuscular mycorrhiza (Newman and Reddell 1987; Gardes and Dahlberg 1996; Kytöviita 2005). These findings supported previous suggestions that low specialisation by mycorrhizal fungi on the plant species of the Arctic could be an adaptive response to low nutrient availability, thereby ensuring nutrient uptake in nutrient-poor environments (Botnen et al. 2014; Abrego et al. 2020). However, from a methodological viewpoint, it is important to emphasise that the fungal ITS2-primers might provide a skewed picture of the fungal diversity due to primer biases (Tedersoo et al. 2015; Botnen et al. 2020). As a result, the detection of primarily ectomycorrhizal taxa in the current study may in part be an artefact, as a different set of primers will be needed to efficiently detect arbuscular mycorrhizal fungi (Lekberg et al. 2018; Rasmussen et al. 2022).

To the group of endophytic-symbiotroph fungi, we assigned two fungi from the genera Phialocephala and Cadophora, both belonging to order Heliotales. These fungi were engaged in predominantly positive associations with selected plant species: Minuartia rubella, Hierochloe alpina, Cassiope tetragona, Bistorta vivipara, Saxifraga cernua, Salix arctica, Poa arctica and Pedicularis hirsuta. Although both mycorrhizal and endophytic fungi may have positive effects on plant fitness, they colonie plant roots differently (Peterson et al. 2008; Andrade-Linares and Franken 2013; Yan et al. 2019). Plants and mycorrhizal fungi communicate chemically during the root colonisation proces, and, to some extent, plants are able to regulate their level of mycorrhizal colonisation (van der Heijden et al. 2015). Concerning endophytic fungi, they can colonie plant parts other than roots and are generally described as generalist opportunistic coloniers (Knapp et al. 2012; Mandyam et al. 2012; Mandyam and Jumpponen 2015). However, the extent to which plants can regulate the range of mycorrhizal and endophytic fungi differently remains to be validated by experimental approaches.

It has been argued that contrasting diversity patterns between plants and bacteria may emerge as a result of inconsistent sampling approaches amongst studies, rather than being attributable to true ecological mechanisms (Bryant et al. 2008). However, in our case, both positive and negative residual taxon-to-taxon association patterns between bacteria and plants remained similar, regardless of whether the plants were detected by eDNA or direct observation. This added credence to the patterns uncovered.

Overall, we find that ecological patterns of plant-plant and plant-microbial association were largely consistent between direct observations and eDNA-based scoring of plants. Differential responses to abiotic conditions amongst plants and soil microbes (e.g. fungi with soil temperature) reinforced the view that aboveground and belowground communities may present decoupled responses to environmental gradients (Classen et al. 2015; Donhauser and Frey 2018; Fei et al. 2022). We reiterate that the residual associations detected in our study might represent biotic interactions, but that such conjectures can only be confirmed through experimental tests (Saine et al. 2020). Furthermore, it is important to remember that the single soil core examined per plot (1 m²) in our study will capture a highly local signal and that increased sampling effort would be required for larger sampling areas (Edwards et al. 2018). Nonetheless, given the intimate relationships of fungi and plants, our results pointed to widespread associations between fungal and plant communities. The overall dominance of positive associations between plants and the fungi which were a priori considered to be mutualistic and the high proportion of positive associations observed between plants and bacteria considered as antagonistic, highlighted the strong potential role of plant-soil microbe interactions as key drivers of Arctic community assembly.

The authors have declared that no competing interests exist.

No ethical statement was reported.

BP and TR were funded by the Academy of Finland (VEGA, grant 322266 to T.R). TR and OO were also funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC-synergy grant 856506—LIFEPLAN). MMJ was supported by the Academy of Finland's ‘Thriving Nature’ research profiling action.

BP and TR designed the study. LS, LP and NMS surveyed vegetation and produced environmental data. BH and SP collected soil samples. BP and HW generated metabarcoding libraries from soil samples. PS trained PROTAX to the reference library and made the taxonomic assignments for plants. EV performed the bioinformatics and taxonomic assignment for fungi and bacteria. BP, OO, MJ and JS analysed the results. BP, HW and TR wrote the first draft of the manuscript. All authors then helped draft the manuscript or provided comments on the final manuscript.

Bastien Parisy https://orcid.org/0000-0001-6785-8544

Niels M. Schmidt https://orcid.org/0000-0002-4166-6218

Helena Wirta https://orcid.org/0000-0002-4667-2166

Laerke Stewart https://orcid.org/0000-0002-8618-6779

Loic Pellissier https://orcid.org/0000-0002-2289-8259

Sam Pannoni https://orcid.org/0000-0003-4456-4839

Panu Somervuo https://orcid.org/0000-0003-3121-4047

Mirkka M. Jones https://orcid.org/0000-0002-8157-8730

Jukka Siren https://orcid.org/0000-0002-2680-0597

Eero Vesterinen https://orcid.org/0000-0003-3665-5802

Otso Ovaskainen https://orcid.org/0000-0001-9750-4421

Tomas Roslin https://orcid.org/0000-0002-2957-4791

The sequence datasets generated during the current study are available in the Sequence Read Archive repository, in the BioProject SUB12531579.