Atmospheric pollution, soil nutrients and climate effects on Mucoromycota arbuscular mycorrhizal fungi.

Fine root endophyte mycorrhizal fungi in the Endogonales (Mucoromycota arbuscular mycorrhizal fungi, M-AMF) are now recognized as at least as important globally as Glomeromycota AMF (G-AMF), yet little is known about the environmental factors which influence M-AMF diversity and colonization, partly because they typically only co-colonize plants with G-AMF. Wild populations of Lycopodiella inundata predominantly form mycorrhizas with M-AMF and therefore allow focussed study of M-AMF environmental drivers. Using microscopic examination and DNA sequencing we measured M-AMF colonization and diversity over three consecutive seasons and modelled interactions between these response variables and environmental data. Significant relationships were found between M-AMF colonization and soil S, P, C:N ratio, electrical conductivity, and the previously overlooked micronutrient Mn. Estimated N deposition was negatively related to M-AMF colonization. Thirty-nine Endogonales Operational Taxonomic Units (OTUs) were identified in L. inundata roots, a greater diversity than previously recognized in this plant. Endogonales OTU richness correlated negatively with soil C:N while community composition was mostly influenced by soil P. This study provides first evidence that M-AMF have distinct ecological preferences in response to edaphic variables also related to air pollution. Future studies require site-level atmospheric pollution monitoring to guide critical load policy for mycorrhizal fungi in heathlands and grasslands.

Introduction

Mycorrhizal fungi enable plants to obtain up to 80% of their nutritional resources, mostly nitrogen (N) and phosphorus (P) otherwise bound in soil, in exchange for photosynthates (Smith and Read, 2008). Factors affecting plant nutrient availability are key drivers of ecosystem processes within heathlands, which are nutritionally poor (Read et al., 2004) and rapidly declining in their British stronghold and across Europe (Diaz et al., 2006). Mycorrhizal species richness within a habitat is directly linked to plant species richness and adaptation to the local environment (Johnson et al., 2005, 2010) and vice versa.

Some host plants demonstrate mycorrhizal fungal specificity or preference (van der Heijden et al., 2015). This is the case of the clubmoss Lycopodiella inundata, a locally rare perennial lycophyte which favours wet heathland habitats. It establishes a mutualism with multiple closely related taxa within the Endogonales (Mucoromycota) arbuscular mycorrhizal fungi (M‐AMF) clade. Many M‐AMF taxa can be harboured within a single host root and thus far, L. inundata has been found to associate primarily with M‐AMF (Hoysted et al., 2019) while other lycopods are colonized predominantly by Glomeromycota AMF (G‐AMF) (Benucci et al., 2020; Rimington et al., 2020). Other plant lineages, however, form endomycorrhizal associations with both G‐AMF and M‐AMF simultaneously (Field et al., 2016; Rimington et al., 2020). In these cases, it is difficult to distinguish microscopically the two endomycorrhizal AM fungal groups within roots, their specific functional roles or responses to environmental variables. Hyphal diameter and vesicle size help distinguish between the two fungal groups. Mucoromycota‐AMF form often‐branching thin hyphae <2 μm (typically 0.5–1.5 μm in diameter), with small (5–15 μm in length) intercalary and terminal vesicles (Hoysted et al., 2019; Kowal et al., 2020a). In contrast, G‐AMF hyphae are coarse, with a larger hyphal diameter >3 μm and longer (20–30 μm) vesicles (Orchard et al., 2017a; Hoysted et al., 2019). Both G‐AMF and M‐AMF form arbuscules in tracheophytes but due to the finer hyphal diameter of M‐AMF, previous studies have referred to these taxa as fine root endophytes or FRE (Orchard et al., 2017a; Hoysted et al., 2019).

Thus, M‐AMF‐dominated L. inundata represents an ideal system to study associations between vascular plants and this ubiquitous, but long neglected, mycorrhizal fungal clade. Moreover, it also allows investigations on how this endomycorrhizal symbiosis responds to changes in environmental factors such as atmospheric CO2 concentration (Hoysted et al., 2019), atmospheric pollution and edaphic variables. Here we focus on how the latter variables affect M‐AMF host plant root colonization, abundance and diversity which remain understudied despite the widespread distribution of M‐AMF across the land plant phylogeny, including food crops (Orchard et al., 2017a; Hoysted et al., 2018; Sinanaj et al., 2020).

Considerable M‐AMF diversity has been discovered recently and delimited into 36 species in Endogonaceae and Densosporaceae (Endogonales, Mucoromycota) (Rimington et al., 2018, 2019, 2020), but most taxonomic levels remain formally undescribed (Bonfante and Venice, 2020) and it is still unknown whether all M‐AMF are Endogonales. The placement of M‐AMF and G‐AMF as subphyla has been the subject of recent taxonomic discussion (Spatafora et al., 2016; Orchard et al., 2017b), furthered by Endogonales systematics examined by Desirò et al. (2017). However, as G‐AMF remain a distinct clade from Endogonales (M‐AMF), we maintain the higher taxonomic order of phyla herein (Tedersoo et al., 2018).

Molecular, structural and functional differentiation of M‐AMF from G‐AMF (Field et al., 2016; Field et al., 2019) have been starting points to unravel their distinct and complementary ecological roles within shared hosts and habitats (Albornoz et al., 2020). However, despite recent findings on FRE (M‐AMF) prevalence and colonization phenology in L. inundata (Kowal et al., 2020a) and the functional role of M‐AMF in L. inundata N uptake (Hoysted et al., 2019) as well as their widespread occurrence in grasses and other vascular plants (Orchard et al., 2017a; Albornoz et al., 2020), little is known about the environmental drivers of M‐AMF community composition and colonization. Understanding the ecological dynamics of M‐AMF is important to maintain diverse mycorrhizal communities supporting host plant and habitat resilience, which in the case of this rare lycopod with declining European populations, are critical.

Plants associate with different mycorrhizal fungi which facilitate nutrient mobilization and uptake of soil N and P (Read and Perez‐Moreno, 2003; de la Fuente Cantó et al., 2020), the two major growth‐limiting minerals required by autotrophic plants. For example, ericoid mycorrhizal fungi (ErM), forming the main endomycorrhizal type on N‐limited heathlands, access organic N for their Ericaceae host plants (Smith and Read, 2008; Leopold, 2016). In P‐limited habitats, G‐AMF are responsible for up to 100% of P uptake in some plant species (Smith et al., 2003) and are functionally crucial. But recent findings point to fundamentally different nutritional functions between M‐AMF and G‐AMF; G‐AMF are more efficient than M‐AMF in P uptake and transfer to liverwort hosts, independent of N availability. Conversely, host plants colonized only by M‐AMF receive substantial transfer of N, including organic N, from their fungal partners alongside P (Field et al., 2019).

Air pollution resulting in excess nutrients in the environment is one of the major threats to biodiversity (CBD, 2019; IPBES, 2019). At the ecosystem level this results in changes in plant and fungal species composition, loss and/or shifts in plant and fungal species diversity, and nutritional imbalances in plants (SAEFL, 2003; Field et al., 2014; Suz et al., 2014; van der Linde et al., 2018). The effects of pollutant critical loads (below which significant adverse effects to the ecosystem do not occur), and critical levels (above which harmful effects may occur) (CLRTAP, 2004), may be detected even once a site is no longer in exceedance, as ecosystem recovery might take time (Suz et al., 2021). Long‐term models suggest that acidic heathland habitats, which predominantly harbour ErM, are highly susceptible as their recovery can be prolonged well after pollution has ceased (Payne et al., 2013, 2017; Stevens, 2016).

Specifically, N deposition pollution affects soil characteristics which are important in shaping ErM diversity in heathlands (van Geel et al., 2020). Excess soil N or changes in pH can also influence plant nutrient availability, increasing prevalence of N‐tolerant plant species, and ultimately altering species composition through shading or competition (Stevens et al., 2018). The link between woodland‐dominant ectomycorrhizal fungi (EcM) and changes in plant nutrient status through N deposition to ecosystems is also well established (van der Linde et al., 2018; Suz et al., 2021), and similar effects have been reported for N additions to G‐AMF (Corkidi et al., 2002; Johnson, 2010; Liu et al., 2012; Jiang et al., 2018; Ceulemans et al., 2019). In this study we investigated abiotic soil and environmental interactions which influence M‐AMF plant root colonization, richness and community composition in heathlands across environmental gradients while generating diversity data for these groups of fungi. We hypothesized that N deposition will be one of the main factors influencing these fungi. We also studied whether modelled air pollution could be related to soil covariates known to affect G‐AMF (Johnson, 2010). Finally, we investigated the presence of G‐AMF in L. inundata roots and explored their potential contribution to host plant nutrition.

Experimental procedures

Study sites, collection periods and sampling of plant root, plant tissue and soil

We studied 12 heathland sites covering climate and air pollution gradients within Britain and northern Europe, seven from southern Britain, three from northern Scotland and two from the Netherlands (Fig. 1). We sampled soil and plants from three 1 m2 subplots per site. At sites where the distribution of the L. inundata population did not allow three 1 m2 subplots, we collected soil and plants from three population clusters as far apart as possible.

Fig. 1

An example of the European Monitoring and Evaluation Programme (EMEP) pollutant data used for statistical analyses. Nitrogen dioxide (NO2) concentrations are shown on a gridded map at a resolution of 0.1° × 0.1° (January 2019). Darker colours represent low concentrations (minimum 0.01 μg m−3). Lighter colours represent high concentrations (maximum 10.43 μg m−3). Lycopodiella inundata sites are represented by red dots. Golden circles around sites show the four concentration gradients, with Scotland having the lowest values followed by Cornwall, and a west to east gradient towards the Netherlands. Source: EMEP (.

Soil and plant data were paired by subplot over the three collection periods T0, T1 and T2, each 5–7 weeks long. For each of the three collection periods, the new pairs of subplots were randomly selected across a grid covering the L. inundata population. If no L. inundata was present, the next nearest subplot was sampled. Attention was placed on minimizing disturbance within the plot. Due to site access restrictions during the Covid pandemic in 2020, three of the 12 sites were not sampled at T2. For all but one subplot (see below) we collected 6–10 healthy‐looking L. inundata plants by selecting bright green stems without signs of decay.

At each subplot we removed at least 500 g of soil (dry weight) from the organic matter (OM) layer avoiding the mineral layer, up to 10 cm deep. The OM layer typically ranged from 0.5 cm to 5 cm, sometimes within a single subplot. If the OM layer of a single subplot was insufficient, the three subplot samples were aggregated. Plant litter was removed before placing each sample into a sealed bag.

Soil was analysed for the following edaphic variables: pH, electrical conductivity [EC1:5 (μS cm−1)], N (% wt./wt.) and C (% wt./wt.) content, C:N, macronutrients [P (mg L−1), K (mg L−1), Mg (mg L−1), Na (mg L−1), S (mg L−1), Ca (mg L−1)] and the micronutrient Mn (mg L−1). We initially used Olsen P index to measure plant available P but added the Mehlich III technique in T2, as this is more reliable for acidic soils (Mehlich, 1984; Wolf and Baker, 1985). Soil extraction methods and chemical analyses are detailed in Methods S1.

In addition, we pooled three soil core samples from 1 m2 subplots located outside of the Lycopodiella zones (named as ‘outlier’ soil samples). In these areas the vegetation composition had unequivocally shifted to grasses or sedges lacking Erica tetralix, Calluna vulgaris and/or Drosera spp.

The roots collected from T0 and T1 were processed for quantification of FRE (M‐AMF) colonization in Kowal et al. (2020a, 2020b) and roots collected during T2 were measured and prepared for microscopy analyses following the same procedures. Within four days from collection, roots were placed in either 70% (vol./vol.) ethanol or cetyltrimethylammonium bromide (CTAB) extraction buffer and stored at 4°C (ethanol) or −20°C (CTAB).

We also measured L. inundata chlorophyll fluorescence, leaf nutrient content and vegetation composition in each site, as detailed in Methods S2. Plant families within the same subplots were linked to mycorrhizal type (colour coded in Table S1) following Brundrett and Tedersoo (2018).

Root colonization of endomycorrhizal fungi

We examined M‐AMF presence and/or absence in T2 root samples at individual root, plot and site levels (T0 and T1 roots were analysed in Kowal et al., 2020a). We further quantified abundance of colonization for T0, T1 and T2 roots by estimating cell coverage by M‐AMF hyphae or vesicles within each individual root using four categories (rare: 3–5 cells; low: ≤10% of cells; medium: 11%–25%; and high: >25%). We also examined presence and/or absence of G‐AMF in each root. Microscopically, M‐AMF and G‐AMF were differentiated by hyphal diameter. When arbuscules were present, they were attributed to either M‐AMF or G‐AMF only when visibly linked to the distinguishing hyphal morphotype as seen in Fig. S1 with G‐AMF.

Molecular identification, community composition and phylogenetic analyses of mycorrhizal fungi

A subsample of roots from nine plants per site was randomly selected (three plants per plot, where plots were designated) at T1 except at Aldershot, where only five plants were collected due to plant conservation measures in place. At T2, we sampled roots from two of the eight sites. For each plant, one half‐root measuring at least 0.6 cm was used for DNA extraction. If no roots >0.6 cm were available, two or three smaller roots were pooled in the same DNA extraction. Total DNA was extracted with the DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions with modifications, as detailed in Methods S3.

Amplification of the fungal 18S rDNA region, PCR product cloning, re‐amplification of Escherichia coli colonies and DNA sequencing were carried out following Rimington et al. (2019) with 8–12 E. coli colonies re‐amplified per sample. Sequences (ca. 600 bp) were assigned to subphyla using NCBI BLASTn (Camacho et al., 2009) and those matching Endogonales (or Mucoromycota) and Glomeromycota were further sequenced using the primers NS3 and NS5 (White et al., 1990). Sequences were edited, assembled into contigs (ca. 1500–1700 bp) and aligned using MAFFT v 7.017 (Katoh et al., 2002) on Geneious v. 8.1.9 (Kearse et al., 2012). The alignment was edited manually and the low‐quality ends were trimmed. The algorithm UCHIME2 (Edgar, 2016) within USEARCH v. 11.0.667 (Edgar, 2010) was used to identify chimaeric sequences. The full‐length sequences (trimmed to 1513 bp) belonging to Endogonales (or Mucoromycota) and Glomeromycota were clustered into Operational Taxonomic Units (OTUs) at a 98% similarity threshold (Rimington et al., 2019) using USEARCH (Edgar, 2010).

A representative DNA sequence of each Endogonales OTU was queried against the SILVA SSU database (Quast et al., 2013) using SINA v. 1.2.11 with EMBL‐EBI/ENA taxonomy as a reference, and renamed according to the best scoring hit. The representative sequence of each Glomeromycota OTU was queried against the AMF fungal DNA database MaarjAM (Opik et al., 2010) and named according to the best scoring hit. Representative sequences of each OTU, defined by the UCLUST algorithm, were aligned with Endogonales DNA sequences from Bidartondo et al. (2011), Desirò et al. (2013), Field et al. (2016), Hoysted et al. (2019), Rimington et al. (2019) and Endogonales sequences from GenBank, using MAFFT v 7.017 (Katoh et al., 2002) on Geneious v. 8.1.9 (Kearse et al., 2012) and used for phylogenetic analysis. A Maximum Likelihood (ML) tree was built using RAxML‐HPC on BlackBox on CIPRES Scientific Gateway (Miller et al., 2010) using 1000 bootstrap iterations. We chose two species as outgroups (Basidiobolus ranarum and Olpidium brassicae; GenBank accessions AY635841 and DQ322624 respectively). Taxonomic notations of Densosporaceae and Endogonaceae were assigned by comparison with Rimington et al. (2019). Bayesian inference was performed using MrBayes v.3.2.7 (Huelsenbeck and Ronquist, 2001) on XSEDE on CIPRES Scientific Gateaway with 1.000.000 generations and both trees were visualized using FigTree v.1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Accumulation curves of the OTUs belonging to M‐AMF were constructed in R using the function specaccum (Oksanen et al., 2020). A Bray–Curtis dissimilarity matrix was generated based on the relative abundance of each OTU in each site and we used the Hellinger transformation (Legendre and Gallagher, 2001). Non‐metric multidimensional scaling (NMDS) ordination was used to visualize the main factors affecting M‐AMF community dissimilarities among sites and 12 environmental variables were fitted using the ‘envfit’ function in the vegan package (Oksanen et al., 2020).

Air pollution and climate data

Pollutant concentrations and depositions were compiled for each of the 12 L. inundata sites, based on modelled data from the Norwegian Meteorological Institute (MET) (EMEP, 2021). For each site concentrations in the air [μg (S or N) m−3] of SO2, NO2, NH3, NO (NO2 + NO), NH (NH3 + NH4 +) and OxN (oxidized reactive N), and deposition [mg (S or N) m−2) (dry and wet) of oxidized S, oxidized N and reduced N were extracted from the dataset (Table S2). Data were calculated on a grid, at a resolution of 0.1° × 0.1°. Further details are described in Methods S4. We extracted mean monthly precipitation and air temperature data for 2018, 2019 and 2020 from www.worldweatheronline.com.

Statistical analyses

All regression models and correlations were carried out in R v4.0.2 (R Core Team, 2020). We fitted in the models three root fungal colonization response variables measured at T0, T1 and T2: % M‐AMF colonization; % G‐AMF colonization; and M‐AMF colonization density per root scored using the rare/low/medium/high categories detailed above. Linear regression models were fitted to all responses with either a binomial (for % colonization responses) or normal distribution assumed (for colonization density).

Explanatory variables were selected from the soil, climate, and atmospheric pollution data, centred, and scaled prior to analysis. Our initial data exploration yielded two main models, bulk density (BD) model and the All‐sites model. Data were analysed at both subplot level (n = 48, 58) and at site level (n = 24, 28) respectively, with averages taken across subplots. Ordinal colonization intensity data were converted to a numeric scale by assigning each root fragment the class midpoint, averaged at plot or site level, and analysed after square root transformation. Additional assumptions used for model development as well as limitations are detailed in Methods S5.

We also tested whether presence of G‐AMF correlated with M‐AMF presence using Spearman's rank‐order test and whether G‐AMF structures were more likely to be alone or present with M‐AMF using a Chi‐square test.

Fungal community composition and OTU richness were assessed at T1 only; therefore, it was not possible to fit a multiple regression with all predictors due to the reduced number of observations. Instead, Spearman's correlations were calculated between fungal OTU richness and each of the environmental covariates at site and plot level. Pairwise correlations were also calculated between M‐AMF colonization and soil, climate, and air pollution variables, chlorophyll fluorescence and leaf chemistry. Correlation coefficients of r > 0.70 were considered important.

Results

Root endomycorrhizal fungal colonization and community assessments

Microscopic assessments of root fungal colonization

Quantification of M‐AMF and G‐AMF structures observed microscopically in L. inundata roots for each of the sampling points of spring 2019 (T0), autumn 2019 (T1) and spring/summer 2020 (T2) are summarized in Fig. 2 and Table S3 (T0 and T1 were reported in Kowal et al., 2020a). See Table S4 for site names, abbreviations, locations and sampling dates. Colonization of L. inundata roots by M‐AMF significantly differed across seasons T0 and T1, but in T2 the percentages of colonization were intermediate between T0 and T1. There were also significant differences between individual root M‐AMF colonization densities across sampling seasons with significantly higher intensity observed in T1 and T2 than in T0 (Table S5).

Fig. 2

Mucoromycota arbuscular mycorrhizal fungi (M‐AMF) and Glomeromycota arbuscular mycorrhizal fungi (G‐AMF) root colonization (logit transformed) across sites, by sampling point. T0 = spring 2019, T1 = autumn 2019, T2 = late spring/summer 2020. Site NLB was not sampled in T0 and AL, MH and SH were not sampled in T2. Data from T0 and T1 are published in Kowal et al. (2020a).

We found more G‐AMF hyphae and vesicles in T2 roots, in approximately 15% of the roots overall and in up to 25% of the roots at one site (TH) (Table S3). Arbuscule‐like structures were rarely observed (see Fig. S1 for an example of these structures). Colonization of M‐AMF and G‐AMF were only weakly correlated at plot (rho = 0.36) and site (rho = 0.30) levels. However, G‐AMF were more likely to be present when roots were also colonized by M‐AMF (x 2 = 16.07, df = 1, p < 0.001).

Molecular identification of root fungal communities

Across the 12 sites, a total of 190 18S rDNA sequences (182 M‐AMF and 8 G‐AMF) were obtained (NCBI accessions OM214587–OM214776, Table S6a) resulting in 46 OTUs (39 M‐AMF and 7 G‐AMF), of which 30 were singletons and 16 included more than two sequences (Table S6b).

The 39 M‐AMF taxa were detected in all sites across 92.5% of the samples, whereas G‐AMF were only found in five sites (CE, HB, MP, ST, TH) in 7.5% of the samples. We observed G‐AMF structures microscopically, but we did not detect DNA of G‐AMF structures in roots of L. inundata in the Netherlands samples.

The number of OTUs per site ranged from 4 (HB) to 13 (BD). Three M‐AMF OTUs (OTU1, OTU2 and OTU3) accounted for 51.6% of the total sequences. The OTU1 was present at all but one site (AL) and was the most abundant OTU at eight sites including all regions, i.e. Scotland, Netherlands, Cornwall and other southeast England sites. At the remaining sites where OTU1 did not dominate, OTU2 or OTU3 were the most abundant. There were three OTUs unique to Scotland and one OTU unique to the Netherlands. No clear geographical pattern was identified, with sequences from England, Scotland and Netherlands evenly distributed across the phylogenetic tree (Fig. 3). The OTU accumulation curves varied across sites, with some sites closer to reach an asymptote (e.g. HB, NLA, Fig. S2a). The NMDS ordination showed that soil P influences M‐AMF community dissimilarities across sites (Fig. S2b).

Fig. 3

Maximum likelihood phylogeny of Endogonales (Mucoromycota) and Glomeromycota OTUs colonizing Lycopodiella inundata roots (this study) and including previously published 18S M‐AMF sequences. Support values derived from 1,000 bootstrap iterations are shown if greater than 50. Each OTU is represented by one centroid sequence as defined by the UCLUST algorithm and the number of sequences assigned to each OTU is shown in brackets. Singleton OTUs are highlighted in red font. Coloured squares represent the geographical origin of sequences, blue = Scotland, red = England and orange = the Netherlands.

Maximum likelihood (ML) (Fig. 3) and Bayesian (Fig. S1) analyses of the 18S rDNA sequences produced similar tree topologies at the main branches. Glomeromycota‐AMF and M‐AMF OTUs clustered in two distinct clades and 13 distinct clades were identified in the Endogonales (Fig. 3). Endogonales M‐AMF OTUs were distributed across seven clades within Densosporaceae (three clades) and Endogonaceae families (four clades), with a larger number of sequences clustering in the Densosporaceae (clades I, II and IV; Fig. 3). The presence of G‐AMF OTUs was negatively correlated with M‐AMF OTU richness (r = −0.77).

Vegetation surveys and soil characteristics of Lycopodiella inundata subplots

There were no significant differences in dominant plant species composition across sites and all 30 subplots surveyed per season contained ErM, G‐AMF and non‐mycorrhizal (NM) associated plants (Table S1). One to three subplots contained EcM‐associated seedlings (e.g. pine and birch). Calluna vulgaris and/or Erica tetralix dominated the ErM shrub community cover and Molinia caerulea the G‐AMF‐associated grass community. Drosera spp. (NM) along with Juncaceae and Cyperaceae species (NM/G‐AMF) were present at all sites (although not in all the subplots surveyed) and Sphagnum spp. (NM) were present in half of the sites. Lycopodiella inundata population cover varied across sites, subplots and seasons but formed denser carpets at PL (T1) than elsewhere. Overall mean strobili presence was greater in T1 than T2, 32% (in 30 subplots) compared with 3% (in 25 subplots) respectively. We did not find any correlation between M‐AMF colonization and strobili counts across both sampling times.

Soil chemical characteristics for the three sampling periods are summarized in Table S7. There were many strong correlations (r > 0.70) between soil and air pollution variables, as well as temperature and precipitation at the plot and site levels (Table S8a and b).

Model relationships to environmental covariates

Mucoromycota‐AMF colonization

At plot level (Table 1; Fig. 4), several covariates were significant in both models. Soil C:N, S and P were negatively related with M‐AMF colonization (p < 0.05, <0.05, 0.01 respectively), while soil electrical conductivity and Mn were positively related with M‐AMF colonization (p < 0.05, 0.001 respectively).

Table 1

Plot‐ and site‐level results of All‐sites and bulk density (BD) models with percent of roots colonized by Mucoromycota arbuscular mycorrhizal fungi (M‐AMF) and Glomeromycota AMF (G‐AMF) as response variables using quarterly European Monitoring and Evaluation Programme data.

	M‐AMF				G‐AMF
	Site level		Plot level		Site level		Plot level
	All‐sites model	BD model	All‐sites model	BD model	All‐sites model	BD model	All‐sites model	BD model
Soil characteristics
pH	−0.177 (0.17)	−0.91 (0.23)***	0.21 (0.15)	0.08 (0.18)	0.368 (0.47)	1.61 (0.69)*	0.41 (0.37)	0.73 (0.45)
Conductivity	−0.045 (0.22)	0.06 (0.24)	0.52 (0.23)*	0.51 (0.25)*	−1.051 (0.54)*	−2.17 (1.45)	−1.08 (0.59)	−1.15 (0.65)
C:N	−0.147 (0.09)	−0.46 (0.11)***	−0.44 (0.18)*	−0.52 (0.20)*	0.051 (0.23)	−0.05 (0.31)	−0.71 (0.34)*	−0.53 (0.36)
Total N	−0.158 (0.17)	0.10 (0.23)	−0.23 (0.13)	−0.10 (0.21)	0.355 (0.39)	0.09 (0.89)	0.04 (0.38)	−0.2 (0.5)
Mg	0.186 (0.21)	−0.45 (0.24)	n/a	n/a	1.732 (0.64)**	n/a	n/a	n/a
Ca	0.104 (0.21)	0.59 (0.20)***	0.21 (0.20)	0.32 (0.18)	−3.21 (1.88)	−0.24 (1.01)	−0.15 (1.79)	−0.71 (2.38)
Mn	0.711 (0.12)***	1.00 (0.14)***	0.34 (0.11)***	0.39 (0.11)***	−0.104 (1.86)	−0.37 (1.01)	−0.53 (0.90)	−0.32 (0.78)
S	−0.33 (0.12)**	−0.76 (0.21)***	−0.39 (0.16)*	−0.54 (0.19)***	0.171 (0.25)	1.29 (0.74)	0.24 (0.23)	0.52 (0.34)
P	−0.214 (0.08)**	0.00 (0.11)	−0.31 (0.10)**	−0.24 (0.11)*	0.504 (0.26)	−0.1 (0.31)	0.07 (0.19)	0.08 (0.21)
Bulk density	Not included	0.63 (0.20)***	Not included	0.32 (0.18)	Not included	−1.16 (0.52)*	Not included	−0.55 (0.41)
Climate
Precipitation	0.145 (0.14)	0.46 (0.28)	0.12 (0.11)	0.24 (0.19)	−1.207 (0.38)**	−1.68 (0.67)*	−1.09 (0.34)**	−1.35 (0.45)***
Atmospheric pollution
NO x concentration	0.004 (0.18)	−0.45 (0.26)	0.22 (0.21)	0.15 (0.28)	0.997 (0.44)*	0.92 (0.89)	0.13 (0.36)	0.35 (0.52)
Total N deposition	0.282 (0.16)	−0.70 (0.19)***	−0.13 (0.18)	−0.40 (0.22)	−0.054 (0.36)	0.5 (0.46)	−0.48 (0.33)	−0.13 (0.41)
Sampling period
Time T1	4.728 (0.56)***	4.82 (0.60)***	4.60 (0.47)***	4.88 (0.55)***	9.964 (4.23)**	1.83 (2.63)	2.55 (2.12)	2.79 (2.79)
Time T2	4.12 (0.47)***	3.81 (0.46)***	4.84 (0.47)***	4.84 (0.58)***	8.219 (4.11)*	1.01 (2.28)	1.74 (2.21)	1.64 (2.85)

Values are coefficient estimates with standard errors in brackets. Asterisks denote significance at p < 0.05 (*), p < 0.01 (**) and p < 0.001 (***). Significant relationships found in both models coloured green, and significance in only one of the two models, orange.

Fig. 4

Intersection of significant relationships between Mucoromycota arbuscular mycorrhizal fungi (M‐AMF) colonization of Lycopodiella inundata roots and explanatory covariates tested in both linear regression models (All‐sites and Bulk Density). (Aboveground) N = Nitrogen deposition (image: L. inundata host plants in the wild; bar = 1 cm); (Belowground) Soil chemistry covariates (image: L. inundata ink‐stained roots showing typical M‐AMF colonization, bar = 20 μm). Diagonal arrows indicate a positive or negative relationship with M‐AMF colonization. Green‐bound circles indicate strong relationships occurring in both linear regression models at the site or plot levels; green‐filled circles indicate strong relationships observed in both models and both levels. Relationships occurring at only one level, or one model, are circled orange. Explanatory coefficient estimates and p‐values are detailed in Table 1. Atmospheric pollution data are based on the last quarterly averages prior to harvesting the roots (Table S2).

At site level (Table 1; Fig. 4) the All‐sites model (which excluded bulk density, see Statistical analyses) showed that M‐AMF colonization was strongly related to several soil parameters. Soil S and P were negatively related (p < 0.01) and Mn positively related (p < 0.001) to M‐AMF colonization. In contrast, the model including the soil covariate bulk density (‘BD model’) showed a strong negative relationship between M‐AMF colonization and total N deposition, soil pH, C:N ratio and S (p < 0.001); and a strong positive relationship to Ca, Mn and bulk density (all p < 0.001). Despite strong evidence of K relationships with plot level M‐AMF colonization, as it was highly collinear at site level with soil Mg and Mn, it was not included in the models (Methods S5).

Both site and plot level models indicate colonization density (rare–low–medium–high, per individual root) was negatively related to soil S (p < 0.01, 0.05) respectively. At the plot level, this occurred with both models (Table S9). Magnesium was also negatively related at the site level only (p < 0.05).

Sampling season (as categorical variable) was a significant covariate (p < 0.001) in both models; season‐specific analyses may be seen in Table S10. At T1, M‐AMF colonization was correlated to temperature (r = −0.70). We also found a significant negative correlation between M‐AMF colonization and soil S (r = −0.78), P (r = −0.71) and conductivity (r = −0.90).

Mucoromycota‐AMF OTU richness

At T1, we did not find a strong correlation between air pollution covariates and M‐AMF OTU richness (Tables S8a and b), except for a weak negative correlation with NHx concentration and N deposition at the site level. However, soil C:N was negatively correlated with M‐AMF OTU richness (r = −0.75). Root M‐AMF colonization proportions and density per root were strongly correlated (r = 0.92) but both covariates were weakly correlated with M‐AMF OTU richness.

Glomeromycota‐AMF colonization

Response signals were generally less strong regarding G‐AMF colonization (Table 1, Table S8a and b). There was some indication of a positive relationship with NO concentration (p < 0.05) at the site level. Significant negative relationships were measured with soil conductivity and C:N (p < 0.05) while significant positive relationships were measured with pH (p < 0.05) and Mg (p < 0.01). Glomeromycota‐AMF presence was strongly negatively related to mean monthly precipitation in both models at plot and site level.

Plant nutrient tests: leaf chemistry and plant chlorophyll fluorescence

Across sites, mean total leaf C was significantly lower in T0 than in T1 while we did not find differences in mean total leaf N (Table S11). Percentage of roots colonized by M‐AMF was significantly correlated with leaf C content at the site level (r = 0.75) (Table S8b). We found no significant correlations between field measurements of L. inundata chlorophyll fluorescence (Fv:Fm) and its root colonization by M‐AMF and G‐AMF (Fig. S4).

Discussion

Lycopodiella inundata consistently hosts M‐AMF

Phylogenies and OTU assemblages using near‐complete 18S DNA sequences indicate that L. inundata hosts at least 39 M‐AMF OTUs distributed within seven taxonomic clades across 13 Endogonales clades. This represents a significant advance in our knowledge of M‐AMF diversity among early‐diverging vascular plants. Earlier phylogenies including M‐AMF detected in roots of L. inundata were either poorly resolved within Densosporaceae (Rimington et al., 2015) or based on partial (400–700 bp) 18S rDNA sequences (Hoysted et al., 2019). Furthermore, our OTU accumulation curves indicate that there may be yet more diversity to uncover (Fig. S2a). The presence of three main M‐AMF OTUs across most sites highlights the ubiquitous presence of these fungal lineages and their consistent association with L. inundata. We only detected G‐AMF in 7.5% of samples across the 12 sites. These T1 molecular findings agree with the root colonization observed across sites, except for the Netherlands, where G‐AMF structures were observed microscopically but not detected in the DNA analyses.

Environmental factors affect M‐AMF colonization, community richness and composition

Our analyses show that atmospheric pollution may be affecting M‐AMF colonization levels and OTU richness, at least indirectly. This is indicated by the negative relationship between percentage of roots colonized and total N deposition. Also, our model indicates that soil C:N is a key variable affecting M‐AMF colonization negatively. This is similar to other studies focused on G‐AMF (Johnson et al., 2010; Tedersoo and Bahram, 2019). In fact, soil C:N has been shown to be affected negatively by N deposition across Europe (Mulder et al., 2015), in both grasslands and moorlands (Evans et al., 2006; Volk et al., 2016).

We found indirect links between air pollution variables and M‐AMF OTU richness. Soil C:N was negatively correlated with M‐AMF OTU richness, raising the possibility that higher N atmospheric pollution could indirectly lead to a decrease in M‐AMF symbionts and changes in diversity and composition of AMF functional groups as previously observed in soil carbon richness gradients (Johnson et al., 2013). In forest systems, atmospheric N negatively affects diversity and composition of mycorrhizal fungi (Lilleskov et al., 2019). However, M‐AMF's specialization for providing N from organic sources to host plants, at least in microcosm experiments with dual G‐AMF/M‐AMF host liverworts (Field et al., 2019), may be critical for host plant resilience when changes in levels of inorganic N from air pollution results in an imbalance of N resources. Nonetheless, given that little is still known with respect to C:N dynamics in G‐AMF (Corrêa et al., 2015) and M‐AMF, field manipulation studies are needed to improve predictions of these feedback cycles. Thus far, the direct impact of N deposition on M‐AMF diversity and the role of certain M‐AMF OTUs in mining organic N in heathlands remain to be tested. The unique M‐AMF OTUs identified from Scottish root samples, where atmospheric N is lowest, may provide a clue where to begin investigations.

The strong relationships we found between several soil variables and the extent of M‐AMF colonization are consistent with similar research on Trifolium subterraneum colonized by M‐AMF (Albornoz et al., 2020). We also found distinct soil environmental niches between M‐AMF and G‐AMF presence regarding soil pH. Our BD model (site level) shows a negative relationship between pH and M‐AMF colonization, but a positive relationship with G‐AMF. This agrees with Tedersoo et al. (2020) where they found contrasting ecological preferences between M‐AMF and G‐AMF in extensive soil environmental DNA sampling across habitats in Estonia and North Latvia, with M‐AMF preferring acidic soils. In that study, pH had the strongest effect on the diversity of fungi. This finding is also consistent with Albornoz et al. (2022) which found M‐AMF preference for acidic soils across a wide sampling of agricultural sites in Australia. Still, more studies have focused on G‐AMF ecological requirements, without differentiating M‐AMF. For instance, in semi‐natural plant–soil feedback systems, soil pH is the principal driver affecting G‐AMF community composition (Dumbrell et al., 2010). Similarly, presence of keystone G‐AMF taxa in agroecosystems is best explained by soil pH, P levels, bulk density and salinity (Liu et al., 2014; Banerjee et al., 2019), but these studies ignored M‐AMF. Our NMDS analysis shows a similar influence of soil P in M‐AMF composition across sites.

Both models indicated that M‐AMF colonization is related to soil S (negatively) and Mn (positively) at both plot and site levels and M‐AMF density per root is also negatively related to both S and Mg. Soil sulfate is linked directly to the atmospheric concentrations of sulfur dioxide (Feinberg et al., 2021). Reductions in SO2 levels over the last decades in both the United Kingdom and the Netherlands may be relevant to the observed pH effect on colonization in this study. In a related heathland manipulation study, Tibbett et al. (2019) found that elemental S additions were the primary factor affecting soil pH and a negative G‐AMF colonization response. Our temporal data suggest significant negative correlations between SO deposition and M‐AMF colonization in spring, followed by negative correlations between soil S and M‐AMF colonization in autumn. Thus, a reduced SO deposition followed by limited availability of soil S could be affecting root colonization by M‐AMF and heathland recovery in general, but further work is needed to assess this M‐AMF specific response. It is possible that in acidic soils, such as in heathlands, high levels of these micronutrients are required because of their poor solubility (Millaleo et al., 2010). The precise role of micronutrients such as Mn and Mg has been less tested than other soil variables and nutrients but our study indicates they are important indirect factors likely affecting soil pH, at least in heathlands.

Presence of Glomeromycota in Lycopodiella inundata roots

Four of the eight G‐AMF DNA sequences were from roots also colonized by M‐AMF, confirming that some plants of L. inundata are co‐colonized by both groups of fungi. This finding is in keeping with the association between G‐AMF and other Lycopodiopsida (Rimington et al., 2015) and dual colonization by M‐AMF and G‐AMF across different plant lineages as seen with mutualisms in liverworts (Field et al., 2016; Rimington et al., 2020), grasses (Hoysted et al., 2019) and angiosperms (Orchard et al., 2017a). Nonetheless, the strong preference for M‐AMF by L. inundata is certainly consistent.

We found a positive relationship between soil P and the presence of G‐AMF in some L. inundata roots. Experimental microcosms using dual G‐AMF and M‐AMF host plants show that G‐AMF may be more efficient than M‐AMF in supporting plant P acquisition (Field et al., 2019; Hoysted et al., 2019). While not directly comparable, this could help explain the rare presence of G‐AMF. Furthermore, we observed that the rare G‐AMF structures within L. inundata were more likely to be present in roots also colonized by M‐AMF rather than occurring on their own. This suggests that L. inundata largely relies on M‐AMF for its P requirements as previously assumed (Hoysted et al., 2019), but may also recruit G‐AMF symbionts under certain opportunistic conditions, such as when N deposition is high and/or M‐AMF OTU richness is lower. This differs from previous studies (Orchard et al., 2017a) which suggested that M‐AMF enhance host plant P uptake rather than providing primary access to P. Further field studies focusing on the functional role of these fungi are needed.

Our model also showed a positive relationship between NO and G‐AMF. However, despite data suggesting some Glomus (G‐AMF) species might be N‐tolerant, responses to N deposition by G‐AMF can be variable (Treseder et al., 2007, 2018) and given their rarity in L. inundata roots, extensive root DNA sequencing across pollution gradients would be required to confirm this relationship.

Limitations of local and national modelling interactions and grid resolutions

The models may be underestimating or masking the relationship strength between M‐AMF colonization and N deposition due to limitations relating to EMEP model resolution and collinearity among the N covariables (Methods S4 and S5). Work carried out as part of the UK Joint Nature Conservation Committee‐led project Nitrogen Futures has demonstrated the variability within grids when modelled at different resolutions. Mean NH3 concentration and N deposition across all locations in a grid was higher with a resolution of 1 × 1 km2 compared to 2 × 2 m2. Conversely, the maximum for any location in the grid was higher at the lower resolution (Thomas et al., 2020).

Future directions

The season‐specific correlations observed provide some evidence that relationships with covariates may not be consistent over time due to variation in climate factors and changes during the growing season. Therefore, future work should investigate potential interactions between time of year and abiotic drivers.

Leaf C content significantly correlated with M‐AMF colonization, suggesting there may be a link between root colonization and plant tissue C content, as previously shown (Zhu et al., 2014; Mathur et al., 2018), and this may provide a simple non‐invasive tool to infer relative host C allocation to these fungi.

We observed that one site in southern England with low M‐AMF colonization (Aldershot), which is also the L. inundata population under greatest decline, lacked the dominant OTU1, present however in all other study sites. Roots in neighbouring sites where OTU1 was present had much higher colonization despite having similar N deposition values. This could suggest that OTU1 may provide their host plants with coping mechanisms to N deposition stress. It is also possible that the widely dispersed OTU1 is present but not yet detected in Aldershot by our sampling effort or it may be unable to compete with the vegetation changes occurring near this population. Transfer experiments of plants hosting OTU1 from nearby thriving populations in southern England and monitoring whether this facilitates population stability and growth over time would allow testing this hypothesis.

We expected to find G‐AMF more commonly and opportunistically colonizing L. inundata where G‐AMF plants (e.g. Molinia caerulea) were more dominant and less where ErM, NM and/or EcM plants were more frequent. However, given the low variability in vegetation composition across our sites, we were unable to test vegetation as a categorical predictor of M‐AMF and G‐AMF colonization of L. inundata roots. Molecular analyses coupled with experimental microcosms – testing donor and target plants – would help disentangle this putative association.

This study found that several soil characteristics influence M‐AMF colonization and richness. Furthermore, our analyses also indicate that atmospheric pollution may indirectly interact with these same soil characteristics, and therefore indirectly influence M‐AMF colonization and community composition. However, higher‐resolution air pollution monitoring is needed, at the field experiment scale, to couple air pollution monitoring data with M‐AMF resilience and diversity measures. Without such an investment to test and set air pollution critical load and levels specifically for mycorrhizal fungi in vulnerable habitats such as heathlands, we may be overlooking irreversible ecosystem changes occurring belowground.

Appendix S1. Supporting Information.

Click here for additional data file.

Fig. S1. Images of an ink‐stained root showing arbuscule‐like structures associated with Glomeromycota arbuscular mycorrhizal fungi (bar = 80mm).

Fig. S2. Endogonales (Mucoromycota) arbuscular mycorrhizal fungi (M‐AMF) community analyses. (a) M‐AMF Operational Taxonomic Unit (OTU) accumulation curves per site. b) M‐AMF community composition of the 12 sites displayed using NMDS.

Fig. S3. Bayesian phylogenetic inference (1,000,000 generations) of Endogonales (Mucoromycota) (M‐AMF) and Glomeromycota OTUs colonising Lycopodiella inundata roots (this study) and previously published 18S rDNA M‐AMF sequences.

Fig. S4. Plant tissue fluorescence (Fv:Fm) by site and sampling point comparing Lycopodiella inundata (L) and Molinia caerulea (M), sampling periods T1 and T2.

Table S1. Compilation of site vegetation surveys by plot (T1 and T2) including plant mycorrhizal type.

Click here for additional data file.

Table S2. Pollution data (source: EMEP 2021).

Table S3. Summary of root colonisation across sites, by sampling season.

Table S4. Site names, sampling dates and number of model observations.

Table S5. Mucoromycota arbuscular mycorrhizal fungi (M‐AMF) colonisation densities by site.

Table S6a. Summary of full‐length 18S rDNA sequences obtained, distributed among the 12 sites.

Table S6b. Summary of Mucoromycota and Glomeromycota sequences and Operational Taxonomic Units per site (T1 only).

Table S7. Soil chemistry parameters (mean + SD) from Lycopodiella inundata subplots for all seasons combined (T0,T1,T2) compared with soil collected away from L. inundata plots.

Table S8a,b. Pairwise Spearman correlations at plot‐ and site‐level between root colonisation, soil chemistry and pollution.

Table S9. Summary of models evaluating Endogonales (Mucoromycota) arbuscular mycorrhizal fungi colonisation density per root segment using ordinal categories of none, rare, low, medium or high colonisation and relationships with covariates.

Table S10. Season‐specific models: site level correlations with Mucoromycota arbuscular mycorrhizal fungi colonisation.

Table S11. Comparison of leaf nutrition between T0 and T1.

Table S12. Covariate sets included and excluded in the models for both plot and site levels.

Method S1. Soil chemical extraction methods and additional notes on soil and root sampling.

Method S2. Other field measures and data collection.

Method S3. DNA extractions.

Method S4. Air pollution data.

Method S5. Model development assumptions.

Click here for additional data file.