Microbial community composition unaffected by mycorrhizal plant removal in sub-arctic tundra

Vegetation changes in a warming Arctic may affect plant-associated soil microbial communities with possible consequences for the biogeochemical cycling of carbon (C) and nitrogen (N). In a sub-arctic tundra heath, we factorially removed plant species with ecto-and ericoid mycorrhizal associations. After two years, we explored how mycorrhizal type-specific plant removal influences microbial communities, soil and microbial C and N pools, and extracellular enzymatic activities. Removal of ecto-and ericoid mycorrhizal plants did not change the soil fungal or bacterial community composition or their extracellular enzyme activities. However, ericoid plant removal decreased microbial C:N ratio, suggesting a stoichiometric effect decoupled from microbial community composition. In other words, microbial communities appear to show initial plasticity in response to major changes in tundra vegetation. This highlights the importance of longer-term perspectives when investigating the effects of vegetation changes on biogeochemical processes in Arctic ecosystems.


Introduction
Ongoing changes in climate are more pronounced in high-latitude and high-elevation ecosystems than elsewhere in the world (IPCC, 2021;Rantanen et al., 2022) and can alter plant community composition.Shifts in vegetation in the northern hemisphere, such as tall shrub encroachment and advancing tree lines across the tundra biome, are already visible and predicted to continue (ACIA, 2004;Åkerman and Johansson, 2008;Mekonnen et al., 2021;Myers-Smith et al., 2020).These vegetation shifts influence the soil food web and consequently biogeochemical cycles, which are mediated by the soil microbial community (Kardol and De Long, 2018;Wardle, 2006;Zheng et al., 2019).This relationship is evidenced by studies linking soil communities and their functioning to vegetation composition (Tedersoo et al., 2020).Plant roots and mycorrhizal fungi in the soil form specific mutualistic associations (Read, 1991).While plants associated with arbuscular mycorrhizal fungi are often dominant in low latitude ecosystems, shrubs and trees associated with ecto-and ericoid mycorrhizal fungi dominate many northern ecosystems (Read, 1991).Specifically, associations with ectomycorrhizal fungi are found on many trees and tall shrubs, while ericoid mycorrhizal fungi associate with dwarf shrubs, which are dominant in tundra ecosystems (Clemmensen et al., 2021;Iversen et al., 2015;Tedersoo et al., 2020).Through mycorrhizal associations, plants directly affect microbial communities in the soil (Fitter and Garbaye, 1994;Marschner et al., 2004), which can indirectly affect the soil food web further in higher trophic levels (Bomberg et al., 2010;Uroz et al., 2007).However, the ways in which changes in plant species composition affect soil microbial communities and functions are poorly understood and often context-dependent (Wardle et al., 2004).
Changes in relative abundance of plant species with different mycorrhizal associations have been observed to alter the relative abundance of functional groups of fungal species (i.e., fungal guild; Defrenne et al., 2023;Fanin et al., 2019;Martínez-García et al., 2015).Therefore, removal of plants with certain mycorrhiza-fungal associations may decrease the relative abundance of that fungal guild in the soil and consequently create an opportunity for other guilds to occupy that niche (Fanin et al., 2022;Gadgil and Gadgil, 1971).For example, removal of ericoid mycorrhizal shrubs has been shown to favor ectomycorrhizal fungi in boreal forest soils (Fanin et al., 2022), and to increase relative abundance of saprotrophic soil fungi, while decreasing overall fungal biomass (Fanin et al., 2022;Gadgil and Gadgil, 1971).However, functional niches of fungal guilds at least partly overlap (Bödeker et al., 2016;Lindahl and Tunlid, 2015) and might be less host-dependent than previously assumed (Selosse et al., 2018).For example, some ectomycorrhizal and saprotrophic fungi exist on a lifestyle continuum (Lindahl and Tunlid, 2015) and fungal species from contrasting guilds can have partially overlapping functions (Müller et al., 2020).Thus, how long it would take before these functional niches are completely taken over by other fungal guilds after abrupt vegetation shifts, and whether this indeed results in fungal community shift, remains unclear.Furthermore, although these effects are documented for temperate and boreal forest ecosystems dominated by ectomycorrhizal trees, it is yet to be determined whether similar mechanisms apply in a system that is more heavily dominated by ericoid mycorrhizal shrubs.
A prominent example of vegetation shifts in arctic tundra ecosystems, is the expansion of ectomycorrhizal deciduous tall shrubs into otherwise evergreen dwarf-shrub-dominated tundra vegetation which generally form ericoid mycorrhizal associations (Callaghan et al., 2013;Mekonnen et al., 2021).This expansion, also termed shrubification, is expected to affect soil organisms indirectly by changes in litter input quality and soil microclimate (Myers-Smith et al., 2019), but also directly by shifts in plant-fungal associations.Typical dwarf shrub species in low tundra, such as crowberry (Empetrum nigrum) and bog bilberry (Vaccinium uliginosum), associate with ericoid mycorrhizal fungi (Read et al., 2004).Trees and tall shrubs in subarctic forests, such as Betula pubescens and Betula nana, associate with ectomycorrhizal fungal partners (Read, 1991).Both ecto-and ericoid mycorrhizal fungi can efficiently degrade organic matter in the soil and thus help their host plants access nitrogen (N) and phosphorus in nutrient-poor, high-latitude environments (Kotilínek et al., 2017;Read and Perez-Moreno, 2003).In addition, both guilds can showcase a facultative saprotrophic lifestyle, however, unlike ericoid mycorrhizal fungi, certain ectomycorrhizal fungi can degrade complex organic compounds using specific oxidative extracellular enzymes (i.e.Class II peroxidases; Bödeker et al., 2009;Lindahl and Tunlid, 2015;Shah et al., 2016).Thus, ectomycorrhizal fungi may accelerate decomposition and consequently accumulate lower C stocks in the soil relative to ericoid mycorrhizal fungi (Clemmensen et al., 2015;Parker et al., 2015).Furthermore, a decline in abundance of ecto-and ericoid mycorrhizal fungi in favor of truly saprotrophic fungi may remove the control of mycorrhizal fungi on N supplies leading to enhanced soil N availability (Fanin et al., 2022;Fernandez and Kennedy, 2015;Lindahl et al., 2007.However, current empirical evidence of the long-term consequences of shrubification for the belowground community composition and functioning is inconclusive (e.g., Ahonen et al., 2021;Hagedorn et al., 2019).It remains unclear if and how fast changes in tundra vegetation cascade down to mycorrhizal and soil microbial communities and how potential changes will affect the belowground functions they mediate.On the one hand, increasingly shrubby plant communities may sequester more C in their living biomass than the current dwarf shrub tundra, as larger shrubs are more productive.On the other hand, the mycorrhizal association and rhizodeposition of taller shrubs may promote soil respiration and therefore soil C losses, resulting in an unknown ecosystem C net balance (Mekonnen et al., 2021;Parker et al., 2015).
Here, we aim to gain insights into the consequences of shrubification in sub-arctic tundra (Sweden) for plant-microbial interactions and, ultimately, soil C dynamics.We explore how soil fungi and bacteria covary with plant communities within a tundra heath ecosystem and experimentally test how changes in the plant community composition affect soil microbial communities and biogeochemistry.We investigate this by using a two-way factorial plant removal field experiment combined with a non-species-specific removal gradient to control for potential disturbance effects.The experiment simulates the shift in plant mycorrhizal type dominance expected from shrubification (from ericoid to ectomycorrhizal fungal association) by selectively removing plant species associated with ecto-or ericoid mycorrhizal fungi.Specifically, we measured initial soil bacterial and fungal community composition (i.e., before the plant removal), and investigated soil fungal and bacterial biomass and community composition, extracellular enzymatic activities, and dissolved and microbial pools of C and N after two years of mycorrhizal plant removal.We hypothesize that removing host plants of ecto-or ericoid mycorrhizal fungi will: 1. Change the soil microbial community composition by decreasing the respective mycorrhizal fungal guild in favor of bacteria and the other fungal guilds.2. Decrease the overall biomass of the soil fungal community, as previously observed (Fanin et al., 2022), despite an expected increase in the biomass of saprotrophic fungi.3. Decrease enzyme activities overall.More specifically, we expect oxidative enzyme activities to decrease after removing ectomycorrhizal plant species.4. Cause early signs of long-term changes to the soil C and N pools mediated by soil fungi and bacteria.
Testing these hypotheses experimentally will elucidate to which extent changes in tundra vegetation, as observed in ongoing shrubification, might sequentially affect plant-associated microbial communities and concomitant biogeochemical processes on a short timescale.

Experimental design
With this experiment we investigated the effects of dominant mycorrhizal associations on belowground processes in a sub-arctic alpine forest-tundra ecotone.The study was located on a south-facing slope at the forest-tundra ecotone in northern Sweden, near Abisko Scientific Research Station (ANS; 560 m a.s.l., 68 • 19′38.1"N18 • 50′53.7"E;Fig. 1).The vegetation at the study site was dominated by ericaceous dwarf shrubs such as Empetrum nigrum, Vaccinium uliginosum and to a lesser extent V. vitis-idaea, Andromeda polifolia and Cassiope tetragona, which associate with ericoid mycorrhizal fungi.Plant species with ectomycorrhizal fungal association were deciduous shrubs mostly represented by Betula nana and various Salix species.In August 2018, five spatially replicated blocks were established across approximately ha and located between 20 and 100 m apart.Within each block, seven × 2 m plots (N = 35) were placed at least 1 m apart (Fig. 1).Four plots in each block were randomly assigned one of the mycorrhizal plant removal treatments: no removal control (CTL), removal of plant species with ectomycorrhizal association (-EcM), removal of plant species with ericoid mycorrhizal association (-ErM), or combined removal of plant species with ecto-and ericoid mycorrhizal association (-EcM & -ErM).The remaining three plots per block were used for a non-specific biomass removal gradient (described below).
We measured community composition for plant species with ectoand ericoid mycorrhizal associations in all mycorrhizal removal plots before the start of the treatment.We quantified cover percentage of the to-be-removed species using a frame (1 × 1 m) with 10 × 10 cm grids in each plot (Table S1).In addition, the mycorrhizal status of all plant species observed at the site was determined based on the literature and microscopy observation of roots (Akhmetzhanova et al., 2012;Harley and Harley, 1987).Additionally, roots from plant specimens collected at the site were sampled for fungal community analyses.This confirmed that root-associated fungal communities differed significantly depending on their host plant's assigned mycorrhizal type (Supplementary Methods, Fig. S1).
The removal of plants identified as part of the different mycorrhizal types was carried out for small individuals by pulling up the plant stems, including belowground parts and roots, and for larger, woody individuals by clipping the aboveground biomass at ground level.The removed biomass was collected and dried at 40 • C until stable mass.
Plant removal is destructive, and the plant biomass removal may cause artefacts that could falsely be assigned as real treatment effects induced by our mycorrhizal plant removals.To control for this possible source of error, we established a non-specific biomass removal gradient where 5-100% of the plants were removed without considering the mycorrhizal fungal association.In each block, there was a plot of moderate (5-25 %), medium (30-60 %), and high (65-100 %) levels of removal and disturbance.In this way the 20 non-specific plant removal plots represented a continuous gradient of plant removal while the levels of removal were evenly distributed across blocks.(Monteux et al., 2024).Biomass of all vascular plants was removed based on area, an appropriate proxy for biomass (Monteux et al., 2024), using a grid with 100 cells 20 × 20 cm, where each cell represents 1 % of the plot area.The location of the grid cells to be cleared was determined pseudo-randomly, ensuring that all plant free grid cells were spread across the entire plot (see Fig. 1 for removal level in each plot; Monteux et al., 2024).
The first two blocks (1 and 2) were established in August-September 2018, at the end of the growing season.The remaining blocks (3, 4 and   5), due to practical constraints, were established later in July 2019, at which time regrowth in blocks 1 and 2 was also removed a second time.Since the initial removal at blocks 1 and 2 took place at the end of the growing season and little growth occurs over winter, July 2019 is considered the onset of the experiment.Experimental removal treatments and non-specific biomass removal were repeated yearly, and the edges of the plots were trenched with a spade and breadknife twice a year to prevent the ingrowth of roots or mycorrhiza into the plots.

Soil sampling
Soil samples were taken from all plots (N = 35) on July 18, 2019, and September 9, 2021, using apple corers (1.5 cm diameter) in the organic layer to a depth of 5-8 cm, in four cores per plot and sampling date (i.e., one sample from each square-meter in each plot).The four cores from each plot were pooled and homogenized by hand, and ca. 2 g fresh soil was subsampled to determine soil gravimetric moisture content (105 • C for 48 h) and organic matter content (SOM; loss on ignition at 475 • C for 4 h).Soil from the edges of the four soil coring holes was sampled for DNA analyses using bleach-and ethanol-cleaned tweezers and pooled into a single 1.5 ml microcentrifuge tube per plot.This tube was immediately either snap-frozen in dry ice and stored at − 80 • C (2019) or placed in a cooler box with frozen cooling elements for less than 8 h, then stored at − 20 • C (2021) for 1-3 months until DNA extraction.
To determine soil microbial biomass and activities during peak growing season, we collected another set of soil samples on July 14-15, 2021 from all plots (N = 35).In each plot, we took two soil cores (diameter 1.5 cm) in the middle of each subplot from the organic horizon, which were then pooled.These samples were homogenized (2 mm sieve) within two days and kept cool (+4 • C) until further analysis.

Carbon, nitrogen, and extracellular enzyme activities
There were slight differences in sample handling and analysis of soil chemistry variables between 2019 and 2021, and extracellular enzyme activities were measured only in 2021.In 2019, 3 g of fresh soil was set in 75 ml ddH 2 O and shaken at 100 rpm for 1 h.Water-extractable dissolved organic carbon (DOC) and nitrogen (DN) were measured by hightemperature catalytic oxidation (HTCO) using a Shimadzu TOC-V CPH analyzer with a TN unit (Shimadzu Corporation, Japan) after filtering and acidifying water extracts with 0.45 μm Filtropur S (Sarstedt AG & Co., Germany) and 50 μl 20 % HCl to 20 ml filtrate, respectively, to remove putative carbonates.In 2021, soil extractions were conducted with 20 ml 0.05 M K 2 SO 4 (2 h, 150 rpm) and filtered (Ahlström Munksjö 393), followed by freezing at − 18 • C until analysis.K 2 SO 4 -extractable dissolved organic carbon (DOC) and nitrogen (DN) were measured as described above.For the additional analysis of microbial biomass carbon (MBC) and microbial biomass nitrogen (MBN), the same extractions were performed on a different subsample after fumigation with chloroform (Brookes et al., 1985).MBC and MBN were then calculated as the difference in DOC or DN between fumigated and non-fumigated extracts.MBC and MBN values were furthermore corrected using 0.49 and 0.54 conversion factors for their respective extraction efficiencies and all were standardized as mg g − 1 dry soil.
Extracellular enzymes were extracted within three days of sampling from 2 g fresh soil and subsequently used for the determination of potential hydrolytic and oxidative extracellular enzyme activities following Criquet et al. (1999) and Jassey et al. (2011).Specifically, the potential hydrolytic decomposition of cellulose and xylo-oligosaccharides was analyzed by β-glucosidase (BG) and β-xylosidase (BX) assays.Potential hydrolytic decomposition of chitin was analyzed by N-acetyl glucosaminidase (NAG) assay, and potential oxidative decomposition of lignin and other aromatic compounds by phenoloxidase (POX), Mn peroxidase (MnPER), laccase (LA) and lignin peroxidase (LiPER) assays.The details of assay preparation and fluorescence measurements are presented in Supplementary Methods.

Amplicon sequencing
DNA was extracted using the DNeasy Power Soil Pro Kit according to the manufacturer's instructions (Qiagen, Venlo, The Netherlands).Amplicons were generated from the fungal ITS2 and bacterial 16S rRNA gene V4V5 regions in June 2022 and sequenced on a PacBio Sequel RS1 system for ITS and on an Illumina MiSeq for 16S (Supplementary Methods).Bacterial amplicon sequence variants were inferred using the DADA2 pipeline (Callahan et al., 2016), and their taxonomy assigned based on the SILVA 138.1 database (McLaren and Callahan, 2021).Fungal operational taxonomic units were inferred using the SCATA single-linkage clustering pipeline at 1.5% dissimilarity (Brandström Durling et al., 2011), and their taxonomy assigned based on UNITE 8.2 database (Abarenkov et al., 2020) and an in-house database, then fungal functional groups were assigned based on the Fungal Traits database (Põlme et al., 2020).Quality control of the sequencing and bioinformatics processing was based on negative controls as well as positive control mock communities for bacteria.More details of the molecular workflow and bioinformatic processing are presented in Supplementary Methods.

Quantitative real-time PCR
Fungal biomass was estimated using quantitative real-time PCR of the fungal ITS2 region on a BioRad CFX Connect Real-Time system in duplicates.The DNA concentration in the extracts was quantified on a Qubit Fluorometer (Thermo Fisher Scientific, Waltham, USA), and the extracts were diluted to 0.5 ng DNA μL − 1 .Each reaction contained 2 μl of diluted DNA template, 10 μg bovine serum albumin, and 1x Biorad iQ™ SYBR®Green Supermix (Bio-Rad laboratories, Hercules, USA) in 20 μl reaction volume.The primers used were the same as for amplicon preparation, and qPCR conditions are described in Table S2.The fluorescence signal was read for 5 s at the end of each cycle, at 72 • C. Standard curves (with 0.5 x 10 7 to 0.5 x 10 13 copies per reaction) were obtained by serial dilutions of linearized plasmids of the amplified gene fragment.Data are presented as ITS2 copy numbers per mg dry soil and corrected in each sample for the proportion of non-fungal (i.e., plant) sequences caught by these primers as derived from the amplicon data.Before qPCR, the absence of polymerase inhibitors was ensured in all samples by amplifying a known amount of pGEM-T plasmid (Promega, WI, USA) spiked into PCR reactions containing template and into nontemplate controls by using plasmid-specific primers.No inhibition of the amplification reaction was detected with the amount of DNA used.

Statistical analyses
The two sampling events in 2019 and 2021 represent a description of the initial state of the study system (i.e., before the removal treatment started) and the response to two years of removal treatment, respectively.As such, the two datasets were tested separately in all analyses.To ensure the absence of confounding effects in plot selection on the plant communities before the experiment started, we carried out permutational multivariate analysis of variance (PerMANOVA) with the adonis function in R package vegan (Oksanen et al., 2012) based on a Bray-Curtis dissimilarity matrix (Bray and Curtis, 1957) and using -ErM and -EcM as fixed factors (main and interaction effects) with 999 permutations, after checking for homoscedasticity (function betadisper in vegan).

Non-specific biomass removal disturbance assessment
Before analyzing the effects of the mycorrhizal plant removal treatments on microbial and soil variables, we tested the effect of the nonspecific plant biomass removal.The impact of the non-specific biomass removal on ITS and 16S communities (non-normalized count data) was assessed using a multivariate generalized linear model with negative binomial distribution, implemented in the manyglm function (mvabund, Wang et al., 2012), with the default resampling for ANOVA purposes (999 PIT-trap iterations).Linear mixed models, using non-specific biomass removed as fixed effect and block as random effect, were used for analyzing the remaining soil abiotic variables with the lmer function after checking for data normality by visually assessing the residuals and through Shapiro-Wilk tests.The amount of biomass removed in the non-specific removal plots did not significantly associate with any of the tested variables (Table S3 for microbial and Table S4 for soil variables).In other words, the disturbance of removal of plant biomass (irrespective of mycorrhizal association) did not affect any of our response parameters.We therefore consider that no adjustments of the measured variables were necessary before analysis of the treatment effect in the mycorrhizal-removal experiment (see Monteux et al., 2024).

Mycorrhizal plant removal treatment effects
To analyze the effect of -ErM and -EcM treatments on the microbial communities, we used a multivariate generalized linear model with -EcM and -ErM as fixed factors as above.We included spatial block as a random factor and used the function manyglm implemented in the mvabund package (Wang et al., 2012) with negative binomial distribution and default resampling (999 PIT-trap iterations) for ANOVA purposes.Model fits (residual distribution and normal Q-Q plot) were evaluated visually before ANOVA.In addition, we tested for spatial differences in the microbial communities between the experimental blocks using manyglm with the same parameters as above.This analysis was done separately for both fungal and bacterial communities.Soil fungal and bacterial communities as well as vegetation cover were visualized after ordination by principal coordinates analysis (PCoA) of Bray-Curtis distances for fungal and plant communities, and weighted-UniFrac distances (Lozupone et al., 2011) for bacterial communities.Correlations between soil variables and the ordination spaces were tested with the envfit wrapper (Oksanen et al., 2012).
The effects of -EcM and -ErM treatments on all univariate response variables (fungal ITS gene copy number, DOC, DN, MBC, MBN and C:N ratio, and extracellular enzymatic activities) were tested with linear mixed models including block as a random factor using lmer (Kuznetsova et al., 2017) after checking for data normality by visually assessing the residuals and through Shapiro-Wilk tests.The p values were adjusted for multiple testing (Bonferroni-holm method, N = 17).
We further explored which variables associated with the variation in the different components of our system: bacterial, fungal and plant communities as well as the enzymatic activity profile (as summarized by principal components analysis (PCA) with prcomp function).For each of those, we used the scores along the first two principal coordinate axes (PCA for enzyme activities) and assessed which variables were the most informative at explaining them among soil chemistry, enzyme activities (PCA scores), bacterial and fungal community composition (PCoA scores), and the amount of biomass removed.The large number of possible predictors in comparison to the number of samples required using a regression tree approach rather than, e.g., multiple regression or path analysis.To this end, we used the VSURF algorithm (Variable Selection Using Random Forests, Genuer et al., 2015) based on random forests (Breiman, 2001;Hastie et al., 2009).More specifically, we used the relative importance scores at the interpretation phase of a VSURF run carried out using 2000 and 25 forests for the thresholding and interpretation phases, respectively, with 10000 trees each, and the default mtry parameter for variables to compare at each branch (p/3, where p is the total number of explanatory variables).To assess the robustness of findings derived from VSURF, variable selection algorithm on partial least squares regression results was carried out using PC1 and PC2 of either 16S or ITS datasets as response variables and using block and all soil and other community variables as predictors (plsVarSel, Mehmood et al., 2012).We acknowledge the limits of this approach, especially since we are analyzing a low number of data points, and thus only use it to suggest informative variables among multiple collinear variables.
All computational work, unless specified otherwise, was done in R version 4.2.2 (R Core Team, 2022) and all plotting with package ggplot2 (Wickham, 2016).The code and processed data are stored at the Bolin Centre Database (Kirchhoff and Monteux, 2024).All sequence data are deposited at ENA under project number PRJEB64296.

Initial plant, fungal and bacterial communities
Experimental plots had homogeneous plant and microbial communities before the onset of the experiment, in accordance with the random treatment allocation within blocks.However, mycorrhizal plant and especially fungal and bacterial communities significantly differed between blocks (manyglm, block main effect: p < 0.01, Fig. S2).Random Forest analyses suggested that fungal community was the best predictor of bacterial community and vice versa, and plant community was a better predictor than spatial variation (i.e., block) for both fungal and bacterial communities.Likewise, the Partial Least Squares approach ranked fungal and bacterial PC1 first in predicting the respective other, while plant PC1 and block followed in order of importance.We thus find similar patterns of community composition differences in plant and microbial communities within one habitat type, i.e. within the same heavily Empetrum-dominated tundra heath.Furthermore, ordination spaces of both plant and microbial communities associated significantly with dissolved nitrogen (DN), and the fungal community additionally with soil moisture (Fig. S2; Table S5).This reflects the initial differences in soil conditions between blocks (Fig. S3) likely underlying those observed for plant and microbial communities.
Although these results give an indication of potential driving factors shaping bacterial and fungal communities, our experiment was not designed to untangle spatial effects from those of the plant community or abiotic factors (e.g., soil chemistry, snow depth) on the microbial communities.It is, however, reasonable to assume that species presence above-as well as belowground, affect each other (Ma et al., 2019;Tedersoo et al., 2020) beyond abiotic factors.Indeed, soil microbial and plant communities have been shown to be strongly associated (Nielsen et al., 2010;Vries et al., 2012), and plant community composition can be a better predictor for soil microbial community composition than soil abiotic factors (Chen et al., 2017;Mitchell et al., 2010;Zinger et al., 2011).This relationship has been observed over large landscape scales (e.g., Chen et al., 2017;Nielsen et al., 2010) and across different habitat types (Zinger et al., 2011).In this case, the initial state of our study system suggests that even within the same habitat type and at a local scale, soil microbial communities can be very closely associated with plant community composition.Altering plant communities in such a system may therefore have profound impacts on soil microbial communities.

Fungal and bacterial communities after 2 years of plant mycorrhizal type removal
Surprisingly, and despite the observed relationship between plant and microbial community composition, we did not find a significant effect on the microbial community after two years of ecto-or ericoid mycorrhizal plant removal (Figs. 2 and 3; Table 1).Therefore, we found no support for our first hypothesis, that removing host plants of ecto-or ericoid mycorrhizal fungi would shift the soil microbial community composition.This was likely not obscured by artefacts due to manipulation, as we observed no effect of removing biomass in general that could have counteracted a specific removal treatment effect (Table S3; sensu Monteux et al., 2024).In addition, we still observed a significant effect of the experimental blocks in both fungal and bacterial communities, as was the case prior to manipulation (manyglm, block main effect: p < 0.01, Fig. S2).The composition of soil bacterial and fungal communities, therefore, appears resistant to severe manipulations of aboveground plant community composition, at least within a two-year time frame.
Response times of plant removal effects on microbial community composition vary across experiments.In temperate grasslands, microbial community composition changed within months (e.g., Hannula et al., 2021), while studies in northern boreal ecosystems found changes on decadal timescales (Fanin et al., 2019).In contrast, northern alpine grasslands that were previously forested still lacked the expected arbuscular mycorrhizal fungi in the microbial communities after a century (Castaño et al., 2022).Similarly, a study near our site found no effect on the fungal community three years after tree girdling and subsequent decreased C supply by trees to the soil (Parker et al., 2022).These lacking responses in microbial communities to changes in vegetation composition might be owed to the large diversity of soil bacteria and the plasticity of certain plant-associated fungi that can function as saprotrophs in the absence of their host plants (Zak et al., 2019).We too, did not observe an effect of ecto-and ericoid mycorrhizal plant removal on soil microbial community composition after two years.This finding could indicate that soil microbial communities are not affected by changes in the abundance of these plant groups.Therefore, shrubification could have limited consequences on soil microbial community composition.
However, we cannot fully endorse this due to the following methodological and ecological reasons.Firstly, two potential methodological caveats prevent us from drawing this conclusion: (a) our potential inability to detect fine-scale changes in a system with high spatial variability, and (b) the stability of DNA in soils at the time scale considered here.
(a) Both fungal and bacterial communities were found to be significantly different across space (i.e., between experimental blocks).Spatial variability of bacterial communities has been explored in heathlands (Hill et al., 2016), suggesting little variability at small distances (<100m) but high variability beyond this, which corresponds to our observed differences across blocks.To our best knowledge, a formal power analysis framework was not available for the analyses we carried out, but only for distance-based methods (Kelly et al., 2015).We therefore cannot rule out that a subtle treatment effect was already present after two years but that it would have required higher Fig. 2. Relative abundances of fungal genera in each of the experimental plots before (2019) and after two years of plant removal treatment (2021).Colors indicate the respective fungal functional group or primary lifestyle (Põlme et al., 2020).Genera with an abundance <5% are combined into one category.(b) Another, likely concomitant, explanation is that microorganisms detected by our DNA amplicon approach were no longer active, or even represented necromass.Indeed, DNA is stable over long-time scales in the cold conditions of arctic soils, particularly in acidic soils with low cation exchange capacity (Carini et al., 2016).Excluding relic DNA (Carini et al., 2016) or analyzing rRNA, which degrades more quickly (Schippers et al., 2005;Weller and Ward, 1989), could provide a more adequate picture of the microorganisms that are active and alive.This might have revealed effects that were transiently unobservable in the DNA fraction after only two years of vegetation manipulation.
Additionally, a lag between above-and belowground changes is reasonable to expect from an ecological point of view (Castaño et al., 2022;Fanin et al., 2019).Only those biotic or abiotic components of the plant-soil continuum directly within the (myco-)rhizosphere and therefore closely linked to plants may initially be affected by aboveground changes.More generally, belowground ecosystem responses may only appear after longer periods of aboveground manipulations (Fanin et al., 2019;Geml et al., 2015), especially in tundra heath communities dominated by fungal networks with slow turnover rates (Gavazov et al., 2022).Importantly, the potential for ericoid mycorrhizal fungi to adjust their lifestyle towards saprotrophism, could lead to functional changes in the community without being reflected in the fungal community composition.Considering these potential explanations, we cannot conclusively demonstrate whether the absence of change in microbial community composition reflects the lack of an effect or a transient response.Unambiguously answering this question will require to re-assess this effect over the longer term and with RNA-based approaches.
While we did not see direct effects on the microbial community composition of the removal treatment, the relationships between microbial communities and soil properties changed.After the mycorrhizal type removal the fungal community composition changed from being correlated with DN and soil moisture in 2019 (i.e.before the removal treatment) to a correlation with soil moisture and DOC in 2021 (after 2 years of removal; envfit p < 0.05; Fig. S2; Table S5).The bacterial community changed from having a significant relationship with DN in 2019 to being correlated with soil moisture and SOM in 2021, respectively (envfit p < 0.05; Fig. S2; Table S5).Among the variables that were only measured after two years of manipulation, both fungal and microbial communities associated with microbial N and the extracellular activities of phenoloxydase, β-xylosidase and β-glucosidase, while only bacterial communities associated with Mn-phenoloxydase activity (Fig. 3, Table S5).However, most of the variables also exhibited trends (p < 0.1) or significant associations with microbial community composition in the non-specific biomass removal plots both before and after the removal treatments (Table S5).As both the mycorrhizal and non-specific removal treatments did not directly affect the soil properties (Tables 2 and S4) or microbial community composition (Table 1), changes in correlation between them likely stem from spatial variability or from the disturbance rather than the treatments.Removal of aboveground plant parts creates not only dead plant roots but also dead mycelium in the soil.This fungal necromass can form a new habitat for example for saprotrophic fungi (Brabcová et al., 2016;Lindahl et al., 2010).Once the fungal necromass originating from the removal of mycorrhizal plants will be depleted, the composition of microbial communities and their effect on soil C and N dynamics may exhibit more visible changes, again emphasizing the need for longer-term approaches to arctic plant-soil systems.

Microbial biomass, enzyme activities and C and N pools after 2 years of plant mycorrhizal type removal
The ecto-or ericoid mycorrhizal plant removal did not significantly affect fungal biomass as measured by ITS qPCR or the overall microbial biomass C (Table 2, Figs.S3 and S4).Consequently, we found no support for our second hypothesis which stated that removing plants with ectoand ericoid mycorrhizal association would decrease the fungal biomass.Since the microbial community composition was also not affected by plant removal, this suggests that fungi and bacteria are largely unresponsive to a drastic change in the vegetation, after two years.Furthermore, neither oxidative nor hydrolytic extracellular enzyme activities were affected by the treatments (Table 2), which falsifies our third hypothesis.In fact, bacterial Pco1 was the most important variable to describe the extracellular enzyme activities in the RandomForest analysis, suggesting that fungal community composition was a poorer predictor for enzyme activities overall.
However, removing ericoid mycorrhizal plants significantly decreased the microbial C:N ratio (main effect of -ErM, p = 0.01; Table 2, Fig. 4), since it tended to increase microbial biomass N (main effect of -ErM, p = 0.09 Table 2, Fig. S3).The non-specific biomass removal did not affect microbial C and N or their ratio (Table S4).This partly supports our fourth hypothesis, which suggested that mycorrhizal removal would lead to changes in soil C and N pools.Altogether, these findings indicate that removal of plants with ericoid mycorrhizal associations influences microbial stoichiometry.This rapid change in C:N ratio may reflect functional changes in the microbial community, which are not yet visible through DNA amplicon analysis due to relic DNA as discussed above.
As fungi have a higher C:N ratio than bacteria (Rousk and Bååth, 2007;Strickland and Rousk, 2010), the decrease in C:N ratio could be interpreted as a shift towards bacteria-dominated communities.However, we do not deem this explanation satisfactory, since the removal treatments did not affect the microbial C, fungal biomass (ITS copy number), or microbial community composition.Rather, the decrease in microbial C:N ratio after ericoid plant removal aligns with an expected lower C or higher N demand from ericoid mycorrhizal fungi (Clemmensen et al., 2013;Hobbie and Högberg, 2012).This could mean that the removal of N-demanding ericoid mycorrhizal plants decreased fungal transfer of N to their host plants (Moorhead and Sinsabaugh, 2006), and, as a consequence, fungal N-mining from recalcitrant compounds might have decreased.Consequently, the fungi may need to allocate less N for the synthesis of extracellular enzymes and, instead, store more N in their biomass.This might lead to a decrease in decomposition of recalcitrant organic compounds in the absence of ericaceous shrubs, which could potentially increase soil C stocks in the longer term (Clemmensen et al., 2015).This suggestion remains rather speculative since the potential loss of inputs from ericaceous shrubs would include ericoid mycorrhizal fungi.Even if these fungi adapt, competition from the new vegetation type and associated fungal community may outcompete them.Nonetheless, the observed change in microbial stoichiometry could be an early sign of changing C dynamics in response to vegetation change or tundra shrubification, which is not measurable immediately because it is buffered by the plasticity of the soil microbial community.The long-term consequences of these observations in relation to changes in vegetation composition remain open because of the complex dynamics and potential compensatory mechanisms within the soil microbial community.

Conclusions
Fungal and bacterial communities showed no significant response to the mycorrhizal plant removal treatments after two years.Similarly, extracellular enzymatic activities as well as soil C and N pools remained largely unaffected over the two-year period.Yet, the microbial biomass C:N ratio decreased following ericoid mycorrhizal plant removal, despite no changes in microbial and fungal biomass or community composition.The absence of change in microbial community composition or biomass may stem from technical limitations due to spatial variability or preservation of relic DNA.Conversely, the absence of community changes may be explained by a time lag between above-and belowground responses to environmental changes.The altered microbial C:N ratio in combination with a consistent lack of treatment effects on microbial community composition, microbial biomass, and other soil parameters hints towards physiological adaptations within the fungal and/or bacterial communities, following a release from nutrient limitation in the absence of ericaceous plant species.
Considering the observed and documented links between the structure of plant and soil microbial communities, we expect that, following an initial lag in the turnover of microbial biomass and relic DNA, changes in microbial community composition will eventually become evident.It is nevertheless striking that soil microbial communities exhibit such remarkable buffering capacity against as major aboveground changes as reported here, even after two years.While the implications of shrubification and tree line advances for microbial communities, biogeochemical processes and thereby the C sink and source functions of tundra heath ecosystems are still unclear (e.g.Gagnon et al., 2019;Parker et al., 2021;Zhang et al., 2013), our study shows that changes in the biogeochemical balance mediated by microbial communities could take many years in response to vegetation changes.

CRediT authorship contribution statement
Leah Kirchhoff: designed and carried out the experiment, acquired,

Table 2
Analysis of variance table for linear models of soil biotic (Extracellular enzyme activities, fungal biomass as qPCR) and abiotic (dissolved and microbial C and N).Fungal biomass was tested before (2019) and after two years (2021) of plant removal treatment, all other variables are only tested for the treatment effect after two years of plant mycorrhizal type removal.All models tested the effect of EcM and ErM removal in a two-factorial way and had experimental block as a random factor.Bonferroni-holm adjusted P values (P) below 0.05 are presented in bold, and between 0.05 and 0.1 in italics.DOC, Dissolvable organic Carbon; DN, Dissolvable Nitrogen; MBC, Microbial biomass C; MBN, Microbial biomass N, MC:N, Microbial biomass C:N ratio; qPCR ITS, Fungal biomass.analyzed and interpreted the data., led the writing of the manuscript, All authors contributed critically to the drafts and gave final approval for publication.Konstantin Gavazov: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Gesche Blume-Werry: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Eveline J. Krab: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Signe Lett: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Emily Pickering Pedersen: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Martina Peter: acquired, analyzed and interpreted the data.. Stephanie Pfister: acquired, analyzed and interpreted the data.. Maria Väisänen: designed and carried out the experiment, acquired, analyzed and interpreted the data.. Sylvain Monteux: designed and carried out the experiment, acquired, analyzed and interpreted the data., led the writing of the manuscript, All authors contributed critically to the drafts and gave final approval for publication.

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

Fig. 1 .
Fig. 1.Graphical description of the study location and experimental set-up.Panel A and B show the location of the study in northern Sweden (A) and exact placement of the study site close to Abisko and above the tree line (B).Panel C shows the distribution of the experimental blocks and treatment plots over the study site.Plant mycorrhizal type removal plots are indicated in different colors according to the treatment.Nonspecific removal plots are indicated without color and the exact removal level is given next to each plot.Panel D shows pictures of the study site and the plots after the first round of plant removal treatments for each of the mycorrhizal type and non-specific removal treatments.For the non-specific removal, exemplary pictures are shown for each removal category (i.e., moderate, medium, and high) which are represented once in each block, but with different removal amounts.Specifically, 5-25 % of plant cover was removed in the moderate, 30-60 % in the medium and 65-100 % in the high removal plots, respectively.Background Data for panel A and B are from Google Maps (2023) and Alberti (2019).

Fig. 3 .
Fig. 3. PCoA based on Bray-Curtis distances for fungal (A) and weighted UniFrac distance for bacterial (B) communities following two years of plant mycorrhizal group removal.Ellipses represent 75 % confidence intervals.Amount of variance explained in principal component (PC) axis 1 and 2 is given as respective axis descriptor.Arrows show the relationship between microbial communities and soil variables.Only variables significantly (envfit p < 0.05) associating with ordination space axes are shown.MBN, Microbial biomass N; DOC, Dissolvable organic C; POX, Phenoloxidase; BG, β-glucosidase; BX, β-xylosidase; MnPOX, Manganese-Phenoloxidase.

Fig. 4 .
Fig. 4. Microbial C:N ratio in the experimental plots after 2 years of ecto-and ericoid mycorrhizal plant removal treatment.Larger dots are mean, and whiskers are SE, smaller dots are individual plots, n = 5.

Table 1
Analysis of deviance table of the multivariate linear model for the two-factorial effect of Ericoid-and Ectomycorrhizal plant removal on bacterial (16S) and fungal (ITS) community composition before and two years after the onset of the removal experiment.