Community composition of arctic root-associated fungi mirrors host plant phylogeny.

The number of plant species regarded as non-mycorrhizal increases at higher latitudes, and several plant species in the High-Arctic Archipelago Svalbard have been reported as non-mycorrhizal. We used the rRNA ITS2 and 18S gene markers to survey which fungi, as well as other micro-eukaryotes, were associated with roots of 31 arctic plant species not usually regarded as mycorrhizal in Svalbard. We assessed to what degree the root-associated fungi showed any host preference and whether the phylogeny of the plant hosts may mirror the composition of root-associated fungi. Fungal communities were largely structured according to host plant identity and to a less extent by environmental factors. We observed a positive relationship between the phylogenetic distance of host plants and the distance of fungal community composition between samples, indicating that the evolutionary history of the host plants plays a major role for which fungi colonize the plant roots. In contrast to the ITS2 marker, the 18S rRNA gene marker showed that chytrid fungi were prevalently associated with plant roots, together with a wide spectrum of amoeba-like protists and nematodes. Our study confirms that arbuscular mycorrhizal (AM) fungi are present also in arctic environments in low abundance.

INTRODUCTION

Most plants form root symbiosis with fungi, of which mycorrhizal associations are most common (Smith and Read 2008; Brundrett 2009). The ectomycorrhizal (ECM) and ericoid mycorrhizal (ErM) symbioses dominate in northern ecosystems, while arbuscular mycorrhiza (AM) is more widespread in the tropics (Smith and Read 2008; Steidinger et al. 2019). In arctic environments, ECM and ErM fungi are common but limited to a few geographically widespread plant species (Gardes and Dahlberg 1996), while AM has been assumed to be less common. Despite the importance of mycorrhizal symbiosis in comparable nutrient-poor temperate systems, a large portion of High Arctic plants have been reported as non-mycorrhizal (Väre, Vestberg and Eurola 1992; Kytöviita 2005; Newsham, Upson and Read 2009). Newsham, Upson and Read (2009) found evidence for an increase in plant species without mycorrhizal symbioses at higher latitudes. In a microscopy-based survey of root-associated fungi (RAF) of 76 plant species in the High-Arctic Archipelago Svalbard, Väre, Vestberg and Eurola (1992) observed only one single AM spore in a soil sample. Similar results have been reported from several other High-Arctic locations, with no (Bledsoe, Klein and Bliss 1990) or very few occurrences of AM formation (Kohn and Stasovski 1990). However, this is not the case for all location studied; AM structures were frequently found in several plants in the Canadian High-Arctic (Dalpé and Aiken 1998; Olsson, Eriksen and Dahlberg 2004; Allen et al. 2006), and a recent study from Svalbard showed AM structures on multiple plant species (Newsham et al. 2017).

Root endophytes, such as dark septate endophytes (DSE), commonly occur in cold-stressed arctic, sub-arctic, Antarctic and alpine areas on both mycorrhizal and non-mycorrhizal plant roots (Read and Haselwandter 1981; Jumpponen and Trappe 1998; Newsham, Upson and Read 2009; Hill et al. 2019). They may even be the most common type of RAF associations in certain plants (Ruotsalainen, Väre and Vestberg 2002). Studies indicate that DSE can have beneficial effects for the plants (Jumpponen and Trappe 1998; Newsham 2011; Berthelot et al. 2016), such as increased nutrient uptake (Jumpponen, Mattson and Trappe 1998; Upson, Read and Newsham 2009; Hill et al. 2019). However, endophytes are in general less studied than their mycorrhizal relatives, and knowledge on their distribution and ecological roles in the Arctic is limited.

The level of host specificity of root-associated fungi varies between ecosystems and mycorrhizal type. In boreal and temperate areas, it appears to be a higher degree of host-specificity compared to the Arctic, especially for ECM fungi (Ishida, Nara and Hogetsu 2007; Ryberg, Larsson and Molau 2009; Ryberg, Andreasen and Björk 2011; Timling et al. 2012; Bonito et al. 2014; Botnen et al. 2014; Linde et al. 2018). Several studies from arctic and alpine areas have shown that the dominant ECM plants typically share the same pool of ECM fungal symbionts (Ryberg, Larsson and Molau 2009; Ryberg, Andreasen and Björk 2011; Timling et al. 2012; Botnen et al. 2014). However, for non-mycorrhizal arctic plants or plants with other mycorrhizal types, such as AM plants, the host identity can influence the associated fungal community (Becklin, Hertweck and Jumpponen 2012; Fujimura and Egger 2012).

It is estimated that plant roots represent around 80% of the plants’ total biomass in the Arctic (Mokany, Raison and Prokushkin 2006). As such, plant roots represent the largest available pool of carbon for both symbionts, parasites and biotrophic fungi in the Arctic, and may host an extensive part of the total terrestrial diversity. It is, therefore, crucial to investigate this diversity to get a broader picture and a better understanding of the ecosystem as a whole. From other areas, it is known that plant roots may host a rich diversity of other micro-eukaryotes (Durán et al. 2018; Hassani, Durán and Hacquard 2018). These organisms may potentially interact with both plant host and associated fungi, but there has been considerably less focus on these organisms compared to fungi. Some surveys of arctic soils have revealed a high diversity of various protists (Shi et al. 2015; Malard and Pearce 2018), where habitat type seems to be an important structuring factor (Gittel et al. 2014; Shi et al. 2015).

Plant-root-associated microorganisms can be assessed by different means. Microscopy provides direct evidence for whether mycorrhizal functional structures are formed (Väre, Vestberg and Eurola 1992) and much of our basic knowledge about the mycorrhizal status of arctic plants is based on microscopic investigations (Hesselman 1900; Väre, Vestberg and Eurola 1992; Gardes and Dahlberg 1996). However, it can be hard to determine which specific fungi are involved using microscopy. DNA metabarcoding, which is based on high throughput DNA sequencing of amplified markers, is established as an important tool to survey microbial and fungal communities (Lindahl et al. 2013). DNA data does not inform whether mycorrhizal functional structures are formed. Still, DNA metabarcoding can provide a first comprehensive view of which organisms are present and is especially useful in poorly explored ecosystems. In the Arctic, mycorrhizal fungi associated with ECM plants have been widely studied using DNA metabarcoding (Blaalid et al. 2014; Botnen et al. 2014; Mundra et al. 2015a, 2016; Morgado et al. 2016), but to a less extent on tentative non-mycorrhizal plants or plant with unclear mycorrhizal status.

In fungal DNA metabarcoding studies, communities are usually surveyed using the rRNA ITS gene region (Schoch et al. 2012; Lindahl et al. 2013), but some fungal groups are often missed out with this marker because of primer and length biases (Bellemain et al. 2010; Tedersoo et al. 2015). The V4 18S rRNA gene marker is typically used to survey the diversity of micro-eukaryotes (including protists) and general 18S rRNA gene primers also amplify fungal groups often missed out by the ITS marker, such as chytrids, Archeaorhizomycetes and AM fungi (Rosling et al. 2011; Tedersoo et al. 2015 2015; Dunthorn et al. 2017). However, the 18S rRNA gene is a far more conserved marker than ITS, and therefore the18S rRNA gene data cannot be interpreted at lower taxonomic levels (i.e. species, genera).

In this study, we use both ITS2 and 18S rRNA gene markers to survey the diversity of fungi and other micro-eukaryotes in roots of 31 assumed non-mycorrhizal arctic plant species or plants species with unclear mycorrhizal status. We aim to (1) assess which types of fungi (taxonomically and functionally) are associated with roots in common arctic plant species not prevalently regarded as mycorrhizal, and to (2) obtain basic knowledge about tentative root symbionts and their mycorrhizal status. By using the 18S rRNA gene marker we also wanted to (3) obtain information about which other micro-eukaryotes are present, which has largely been neglected in the past. Further, we wanted to (4) assess whether the root-associated fungi are host generalist, as is typically the case for arctic ECM fungi, or whether they show some level of host specificity. In this regard, we (5) investigate whether the phylogenetic distance between the host plant species reflects which fungi are colonizing specific plant species.

MATERIALS AND METHODS

Sampling site

Fieldwork was conducted during July 2013 and 2014 in the High-Arctic Archipelago Svalbard, Norway (Fig. 1 and Table S1, Supporting Information; RIS-ID: 6549 and RIS-ID: 6165). We sampled 31 plant species assumed to be non-mycorrhizal or having unclear mycorrhizal status (1–7 root systems per species) from 11 different locations across the archipelago (Fig. 1).

Figure 1.

Map with overview of sampling locations in Svalbard.

Sampling procedure, DNA extraction, PCR procedure and Illumina sequencing

Whole root systems were dug up and cleaned within 24 h of sampling. The roots were first rinsed in tap water, removing all visible soil, debris and plant-roots not attached to the target plants. The roots were removed from the aboveground structure and finally rinsed in separate petri dishes in Milli-Q (MQ) H2O three times for 30 s. The root samples were then stored at −20°C until DNA extraction. Soil samples for determining soil characteristics were collected from the same hole as the plant roots. The roots were put in 50 mL Falcon tubes with 2000 μL CTAB and homogenized and extracted as described in Botnen et al. (2019). In short, a modified CTAB protocol was used (Murray and Thompson 1980; Gardes and Bruns 1993), and the DNA extracts were further cleaned using the E.Z.N.A soil kit (Omega Bio‐Tek, Norcross, GA). The internal transcribed spacer 2 (ITS2) of the rRNA gene was amplified using the forward primer fITS7a (Ihrmark et al. 2012) and the reverse primer ITS4 (White et al. 1990). Both primers were tagged with 12 base-pairs (bp) long molecular identifiers, and two bp (NN) was added to the 5’ end of both forward and reverse primers to achieve high multiplexing of samples. The PCR protocol and Illumina library preparations are described in detail in (Jacobsen et al. 2017). A total of 10 samples were sequenced twice for within-sample comparisons to ensure sequencing consistency. One library was sent to GATC Biotech for adaptor ligation and Illumina HiSeq Rapid Run 2 × 300 bp paired-end sequencing.

The V4 region of the eukaryotic 18S rRNA gene was amplified with the universal primers TAReuk454FWD1 and TAReukREV3 (Stoeck et al. 2010). The primers were tagged with a 6 bp long molecular identifier, and 1–3 bp (Ns) was added to the 5’ end of both forward and reverse primers. PCRs were performed in 26.5 μL reaction containing: 15.6 μL MQ H2O, 2.5 mL GeneAmp PCR Buffer II, 2.5 μL MgCl2, 0.2 μL dNTPs (25 nM), 0.125 μL AmpliTaq Gold DNA polymerase, 1.5 μL reverse and forward primers (10 μM), 1 μL BSA (20 mg/mL) and 1.5 μL DNA template. The PCRs were run under the following conditions: 98°C for 7 min, then 15 cycles of 95°C for 30 s, 53°C for 30 s and 72°C for 45 s, followed by 15 cycles of 98°C for 30 s, 48°C for 30 s and 72°C for 45 s, ending with 72°C for 10 min. Amplicons were checked on 1.5% agarose gel for successful amplification and cleaned using the ZR-96 DNA Clean-up Kit (Zymo Research, Irvine, CA). The DNA-concentrations were then quantified using Qubit dsDNA BR Assay, and pooled in two equimolar libraries. A mock community including four known fungi (Mycena bellarium, Pycnoporellus fulgens, Pseudoinonotus dryadeus and Serpula similis) was included in each of the libraries, and as for the ITS, nine samples were sequence twice to estimate reproducibility of the approach. The libraries were then further purified and concentrated using DNA clean and concentrator (Zymo Research) before being shipped for 2 × 250 PE sequencing on an Illumina MiSeq sequencer at Fasteris (Switzerland).

Soil samples were freeze-dried, homogenized and sieved (2 mm). Carbon (C) and nitrogen (N) contents were determined using a flash elemental analyser (Thermo Finnigan Flas EA 1112, ThermoFisher Scientific, Waltham, MA), and phosphorus (P) content was determined using a segmented flow analyser (SEAL AA3 HR AutoAnalyser, SEAL analytical Ltd, UK) at the University of Oslo. pH was measured in a 10:1 solution of deionized water and soil using LAQUA Twin pH Meter (Spectrum Technologies Inc, Aurora, IL) following the manufacturer's protocol.

Bioinformatics

Raw reads were demultiplexed using cutadapt (Martin 2011), with the following settings: anchored primers, and minimum length of 100 bp. Further quality filtering, error-correction (denoising), and merging of PE reads were conducted using DADA2 (Callahan et al. 2016) in the statistical environment R. Quality filtering was done with the following settings differing from default: maxEE = 2, truncQ = 11, minLen = 50. For the ITS2 data, a minimum of 50 bp overlap was required for merging of reads. For the longer V4 18S rRNA gene region, a minimum overlap of 15 bp was used. ITSx was used to extract the ITS2 region (Bengtsson‐Palme et al. 2013). Because of the intraspecific variability in the ITS region, the ITS2 reads were clustered at a 97% similarity level into operational taxonomic units (OTUs) using VSEARCH v2.11.1 (Rognes et al. 2016). To achieve an improved OTU delimitation, LULU (Frøslev et al. 2017), was used to combined genetically similar OTUs with similar ecological signals in the resulting OTU table. Overview of the number of reads passing the different filtering steps can be found in S2. Taxonomy of the ITS2 OTUs was assigned in UNITE v8 (UNITE Community 2019), and their functional guilds and trophic modes were assigned using FUNGuild (Nguyen et al. 2016). For the OTUs assigned as possible ECM fungi+other trophic mode (i.e. Ectomycorrhizal-Fungal Parasite-Plant Pathogen-Wood Saprotroph) in FUNGuild, the species hypotheses (SH) in UNITE was manually inspected. If the OTU was annotated to an ECM group (based on (Tedersoo, May and Smith 2010)), it was assigned as ECM, otherwise as dubious. OTUs not identified as fungi in UNITE was discarded from the dataset. For the 18S rRNA gene, we used the sequence variants from DADA2 (i.e. tentative haplotypes) as OTUs for further analyses. This is due to lower intraspecific variation in the 18S rRNA gene region. They were taxonomically assigned using a modified PR2 database (Guillou et al. 2013), were the V4 region matching the primers used had been extracted. Only 1.5% of the 18S rRNA gene reads were non-plant and 80% of those reads were from the mock fungal community, resulting in 33 677 reads for further analyses.

Phylogenetic analysis

Phylogenetic trees were constructed to verify the taxonomic assignment of the microeukaryote OTUs (i.e. non-plants) obtained from the 18S rRNA gene V4 region. All the non-plant OTUs were blasted against the non-redundant nucleotide database at NCBI and PR2 and the top five hits for each OTU were added to a reference alignment of 365 full-length 18S rRNA gene sequences covering the major branches of eukaryotic diversity. The sequences were aligned using the L-INS-i algorithm in MAFFT v7.427 (Katoh and Standley 2016) and a preliminary tree was built with FastTree2 (Price, Dehal and Arkin 2010) implemented in Geneious Prime 2019 (https://geneious.com). Based on the preliminary tree the alignment was split into separate subsets for each of the main branches of eukaryotes: Alveolata, Rhizaria, Stramenopiles, Hacrobia, Fungi, Metazoa and other Podiata (i.e. without Metazoa and Fungi). Fungi and the other podiates are hard to get good results with 18S rRNA gene on a species level, which can be related to several of the groups being polyphyletic (e.g. Amoebozoa and Apusozoa), poor taxon resolution, or potential incongruence with other common phylogenetic markers. Each of the subsets was realigned with L-INS-i in MAFFT v7.427 (Katoh and Standley 2016), and ambiguously aligned characters removed with trimAl v1.2rev59 ((Capella-Gutiérrez, Silla-Martínez and Gabaldón 2009), setting the gap threshold (-gt) to 0.3 and the minimum average similarity allowed (-st) to 0.001). Finally, maximum likelihood phylogenetic trees were constructed for each group separately with RAxML v8.0.26 (Stamatakis 2014) using the rapid bootstrapping and subsequent ML search and the GTRGAMMA model. The number of bootstraps was automatically calculated based on the MRE bootstopping criterion (Pattengale et al. 2010). Trees were visualized with FigTree 1.4.4 (Rambaut 2018).

Statistical analyses

If not otherwise specified, the statistical analyses were conducted in R (R Development Core Team 2010). For 18S rRNA gene OTUs, some of the technical replicates (i.e. samples sequenced twice) were slightly divergent, and some samples included only plant reads. Due to this, as well as the low amount of 18S rRNA gene data obtained, compositional analyses of the 18S rRNA gene data were not deemed valid. However, the included mock communities were identical in the 18S rRNA gene data, as well as close to identical in proportion of reads (0–1.3% differences). Thus, we only report overall 18S rRNA gene taxonomical composition and proportional non-plant reads for 18S rRNA gene data. In the ITS2 data, the technical replicates were similar in abundance and frequencies, and though not identical, they showed more inter-sample than intra-sample differences (GNMDS, Figure S1, Supporting Information). Additionally, paired t-tests of abundances were conducted on of each of the replicated sample pairs (raw reads), and no statistically significant differences in between repeat abundance measures of individual OTUs (P > 0.05) were found. Supported by this initial evaluation, we conducted compositional analyses of the ITS2 data, including comparisons of abundance between samples.

To ensure that the abundance data did not significantly affect statistical conclusions, both abundance and occurrence data were analysed in the form of: (1) a matrix rarefied to 1527 reads per samples transformed to presence/absence data and (2) a matrix of abundance data transformed to relative abundance (proportions) of samples with > 1500 reads. Procrustes analyses showed consistent results (P = 0.001) of the ordinations structures derived from the two different matrices.

The community structure of the root-associated fungi was analysed using global non-metric multidimensional scaling (GNMDS). Detrended correspondence analyses (DCA) were performed in parallel to ensure the reliability of the results. Correspondence in the results between the GNMDS and DCA and absence of artefacts, such as arch-effect, tongue-effect and extreme outliers, were interpreted as reliable gradients. The GNMDS settings were as recommended by Liu et al., (2008), and the DCA was run using default settings, both using the vegan package (Oksanen et al. 2012). The envfit function in vegan was used to fit vectors of the environmental variables: host plant species; host plant family; host plant order; longitude; latitude; mean July temperature and precipitation; pH and sampling locations on the ordination. To confirm the expected taxonomical structure, the potential clustering of plant families was visualized by standard error (SE) and standard deviation (SD) of their centroids using the ordiellipse and ordispider functions in vegan. To quantify the components of variation in the community composition explained by the variables mentioned above, variation partitioning, with forward selection, was performed using canonical correspondence analyses (CCA) with 999 permutations, as implemented in vegan.

To test for correlation between genetic relationships of hosts and the fungal community, distance matrices of the host DNA were constructed. The ribulose bisphosphate carboxylase large chain (rcbLA) region of the investigated host plants was downloaded from The Barcode of Life Data System (BOLD; Ratnasingham and Hebert 2007), and chloroplast trnL intron sequences were downloaded from NCBI (NCBI Resource Coordinators 2018) based on Sønstebø et al. (2010). The sequences were then aligned using the L-INS-i algorithm in MAFFT (Katoh and Standley 2013). Genetic distance matrices were constructed using the ape and adegenet package in R (Jombart 2008; Popescu, Huber and Paradis 2012). Correlation was tested by a Procrustes correlation test with 999 permutations between the distance matrix of the hosts DNA and the distance matrix of the fungal communities (Bray-Curtis distances, used in GNMDS analyses above).

RESULTS

Data characteristics, taxonomy and functional roles

A total of 847 fungal OTUs appeared in the ITS2 dataset, distributed across the 31 plant species root systems. Each root systems hosted on average 22.9 fungal ITS2-based OTUs, ranging from 6 in Eriophorum scheuchzeri to 54 in Carex rupestris (Fig. 2). Helotiales was the most abundant order, both in terms of number of reads (Fig. 3A) and OTUs (Figure S2, Supporting Information), followed by Pleosporales and Agaricales. Within the Helotiales, the majority of the ITS2 reads could not be assigned with confidence at a lower taxonomical level.

Figure 2.

Table showing number of samples, mean number of OTUs/reads and proportion of reads assigned as mycorrhizal for ITS2, mean number of ASVs and prop. of non-plant reads x1000 per host. The different colors represents different plant families. Schematic phylogenetic tree showing relationships between sampled hosts is based on: Saarela et al. 2018 (within Poaceae); Tkach et al. 2015 (within Saxifragaceae); Jordon-Thaden et al. 2010 (within Draba) and Soltis et al. 2011 (families within Angiosperms).

Figure 3.

Taxonomy (A) and functional guilds (B) on a genus level based on read abundance data of the ITS2 OTUs.

The 18S rRNA gene data was highly dominated by DNA sequences from the host plant; in all but one sample the host plant sequences made up from 96.2 to 100% of the reads (one outlier sample of Juncus biglumis had 83.9% plant DNA), which resulted in 248 non-plant OTUs. Of the non-plant reads, Fungi, Metazoa, Cercozoa and Stramenopiles made up most of the sequences (Fig. 4). Within Fungi, chytrids were the most abundant subgroup, largely missed out in the ITS2 dataset. Within Metazoa, nematodes were highly prevalent in the 18S rRNA gene data (∼12% of the non-host 18S rRNA gene OTUs), while Oomycota was the dominating stramenopile group. The groups with highest phylogenetic diversity were the cercozoan subphylum Filosa (35 OTUs), Amoebozoa (29 OTUs), nematodes (30 OTUs) followed by ciliates (19 OTUs) and oomycetes (18 OTUs) (Figures S3–S9, Supporting Information).

Figure 4.

Overview of the total taxonomy based on the 18S OTUs. Number of OTUs are indicated in the circle around the pies. The left-hand pie represents approx. Kingdom level taxonomy, while the right-hand pie approx. phylum.

In the ITS2 dataset, reads annotated as fungal saprotrophs or pathogens dominated in most plants (46.7%; Fig. 3B). The most common ITS2 OTU both in terms of reads and occurrences (266 774 reads, 87 occurrences), was affiliated with Rhexocercosporidium panacis (97.9% identity), a plant root pathogen. Likewise, a root system of Draba cf. corymbosa was dominated by one ITS2 OTU with the closest hit (96%) to the plant pathogen Olpidium brassicae, which may indicate an infection. The proportion of reads assigned as mycorrhizal fungi varied greatly, from absence in Festuca and Eriophorum scheuchzeri to 17.7% in Oxyria digyna (Fig. 3). In total, 16 ITS2 OTUs with affiliations to Glomeromycota appeared across 18 samples and eight locations in roots of the known AM hosts: Taraxacum arcticum (all samples), Ranunculus sp. (all but one sample), Potentilla puchella (all but one sample) and Micrantes nivalis (one sample), Papaver dahliana (one sample) and Trisetum spicatum (one sample). The five 18S rRNA gene OTUs with affiliation to Glomeromycota were detected in the same plant genera as for the ITS2 OTUs. Glomeromycota represented 2% of the ITS2 OTUs, and 6.6% of the fungal 18S rRNA gene OTUs. A total of 97 ITS2 OTUs were assigned as ECM fungi. The plant species hosting most potential ECM fungi, in terms of proportional read abundance, were O. digyna, Luzula confusa and Carex rupestris. There were also some ITS2 OTUs classified as orchid mycorrhizal (OM), and they all had closest hit to Serendipita sp. (86–100% identity). Only nine ITS2 OTUs were classified as DSE. The most common OTU found in Deshampsia alpina had its closest hit to Cadophora malrum (93.8%), which was classified as an endophyte.

Host identity and drivers of community structure

The composition of the fungal communities (based on ITS2), was highly structured by the host plant identity, both at species, family and order level (Fig. 5 and Figure S10, Supporting Information). In correspondence with this, we observed a correlation between phylogenetic distance among host plants, assessed both by the rbcLa (Procrustes correlation: 0.41, P = 0.001) and trnL (Procrustes correlation: 0.37, P = 0.001) cpDNA markers, and fungal community distance between the host plants. Furthermore, the variation partitioning analysis indicated that host species explained 33% of the total inertia (TI). The families Poaceae (six species) and Polygonaceae (one species), formed particularly clear clusters with no overlap in SE with other families (Fig. 5). A certain degree of clustering was also observed for Orobanchaceae (Pedicularis sp.), although some overlap with other families was found. The families Caryophyllaceae (mainly Cerastrium) and Juncaceae clustered closely together. The rest of the families cluster together on one side of the ordination, with little difference separated along the first ordination axis (Fig. 5).

Figure 5.

Global non-metric multidimensional scaling (GNMDS) ordinations based on proportional data of ITS2 operational taxonomic units in samples. Ellipses represents the standard error of the point scores of samples based on plant host families in the ordinations, lines (crosses) represents the standard deviation of the scores in the same groups. The colors represent the different plant host families. Arrows represent the direction of maximum increase for the mean July temperature, mean July precipitation, soil N, soil pH, longitude and latitude.

In addition to host plant identity, we also observed that the recovered community structure was related to certain climatic factors (July temperature and precipitation), soil pH, sampling location, longitude and latitude (Table 1). The variation partitioning indicated that sampling location explained 13.5% of TI while interaction effects between different factors accounted for a further 5.5% of TI.

Table 1.

Overview of environmental, geographical and plant taxonomy vectors fitted to the GNMDS by the envfit function in vegan, with their corresponding coefficient of determination (R2) and P-value.

Variable	R 2	P-value
Plant genus	0.4529	0.001
Plant species	0.5460	0.001
Plant family	0.3109	0.001
Plant order	0.1611	0.015
Latitude	0.1087	0.015
Longitude	0.1716	0.001
Mean July temp	0.2632	0.001
Mean July prec	0.1052	0.008
Soil pH	0.4213	0.001
Soil N%	0.0883	0.028
Sampling location	0.3035	0.001

DISCUSSION

Our study revealed a rich diversity of fungal groups associated with roots of arctic plants species expected to be non-mycorrhizal or with unclear mycorrhizal associations. The community composition of these fungal groups was strongly related to plant host identity. Although we obtained limited non-plant sequence data from the 18S rRNA gene marker, these data indicate that ITS2 data provide a highly biased view of the fungal communities, especially missing out Chytridiomycota and Mucoromycota, which were among the most abundant fungal groups in the 18S rRNA gene data.

As revealed by the ITS2 marker, different ascomycetes orders were dominant in the investigated plant roots. These results are in stark contrast to what has been observed associated with arctic ECM plants at Svalbard, where ECM basidiomycetes typically dominate (Bjorbækmo et al. 2010; Blaalid et al. 2014; Botnen et al. 2014; Davey et al. 2015; Mundra et al. 2015a, b, 2016). The observed taxonomic composition is more similar to what has been observed as aboveground endophytes (Zhang and Yao 2015) in Svalbard. We observed relatively few DSE. DSE are frequently reported from plants in alpine and polar areas, where they have been shown to play an important function in nutrient uptake (Jumpponen, Mattson and Trappe 1998; Upson, Read and Newsham 2009; Hill et al. 2019). The low proportion of DSE could partly be due to insufficiently annotated reference sequences in the databases; several of the Ascomycetes annotated with unknown functions may be DSE.

While the ascomycete orders Helotiales, Pleosporales, Chaetothyriales and Sordariales dominated in most of the 31 plant species in the ITS2 data, O. digyna represented an exception. There was a relatively high abundance of ECM basidiomycetes associated with O. digyna roots, such as Thelephorales (Tomentella), Inocybe and Cortinarius of Agaricales. Oxyria digyna belongs to the same family as Bistorta vivipara, which is a widespread ECM forming plant in the Arctic (Hesselman 1900; Gardes and Dahlberg 1996) and is well known for its extensive root structures (Maessen et al. 1983). Cripps and Eddington (2005) did not observe any signs of mycorrhiza in alpine O. digyna. ECM fungi have earlier been observed inside roots of non-ECM-forming plant species without typical ECM structures (Vrålstad, Schumacher and Taylor 2002; Vrålstad 2004; Smith and Read 2008; Schneider‐Maunoury et al. 2020). Thus, many ECM fungi may live inside the plant roots without forming ECM symbiosis. Alternatively, based on our DNA data, we may hypothesize that O. digyna forms ECM-like symbiosis in the Arctic. This must be confirmed with re-synthesis experiments and microscopy. Still, the majority of the ITS2 reads in O. digyna were classified as belonging to Pathotroph–Saprotroph fungi by FUNguild. Hence, if a mycorrhizal association is present, it is likely not the dominating root symbiosis in this species.

In the ITS2 data, we also observed 19 OTUs (26 occurrences) affiliated with Serendipita vermifera, which was functionally annotated as orchid mycorrhizal (OM) by FUNguild. However, in agreement with our results, Serendipita has been observed with a global distribution in numerous plant species not belonging to Orchidaceae (Oberwinkler et al. 2013; Botnen et al. 2014; Ray et al. 2018; Sarkar et al. 2019; Thoen et al. 2019), suggesting that different types of root associations and mycorrhiza exist in this species complex. This highlights that the automatic functional annotation through FUNguild must be critically evaluated.

We further observed AM fungi (Glomeromycota) in seven plant species. This includes some of the plant species Newsham et al. (2017) recently found to form AM in Svalbard, including the genera Taraxacum and Ranunculus. Still, a varying degree of AM has been observed in the Arctic; from no sign of AM in the plant roots studied by Väre, Vestberg and Eurola (1992) and Bledsoe, Klein and Bliss (1990) to being prevalent in other studies/sites (Dalpé and Aiken 1998; Olsson, Eriksen and Dahlberg 2004; Allen et al. 2006). Allen et al. (2006) argues that the difference observed might reflect microclimate, or a combination of biotic and abiotic factors. We did, however, detect Glomeromycota in eight of 11 locations, suggesting the AM fungi are relatively widespread in Svalbard, but likely in low abundance.

Ordination analyses of the ITS2 data clearly showed that host plant identity was the most important factor structuring the root-associated fungal communities. In this respect, our results contrast previous studies of ECM and ErM plants in the Arctic, where a low degree of host preference has been observed (Ryberg, Larsson and Molau 2009; Ryberg, Andreasen and Björk 2011; Walker et al. 2011; Timling et al. 2012; Botnen et al. 2014). The strong host preference revealed in this study may be linked to function; a large part of the OTUs was classified as putative biotrophs such as endophytes and parasites. Biotrophs might have evolved a stronger host-specificity compared to ECM and ErM fungi due to the host defence mechanisms they must encounter. Likewise, in a recent study of endophytic fungi in aboveground structures of arctic plant species, a high degree of host preference was observed (Zhang and Yao 2015). In correspondence with the ordination analyses, we observed a strong positive relationship between the phylogenetic distance of the plant hosts and the fungal community distances of the ITS2 data. Hence, closely related plants host more similar fungal communities than more distantly related plant. Hoksema et al. (2018) observed a similar pattern: They studied mycorrhizal associations in a large meta-study and concluded that the evolutionary history of the plant hosts impacts the strength of the mutualism.

Although host identity accounted for a large part of the compositional variation, the community structure was also related to climatic factors, pH and N. This was expected since previous studies have identified both climate, pH and N as important structuring factors for RAF diversity (Tedersoo et al. 2014; Timling et al. 2014; Steidinger et al. 2019).

The ITS2-primers might provide a skewed picture of the fungal diversity due to primer biases (Bellemain et al. 2010; Rosling et al. 2011; Tedersoo et al. 2015; Nilsson et al. 2019). We therefore implemented the 18S rRNA gene marker using primers assumed to amplify all eukaryotes (Hadziavdic et al. 2014). As expected, host plant sequences made up a dominant proportion of the sequences. Unfortunately, as the plant sequences typically made up 99% of the data, we obtained a rather restricted insight into the micro-eukaryotic and fungal diversity based on the 18S rRNA gene data. In further 18S rRNA gene-based studies of root-associated fungi/micro-eukaryotes, blocking-primers for plant DNA (Arenz et al. 2015) should be considered to avoid dominance of host sequences. Despite this limitation, it is interesting to see that the non-host 18S rRNA gene sequences provided a markedly different picture of the fungal diversity compared to the ITS2 data. Most strikingly, the assumed less biased 18S rRNA gene data suggest that Chytridiomycota were at least as prevalent in the arctic plant roots as Dikarya, followed by Mucoromycota. In correspondence to our 18S rRNA gene results, a relatively high abundance of Chytridiomycota has earlier been reported in arctic tundra when using 18S rRNA gene primers (Shi et al. 2015). Members of the Chytridiomycota are mainly known as decomposers, and a few are known to act as plant pathogens (Longcore and Simmons 2012). Still, this phylum probably consists of a high diversity of unknown species with unknown functions. Mucoromycota was previously thought to mainly include saprotrophic fungi (Moore, Robson and Trinci 2011). However, recent studies have shown that Mucoromycota may form mycorrhizal or mycorrhizal-like mutualism with both liverworts (Field et al. 2015; Orchard et al. 2017), and early-diverging vascular plants (Hoysted et al. 2019). We may speculate that they also can form mutualistic associations with the arctic plants studied here, which should be an interesting topic for further investigations.

We expected to find a higher diversity of Glomeromycota in the 18S rRNA gene data compared to the ITS2 data, but this expectation was not met. This might be related to the low sequencing depth: while Glomeromycota represented 2% of the ITS2 OTUs, it represented 6.6% of the fungal 18S rRNA gene OTUs. The overall mismatch between the ITS2 and 18S rRNA gene data on higher-order taxonomic composition of fungi is probably due to primer mismatches in the ITS2 primers towards the more basal fungal lineages (Bellemain et al. 2010; Rosling et al. 2011; Tedersoo et al. 2015). Although the 18S rRNA gene primers are considered eukaryote-general, we cannot exclude that mismatches also are present in these primer sites.

The 18S rRNA gene data also revealed a high diversity of other microeukaryotes and invertebrates associated with arctic plant roots. We observed e.g. a relatively high diversity of Oomycetes, which consists of several known plant pathogens (Bourke 1964; Govers, Drenth and Pieterse 1997; Thines 2014). Amoebozoa, with 29 OTUs, and Cercozoa, with 39 OTUs, were also relatively prevalent in this study, both of which contain several species commonly observed in soil, e.g. Glissomonadida (Howe et al. 2009), and Cercomonads (Flues et al. 2018) from Cercozoa; and Tubulinea in Amoebozoa (Anderson 2012, 2017). Members of Glissomonadidia have been suggested to be the dominant predators of soil bacteria (Howe et al. 2009). Since plant roots can also host a large diversity of bacteria, we may speculate that some of the Glissomonadidias observed here may feed on root bacteria. Notably, experimental studies have found that presence of mycophagous amoeba (members of Tubulina) have resulted in reduced mycorrhizal colonization of plant roots (Chakraborty, Theodorou and Bowen 1985). From Metazoa, we observed a high diversity of nematodes, which are known to parasitize plant roots (Ingham et al. 1985; Dickie et al. 2011; Kyndt et al. 2013), graze on fungi (Ingham et al. 1985), or be ‘grazed’/trapped by fungi (Bordallo et al. 2002; Yang et al. 2019). It has been experimentally shown that presence of fungal grazing nematodes affects plant host nutrient uptake (Ingham et al. 1985).

A large part of both the ITS2 and 18S rRNA gene OTUs could not be taxonomically assigned at a low taxonomic level (i.e. below order level). Poor taxonomic resolution is in general common in DNA metabarcoding studies of microbes, which is also the case here, indicating presence of many poorly studied fungi and other micro-eukaryotes associated with arctic plant roots. In the mainly desolated and environmentally stressful arctic habitats, where 80% of the plant biomass is allocated to roots (Mokany, Raison and Prokushkin 2006), plant roots represent hotspots of available C and other nutrients and much of the terrestrial diversity may therefore be centred around plant roots. The arctic plant roots might be considered as isolated ecosystems, where complex interactions between diverse communities of organisms occur. As such, arctic plant roots may represent a future avenue for research on complex biotic interactions and ecosystem functions (Iversen et al. 2015).

DATA AVAILABILITY STATEMENT

Sequence data with corresponding mapping files are available at dryad.org: https://doi.org/10.5061/dryad.6wwpzgmwr

Click here for additional data file.