Low impact of Carpobrotus edulis on soil microbiome after manual removal from a climate change field experiment

Synergic effects between climate change and invasive species may alter soil microbial diversity and functioning, as well as cause major shifts in physicochemical properties. Moreover, some of these ecological impacts may manifest even after the removal of the invasive species. We have conducted a field experiment to assess such effects on soil microbial communities (fungi and bacteria) and physicochemical properties seven months after the removal of Carpobrotus edulis (L.), an invasive plant of coastal dune ecosystems. C. edulis grew on the experimental plots for 14 months under current (“in-vasive species treatment”) and increased warming and drought conditions (“combined treatment”). Then, all plant parts (above and belowground biomass) were removed with a non-aggressive eradication method and soil samples were collected seven months later in the experimental and control plots (no invasive species and current climatic conditions). We predicted a general compositional shift in microbial communities in response to the presence of the invasive species. Moreover, given that water is the most limiting factor in this type of ecosystem, we also predicted a more pronounced compositional shift in the treatment combining invader presence and climate change. While species richness was similar amongst treatments, we observed some taxonomic and functional variation in soil microbial communities. Notably, fungal and bacterial communities exhibited contrasting responses. The species composition of bacteria differed significantly between the “invasive species” and “control” treatments, while, in the case of fungi, the most substantial difference occurred between the “invasive species” treatment and the combined treatment of “invasive species and climate change”. Some chemical properties, such as carbon and nitrogen content or pH, strongly differed amongst treatments, with the “invasive species” showing a different response compared to the other two treatments. Overall, our study suggests smaller short-term effects on the microbial community compared to soil chemical properties. Furthermore and contrary to our initial expectations, the potential impact on the soil microbiome seemed to be weaker in the face of rising temperatures and drought conditions predicted by climate change. This outcome highlights the remarkable complexity of the impact of invasive species and climate change on belowground microbial communities.


Introduction
Synergic interactions between climate change and invasive species are considered a major driver of biodiversity loss (Walther et al. 2009), with subsequent potential impact on ecosystem functioning (Pyšek et al. 2020;Isbell et al. 2023).Invasive plant species are detrimental to biodiversity because they reduce native plant diversity (McGeoch et al. 2010;Pyšek et al. 2012;Castro-Díez et al. 2016) and alter soil properties and biota (van der Putten et al. 2007;Vilà et al. 2011;Zhang et al. 2019) through aboveground-belowground feedback loops (Reinhart and Callaway 2006;van der Putten et al. 2013).Plant-soil feedbacks are critical elements of ecosystems, with plants shaping soil physicochemistry, as well as the composition and functioning of microbial communities (Wardle et al. 2004;Ehrenfeld et al. 2005).In turn, soil microbiota can affect plant growth directly through negative (i.e.pathogens and root herbivores) and positive (mutualists) interactions (van der Putten et al. 2007;Ehrenfeld 2010).Hence, disturbance-induced shifts in soil communities, such as those driven by plant invasions, may lead to critical changes in soil structure and functioning despite the seemingly high functional redundancy of soil biodiversity (Custer and van Diepen 2020;Torres et al. 2021).Some of these ecological impacts, mostly those associated with indirect effects on soil and plant communities and/or soil properties, may even persist after the removal of the invasive species (i.e.residual or legacy effects, Corbin and D'Antonio (2012)), as reported by Parsons et al. (2020).Moreover, an additional layer of complexity is introduced by climate change, which will likely further disrupt plant-soil interactions (van der Putten et al. 2016;Pugnaire et al. 2019;Trivedi et al. 2022) and may also affect the magnitude and duration of invasion legacies (Albertson et al. 2024).Therefore, synergic effects between climate change and invasive species may extend beyond promoting their establishment in new ranges (Walther et al. 2009) to influencing ecosystem responses and resilience.In a global change context, understanding invasive species effects on important ecosystem components, such as the soil microbiome, even after their eradication, is critical to forecast their shortand medium-term impacts under current and future climate conditions.
Plant invasions alter soil structure, chemistry and biota by introducing shifts in litter quality, root growth and root exudates (Liao et al. 2008;Custer and van Diepen 2020;Torres et al. 2021).These changes in microbial communities and soil properties are considered important mechanisms that facilitate the establishment and expansion of invasive species in their new ranges (Reinhart and Callaway 2006;Rodríguez-Echeverría 2010;Lorenzo and Rodríguez-Echeverría 2012;Dawson and Schrama 2016).Moreover, climate change also exerts a direct impact on soil microbiomes (Jansson and Hofmockel 2020).Raising temperatures may stimulate metabolic rates, promoting growth and diversity of soil communities, while limited water availability may have the opposite effect, reducing microbial biomass and biodiversity (van der Putten et al. 2016;Pugnaire et al. 2019;Jansson and Hofmockel 2020).Hence, an increase in aridity due to climate change is expected to reduce the diversity and abundance of soil bacteria and fungi (Maestre et al. 2015).However, our understanding on how plant-soil feedbacks may disrupt ecosystem functioning under future environmental conditions is still very limited (Classen et al. 2015;van der Putten et al. 2016;Anthony et al. 2020) and even more in the case of interactions with invasive plant species.Regarding the sole impact of invasive species on soil microbial communities, previous meta-analyses have reported both highly heterogeneous, but limited impacts on soil microbiomes (Custer and van Diepen 2020), an increase in microbial activity (Xu et al. 2022), as well as a tendency of bacterial diversity to increase with the presence of herbaceous invasive species (Torres et al. 2021).To move towards an integrative assessment of the impact of different global change drivers on the soil microbiome, field experiments simulating species invasion under increased warming and drought conditions offer great opportunities to assess the synergic impact of invasive species and climate change on soil microbial communities.
Coastal dune ecosystems are fragile and vulnerable to plant invasions (Castillo and Moreno-Casasola 1996;Muñoz-Vallés and Cambrollé 2013).One of the most aggressive invaders of these habitats worldwide is Carpobrotus edulis (L.) N.E.Br., a mat-forming species that inhibits germination and reduces light, space and nutrients for native plants (Campoy et al. 2018).While the impacts of C. edulis are locality-specific (Novoa et al. 2014;Rodríguez-Caballero et al. 2020), the invasion often results in an increase in organic matter, salinity, humidity and available phosphorus (P) and nitrogen (N) (Santoro et al. 2011;Novoa et al. 2013;Vieites-Blanco and González-Prieto 2018a;Rodríguez-Caballero et al. 2020).Concurrently, soil pH and the availability of other macronutrients, such as calcium (Ca), sodium (Na) or magnesium (Mg), typically decrease in invaded areas, although increases or no significant changes on these soil characteristics have also been reported (Vilà et al. 2006;Conser and Connor 2009;Novoa et al. 2013;Novoa et al. 2014;Vieites-Blanco and González-Prieto 2018a;Caravaca et al. 2022).These changes in soil properties and nutrient cycling have been linked to alterations in the structure and function of microbial biota (Badalamenti et al. 2016;Novoa et al. 2020;Rodríguez-Caballero et al. 2020;Caravaca et al. 2022) and, in fact, plant-soil feedbacks have been identified as an important mechanism promoting C. edulis invasion (de la Peña et al. 2010;Vieites-Blanco and González-Prieto 2018b).Furthermore, the recovery of native vegetation is hampered by legacy effects on soil abiotic properties, such as the presence of allelopathic substances in C. edulis litter (Conser and Connor 2009;Novoa et al. 2012) or the accumulation of carbon and nitrogen compounds, which may favour the establishment of ruderal species (Novoa et al. 2013).However, soil properties might recover to pre-invasion conditions (Novoa et al. 2013) and, given the strong link between soil properties and soil microbiota, it could be hypothesised that this recovery would be also mirrored by microbial communities in the soil.
In this study, we assess the effects on soil abiotic and biotic components seven months after C. edulis removal from a field experiment in which it had grown under current and future climate conditions for 14 months (see Suppl.material 1: appendix 1 for details and photographs of this experiment).A previous study showed that C. edulis growth and physiological performance increased under warming and drought, suggesting that invasion would be promoted under climate change (Campoy et al. 2021).Here, we focus on the effects on soil microbiome and physicochemistry observed several months after C. edulis specimens were removed from experimental plots.We predict a general compositional shift in microbial communities in response to the presence of the invasive species and, in particular, that this shift would be observable after the removal of the invasive species.Moreover, given that water is the most limiting factor due to reduced soil water-holding capacity in this type of ecosystem, we anticipate a more pronounced compositional shift in the treatment combining invader presence and climate change.Thus, we expect the combined effects of the invasive species and climate change would cause a stronger impact on soil microbial communities.

Experimental design
The study was conducted on an experimental field plot located over a secondary dune in the island of Sálvora, within the National Park of the Atlantic Islands (northwest of the Iberian Peninsula).From September 2015 to November 2016 (14 months), a full factorial experiment was performed in forty subplots (1.55 m × 1.75 m) to study the responses of the invasive C. edulis to climate change (Campoy et al. 2021).The duration of the experiment is four times longer than the one considered in a previous greenhouse experiment, which reported effects of C. edulis on the microbiome after three months (de la Peña et al. 2010).In detail, C. edulis was transplanted on to the experimental subplots in September 2015.At the end of this experiment (November 2016), all plant parts (above and belowground biomass) were carefully harvested (i.e.hand-pulling) and no further intervention was carried out until soil sample collection (June 2017).At the time of harvest, C. edulis was the dominant species, in terms of dry mass contribution, in the plant community, both in the invaded subplots [Inv], where it accounted for 28% of the community and in the climate change subplots [Inv_ClimCh]), where it accounted for 47% of the community.For the current study, 24 subplots were considered with the following treatments: eight control subplots (uninvaded), where spontaneous, natural vegetation had grown under current climatic conditions; eight subplots where C. edulis had grown under current climatic conditions (i.e.invasive species treatment [Inv]) and eight subplots where C. edulis had grown under modified temperature and rainfall conditions (i.e.combined treatment: invasive species + climate change [Inv_ClimCh]), according to the predictions for the study area (EEA 2017), specifically 2 °C temperature increase and 33% rainfall decrease.For detailed information of the experimental set-up, see Campoy et al. (2021) and the additional description provided in Suppl.material 1: appendix 1).
To assess differences in soil microbiome amongst treatments, we took five samples from the top 2-10 cm of soil in each subplot with a soil sampler that was disinfected with 10% diluted bleach.The top 2-cm layer was discarded in order to remove recently accumulated organic material not integrated into the soil structure.To avoid cross-contamination, the first sample from each subplot was discarded to ensure that no traces of disinfectant or soil from the previous subplot remained in the soil sampler.The other four samples were mixed to make a composite sample with which we filled two 2 ml sterilised cryogenic vials.The vials were immersed in a container with liquid nitrogen until their arrival at the laboratory, where they were stored at -80 °C until sample processing for genomic analyses.

Analysis of soil properties
In each of the 24 subplots, we also took one soil sample from the top 2-10 cm, as described in the previous section, for soil physico-chemical analyses.Samples were sieved (< 2 mm) and divided into fresh subsamples, which were kept at 4 °C for inorganic N measurements, and air-dried subsamples, which were finely ground (< 100 μm) in a planetary ball mill (Retsch PM100, Germany, with cups and balls of zirconium oxide) for the rest of the analyses.
Soil pH was measured in a 1:2.5 soil:solution ratio, both in water and 0.1 M potassium chloride (KCl), with a pH-meter (Metröhm, Switzerland).Electrical con-ductivity (EC) was determined in soil extracts (1:5 soil:water ratio) with an EC meter (Metröhm, Switzerland).Total soil C and N and the molar 15 N/ 14 N (δ 15 N) and 13 C/ 12 C (δ 13 C) ratios were determined in an elemental analyser (Carlo Erba, Milano, Italy) coupled to an isotope ratio mass spectrometer (Finnigan Mat, delta C, Bremen, Germany).Nitrogen isotope ratios were analysed because they can provide valuable tracers for biogeochemical cycles susceptible to be affected by invasion at local scales (Liao et al. 2008).To this end, δ 15 N has been proposed as an integrator of N cycling and used to show the magnitude of the effect of plant invasion in soil N dynamics and soil microbial communities (Sperry et al. 2006;Hellmann et al. 2016).Organic C and δ 13 C were also measured after calcium carbonate (CaCO 3 ) removal with the 'capsule method', using 20% hydrochloric acid (HCl) (Brodie et al. 2011).Soil 3 (Eurovector, Milano, Italy), an elemental reference material and the isotope standards IAEA-CH-6 and IAEA-CH-7, for δ 13 C and IAEA-N1 and IAEA-N2, for δ 15 N (from the International Atomic Energy Agency, Vienna, Austria) were included in each set of 10 samples to check the accuracy of the results.Nutrients available in the soil and trace elements (Al, Ca, Cu, Fe, K, Mg, Mn, Na, P and Zn) were extracted with a mixture of 1 M NH 4 Ac (ammonium acetate) and 0.005 M DTPA (soil:solution ratio of 1:5).The mixture stirred for 2 h, was filtered through cellulose paper (Filter-Laboratory 1242, Ø 90 mm) and analysed with a simultaneous ICP-OES (Varian Vista Pro, Mulgrave, Australia) (García-Marco et al. 2020).Inorganic N species (i.e.NH 4 + and NO 3 -) were extracted with 2 M KCl (1:5 soil:solution ratio) and measured by microdiffusion, as described by Vieites-Blanco and González-Prieto (2018a).

DNA isolation, amplification and sequencing
DNA was isolated using the DNeasy PowerSoil Pro DNA isolation kit (Qiagen), strictly following the manufacturer's instructions and including an extraction blank to check for cross-contamination.For fungi, a fragment of the ITS1 genomic region (of around 360 bp) was amplified using the primers ITS1F (Gardes and Bruns 1993) and ITS2 (White et al. 1990).For bacteria, a fragment of the 16S rRNA genomic region (of around 300 bp) was amplified using the primers 515F (Parada et al. 2016) and 806R (Apprill et al. 2015).Illumina sequencing primers were attached to these primers at their 5' ends.Supreme NZYTaq 2x Green Master Mix (NZYTech) was used for DNA amplification.The oligonucleotide indices required for multiplexing different libraries in the same sequencing pool were attached in a second limited-cycle PCR.A negative control that contained no DNA (BPCR) was included in every PCR round to check for contamination during library preparation.Libraries were purified using the Mag-Bind RXNPure Plus magnetic beads (Omega Biotek), following the instructions provided by the manufacturer.Then, libraries were pooled in equimolar amounts and sequenced in a fraction of a NovaSeq PE250 lane (Illumina).DNA quantification was done with Qubit (ThermoFisher).A detailed description of the DNA extraction, amplification and sequencing is provided in the Suppl.material 1: appendix 2. DNA metabarcoding analyses were carried out by AllGenetics & Biology SL (www.allgenetics.eu).

Processing of sequencing data
The obtained amplicon reads were processed using QIIME 2 (release 2021.2) (Bolyen et al. 2019). DADA2 (Callahan et al. 2016), implemented in QIIME 2, was used to remove the PCR primers, filter the reads according to their quality, denoise, merge the reads, remove the chimeric reads and cluster the resulting sequences into amplicon sequence variants (ASVs).
Non-biological DNA (primers, indices and sequencing adapters) at the reads ends was trimmed with cutadapt (Martin 2011).To remove low-quality data, fungi reads were further truncated at position 191 (forward) and at position 216 (reverse) after checking the read quality profiles.In the case of bacteria, reads were truncated at position 202 (forward) and 187 (reverse).Error rates were learned from the dataset to denoise, using the parametric error model implemented in DADA2.To reduce computation time, reads were dereplicated and the R1 and R2 reads were merged with DADA2 by overlapping a minimum region of 12 identical base pairs.As a final step, chimeras were removed by the DADA2 pipeline.The number of sequences which passed the different DADA2 processing steps are shown in Suppl.material 1: table S1.
Taxonomy was assigned to ASVs, based on the UNITE reference database (Abarenkov et al. 2020) (updated in May 2021) for fungi and the SILVA reference database (Quast et al. 2013) (updated in August 2020) for bacteria using the pretrained classifier of each respective database and the feature-classifier classify-sklearn approach (Bokulich et al. 2018) implemented in QIIME 2. To avoid misrepresented ASVs and mistagging, singletons (i.e.ASVs containing only one-member sequence in the whole dataset) occurring at a frequency below 0.01% in each sample were removed.Sequences assigned only at the kingdom level ("Fungi" or "Bacteria" and "Archaea", respectively) and the unidentified sequences were removed.ASVs present in the negative controls that were assigned to bacteria/fungal taxa, were also removed.Additionally, in the case of bacteria, non-bacterial ASVs were removed.

Statistical analyses
The final filtered ASV table was converted into a Biological Observation Matrix file (.biom) and directly imported into R 4.1.2(R-Core-Team 2021) using the package phyloseq (McMurdie and Holmes 2013).Sample 34 was removed from further analysis given that preliminary inspection of the data showed its anomalous behaviour in terms of soil properties and community composition.
To test for differences in soil properties amongst experimental treatments, we computed a multivariate analysis of variance (MANOVA) using Pillai's Trace as test statistic with the manova() function.Given the large number of soil variables, we summarised soil variables into Principal Components using function principal() with "varimax" rotation from the psych package (Revelle 2023) prior to MANO-VA analyses.When needed, we transformed variables to approach normality and only PCA components accounting for at least 70% of cumulative variance were retained.Function powerTransform() in car package (Fox and Weisberg 2019) was used to guide boxcox transformations.
To test for differences in the number of ASVs (observed richness) and, independently, in read depth amongst treatments, we fitted Negative Binomial Generalised Linear Models (NB GLM) with function glm.nb() in MASS package (Venables and Ripley 2002) for fungi data and, independently, for bacteria data.We also assessed whether richness relationships may be affected by an additive or multiplicative effect of read depth by sequentially comparing these GLM models with function anova() and a χ 2 test.Additionally, we checked that estimates of rarefied richness to the smallest read depth provided highly similar values to observed richness (Suppl.material 1: fig.S1) and analogous results in the NB GLM models.We computed rarefied richness with function rarefy_even_depth() in phyloseq package and averaged it across 100 resamples.
To test for differences in community composition amongst treatments, we computed a Permutational Multivariate Analysis of Variance Using Distance Matrices (PERMANOVA) with function adonis2() in the vegan package (Oksanen et al. 2023).To quantify differences in community composition, we first centre-log ratio-(CLR) transformed the data using function transform() in microbiome package (Lahti andShetty 2012-2019).With CLR-transformed data, we used a Euclidean distance in the adonis2() function, which is equivalent to Aitchison distance, the metric recommended for this type of data (Gloor et al. 2017).We used function gg_ordiplot() in ggordiplots package to represent differences in community composition amongst treatments and function betadisper() in vegan to assess the homogeneity of dispersions (variances) amongst treatments.To assess which treatments were significantly different, we also conducted pairwise PERMANOVAs, considering only two treatment categories at a time.We complemented the analysis of community composition with a Differential Abundance Analysis (DAA) at the genus level, in order to assess whether some particular taxa differed in their relative abundance amongst treatments.We used function trans_diff() in microeco library (Liu et al. 2020) with the ANCOMBC (Lin and Peddada 2020) and ALDEx2 (Fernandes et al. 2014) methods.Differences in ALDEx2 were assessed with the Kruskal-Wallace and the glm tests.Both ANCOMBC and ALDEx2 take into account the compositional nature of the data and provide consistent results across studies (Nearing et al. 2022).Moreover, ALDEx2 uses the CLR-transformation, as in the PERMANOVA analyses.
To assess whether soil properties would explain differences in diversity (richness and community composition) once the "Treatment effect" is accounted for, we followed a manual stepwise selection procedure on the richness models and, independently, on the community composition models.For richness models, we considered treatment as a random effect and fitted a Generalised Linear Mixed-Effects model (GLMM) for the negative binomial family with function glmer.nb() in lme4 package (Bates et al. 2015).In each step, we assessed which variable had the largest additional significant contribution to the reduced model by using function anova().For community composition models, we used the argument "by = margin" implemented in adonis2() to independently assess the marginal effects of each term to the reduced model.
To assess whether microbial functional profiles differed amongst treatments, we assigned functional traits to the obtained ASVs using FUNGuild (Nguyen et al. 2016) in the case of fungi and FAPROTAX (Louca et al. 2016) in the case of bacteria with function trans_func() in microeco package.FUNguild ranks the results within a confidence ranking as "highly probable", "probable" and "possible".Only those results marked as highly probable or probable were used, as recommended by Nguyen et al. (2016).Fungi were ascribed to one of the three main trophic groups (symbiotroph, saprotroph and pathotroph), as well as minor categories (e.g.dung saprotrophs) and bacteria were classified according to 45 functions related to energy source, N-cycling or C-cycling.Functional groups with an average relative abundance above 1% were tested for differences amongst treatments using MANOVA.
All analyses were performed in R (R Core Team 2021).The R code is provided as supporting information (Suppl.material 2).

Composition of microbial communities
We delimited a total of 1295 ASVs of fungi and 5793 ASVs of bacteria (Table 1).The mean richness per sample was 256 for fungi and 1150 for bacteria.Rarefaction curve plots suggest that samples were complete (Suppl.material 1: fig.S2), despite read depth being significantly lower in the control treatment than in the global change treatments in the case of bacteria (negative-binomial likelihood ratio stat (NB LR stat) = 7.77, p = 0.020; Fig. 1).Marginally significant differences in read depth amongst treatments were also observed for fungi data (NB LR stat = 5.90, p = 0.052).

Differences in soil properties
Variability in soil properties was summarised into six principal components, accounting for 75.5% of variance (Suppl.material 1: table S2).Variable transformations and loadings can be found in the Suppl.material 1: table S2.A MANOVA analysis showed significant differences amongst treatments (Pillai-Bartlett statistic = 1.208;F 12,32 = 4.06, p < 0.001) in PCA components related to carbon and nitrogen compounds, pH, EC and available Ca and Mn (Table 2).In general terms, the invasive species treatment was the one showing a distinct type of behaviour, while the control and the combined treatments (Inv_ClimCh) had similar soil properties (Fig. 3).Differences in raw variables of soil properties are shown in Suppl.material 1: fig.S6.

Differences in species richness and community composition
Species richness was not significantly different amongst treatments in fungi (NB LR stat = 3.48, p = 0.176) and only marginally significant in bacteria (NB LR stat = 5.75, p = 0.056) (Fig. 1).In fact, the treatment effect turned to be not significant when an additive model with sequencing depth was considered in bacteria, as  the number of reads was the only significant predictor in such model (Z = 16.39,p < 0.001, see Suppl.material 1: table S3 for details).Analogous models for averaged rarefied richness provide identical results and are shown in Suppl.material 1: table S3.
We also built a Generalised Linear Mixed-Effects model (GLMM) with treatment as a random effect in order to assess if soil PCA components had a significant contribution when included into this minimal model.In the case of fungi, only RC9 (with Ca and Mn as main variables loadings) had a significant contribution, with an increase in explained variance of 3.2% (Suppl.material 1: table S4).Soil PCA components did not significantly improve model fitting once treatment was included as a random effect in a GLMM of bacteria richness (Suppl.material 1: table S4).
PERMANOVA analyses evidenced differences in community composition amongst treatments both for fungi (F 2,20 = 1.37, p = 0.001) and bacteria (F 2,20 = 1.16, p = 0.036), although the proportion of explained variance was low in both cases (fungi = 12.0%; bacteria = 10.4%, see ordination in Fig. 4).Pairwise PERMANOVA evidenced that, in the case of fungi, the combined treatment (Inv_ClimCh) was significantly different from both the invasive (F 1,14 = 1.60, p = 0.002) and control treatments (F 1,13 = 1.35, p = 0.008), but the latter two did not significant differ (F 1,13 = 1.16, p = 0.175).Regarding bacteria, only the invasive and control treatments were significantly different (F 1,13 = 1.24, p = 0.037), as no significant differences were observed amongst the rest (Control vs. Inv_ClimCh: F 1,13 = 1.18, p = 0.073; Inv vs. Inv_ClimCh: F 1,14 = 1.07, p = 0.216).The data did not show heterogeneity in dispersion amongst treatments (betadisper: fungi: F 2,20 = 1.99, p = 0.163; bacteria: F 2,20 = 0.349, p = 0.710).Despite the observed differences in community composition, differential abundance analyses (DAA) evidenced few cases of significant differences at the genus level.In fact, ALDEx2 did not detect any significant difference, while ANCOMBC showed significant differences in eight fungal genera and three bacteria genera, all of them with low relative abundances (Suppl.material 1: table S5, figs S7, S8).We also conducted a manual stepwise procedure to assess if soil PCA components significantly contributed to explaining differences in community composition once the treatment factor was accounted for.The final model included the RC1 (mainly N and C variables) and RC2 (pH) components in the case of fungi and the RC2 (pH) and RC3 (K, Na) components in the case of bacteria (Suppl.material 1: table S4).These components notably increased the explained variance to 24.5% in fungi and 23.9% in bacteria.
In the analysis of fungal functional diversity, 55.1% (n = 714) of ASVs were assigned to a trophic group with "Highly Probable" or "Probable" confidence.Saprotrophs were the most abundant group in all treatments, followed by pathotrophs and symbiotrophs (Suppl.material 1: table S6).Wood and soil saprotrophs were the most represented guilds within saprotrophs, while plant pathogens and animal pathogens were the most abundant within pathotrophs.We observed significant differences in the proportion of reads assigned to each trophic group in each treatment (Pillai-Bartlett statistic = 1.62;F 26,18 = 2.96, p = 0.010).The relative abundance of saprotrophs was higher in the invasive and combined treatments than in the control (ANOVA F 2,20 = 6.96, p = 0.005, Fig. 5).There was also higher relative abundance of dung saprotrophs in the combined treatment than in the invasive one (ANOVA F 2,20 = 3.9, p = 0.038, Fig. 5).No significant differences were observed in the other categories (Fig. 5, Suppl.material 1: fig.S9).Moreover, two genera identified as having lower differential abundance (ANCOMBC test) in the combined treatment (Inv_ClimCh), Boeremia and Coniella, were classified as plant pathogens (Suppl.material 1: fig.S7).
In the case of bacteria, 48 functional traits were assigned to a total of 1111 ASVs (19.3% of total ASVs), although most functions had a low representation with average relative abundance < 1%, except in the case of chemoheterotrophy (mostly aerobic chemoheterotrophy), nitrate reduction and dark hydrogen oxidation, see Suppl.material 1: table S7).All bacterial functions related to N-cycling were represented albeit in low percentages.Amongst these, bacteria with nitrate reduction ability was the group with highest relative abundance, followed by nitrogen fixation.Fermentation was the most represented function within the C-cycle.No significant differences (Pillai-Bartlett statistic = 0.157; F 8,36 = 0.38, p = 0.923) were observed amongst treatments in any of the main functional categories (with average abundance > 1%, see Suppl.material 1: fig.S10, table S7).

Discussion
Our study shows taxonomic and functional variation in soil microbial communities between invaded and uninvaded treatments once the invasive plant C. edulis was removed from the field experiment with a non-aggressive eradication method.Lack of strong differences in microbial communities seven months after removal contrasted with the ones observed in some chemical properties, such as carbon and nitrogen content or pH.Uncoupled changes between the biotic and abiotic components suggest that, while microbial community shifts usually coincide with changes in physicochemical properties of the soil (Dawkins et al. 2022), the magnitude of change in the microbial community may differ from the one of soil chemical properties once the invasive species is no longer present.A tendency for recovery to pre-invasion conditions has been previously suggested, based on the dynamics of enzymatic activities (Novoa et al. 2013), but, to our knowledge, there is little information about microbiome status once the invasive species has been removed.Species richness (i.e.number of ASVs) was similar amongst treatments, but differences were found in the taxonomic and functional composition of soil microbial communities.Remarkably, fungal and bacterial communities showed contrasting responses, with the invasive and control treatments showing little overlap in species composition of bacteria while, in the case of fungi, the largest difference occurred between the invasive treatment and the combined treatment of invasive species and climate change.This result suggests that, although some species of the genus Carpobrotus are expected to particularly favour fungal communities (Badalamenti et al. 2016), an enhancement of warm and drought conditions seemed to counteract such effect, which was probably driven by the increase in water holding capacity associated with the presence of this mat-forming succulent invader (Novoa et al. 2014;Badalamenti et al. 2016).This would explain why our results contrast with the ones reported by Anthony et al. (2020), who showed that the impact of garlic mustard invasion on soil fungi would be enhanced under warmer conditions.Taken altogether, our findings indicate that C. edulis could have a significant impact on belowground diversity, with potential effects on the soil microbiome, as documented here and in earlier research (Novoa et al. 2020;Rodríguez-Caballero et al. 2020;Caravaca et al. 2022).Nonetheless, it seems that some of these effects may be ameliorated under rising temperatures and drought conditions predicted by climate change.
Shifts in microbial species composition were not associated with changes in species richness and, remarkably, drove little differences in functional diversity, mostly in the case of bacteria.Aggregate metrics, such as species richness, may mask relevant diversity changes (Hillebrand et al. 2018).On the contrary, compositional turnover metrics may detect community shifts even in the absence of changes in the number of species (i.e.species loss).This is particularly important because soil microorganisms may not be active, but persist under adverse conditions, such as drought, and recover and regrow when conditions are favourable again (Jansson and Hofmockel 2020).Hence, community composition assessments based on relative abundance metrics are more appropriate to detect such changes.Nevertheless, caution is needed in hyperdiverse systems, such as microbial communities, as compositional turnover may detect community changes, but these may still have limited consequences for ecosystem functioning.In fact, our results suggest high functional redundancy within the microbial community of this coastal ecosystem and, as such, only small differences in functional diversity of fungi were observed amongst treatments.
Saprotrophic fungi were more abundant in treatments where C. edulis had been present.Previous studies have attributed similar results to the increase in dead organic matter (Caravaca et al. 2022), given that C. edulis invasion leads to higher plant biomass (Campoy et al. 2021) and saprophytic microbes will be driven primarily by the quality and quantity of organic materials.This result would be in agreement with previous meta-analyses reporting that plant invasions tend to increase litter quality and soil N availability (Liao et al. 2008).Notably, we did not observe an increase of plant pathogens associated with C. edulis.This result contrasts to some extent with the findings of Caravaca et al. (2022), who documented an enrichment of native pathotrophs, especially plant pathogens, in the fungal community of the C. edulis rhizosphere, an outcome expected to facilitate its invasion (Caravaca et al. 2022).In general terms, our results show little functional effects of C. edulis in soil microbiomes shortly after removal, which is consistent with the small differences observed in taxonomic diversity and, in general, with the limited impact on soil microbial diversity reported in this study.Nonetheless, these results are unavoidably biased by the Raunkiaeran shortfall (sensu Hortal et al. (2015)) in our knowledge of biodiversity, as functional traits were only assigned to a fraction of taxa, implying that functional differences may exist amongst treatments, but we are not aware of them yet.Despite the major recent advances in the field of functional microbiology, this remains a barrier to our progress in comprehending ecosystem functioning and its response to global change.
While biotic differences were subtle amongst treatments, there were significant variations in soil chemistry between the invasive species treatment and the others, evidencing that some physicochemical changes may persist for a given time after the removal of the invasive species.These findings are consistent with the results of previous studies showing variations in pH, inorganic N and in the availability of several nutrients in previously invaded sites (Vilà et al. 2006;Molinari et al. 2007;Conser and Connor 2009;Novoa et al. 2013;Novoa et al. 2014;Badalamenti et al. 2016;Vieites-Blanco and González-Prieto 2018a) and provide further evidence of the high variability on the magnitude and direction of the effects of C. edulis depending on the region and habitat (Novoa et al. 2014;Vieites-Blanco and González-Prieto 2018a;Rodríguez-Caballero et al. 2020).In fact, such locality-specific responses may explain some of the differences observed here, although the most striking one, i.e. the decrease in soil total N we found in the invasive treatment, could be also attributed to the relatively short duration of the experiment (i.e. 14 months), which would have prevented litter accumulation and decomposition.Most changes in soil chemistry are likely caused by leaching from C. edulis litter (Conser and Connor 2009; Vieites-Blanco and González-Prieto 2018a) but, if plants are removed before litter is accumulated, key nutrients and cations would be removed as well.For instance, C. edulis necromass has been shown to be richer in Mg and Ca than native necromass, leading to lower concentrations of these cations in the soil (Vieites-Blanco and González-Prieto 2018a).More remarkably, our findings indicate that warming and drought may mitigate changes in soil properties driven by the invasive species, as the combined and control treatments were highly similar in soil chemistry.
Invasive plants affect soil biota through litter and rhizosphere pathways (Zhang et al. 2019), but the magnitude of their impacts, including those of C. edulis, can greatly increase with the time since invasion (Corbin and D'Antonio 2012;Novoa et al. 2012).Thus, the small short-term effects observed on the soil microbiome following the removal of C. edulis could be attributed to the relatively short duration of the experiment, even though shorter experiments (i.e.3-months) have been sufficient to alter Carpobrotus' soil microbiome in greenhouse experiments (de la Peña et al. 2010).This interpretation aligns with the idea that soil microbiota recovery is quicker when community shifts are associated with changes in relative abundances rather than the elimination of certain taxa (Corbin and D'Antonio 2012).Moreover, in well-drained and aerated dune soils, the concentration of allelopathic compounds should rapidly decrease and consequently their phytotoxicity (Bonanomi et al. 2006b).Besides, it is also well established that phytotoxic substances released during litter and organic matter decomposition may persist briefly in the soil, since they can be rapidly degraded by physical (temperature and light) and chemical processes as well as by soil microorganisms.Consequently, their inhibitory effects usually last only a few weeks and are limited to short-term phases of the early decomposition stages (e.g.Bonanomi et al. (2006a)).This could also explain why removal of invasive species in coastal ecosystems may reduce some of the effects on soil microbial communities (Parsons et al. 2020).In the case of C. edulis, Novoa et al. (2013) showed that many soil properties, as well as enzymatic activities, tend to recover to pre-invasion conditions after plant removal.Here, we also demonstrate that such recovery tendency may be accelerated in a climate change scenario.Therefore, although climate change is expected to increase growth and physiological performances of C. edulis, thus favouring its invasion (Campoy et al. 2021), the adoption of quick eradication measures may result in smaller soil microbiome shifts than at present climatic conditions.

Conclusion
Contrary to our initial prediction, here we show that climate change may reduce some of the short term effects of C. edulis on soil microbiome after removal from a coastal dune ecosystem.Nevertheless, our results also indicate that an invasion time of just 14 months is sufficient to result into differences on soil chemical properties (e.g.pH, EC and nutrient availability) and soil microbiome between invasive and control treatments seven months after plant removal.Given that effects on soil microbiome were relatively small, our findings support the low impact of commonly-used eradication measures, such as hand-pulling (Campoy et al. 2018).The removal of living biomass and litter by hand-pulling seems to be instrumental for soil microbiome preservation in sites previously invaded by C. edulis.This removal needs to eliminate both living parts and litter (Novoa et al. 2013), as well as prevent plant re-rooting (Núñez-González et al. 2021) and soil erosion (Chenot et al. 2018) to ensure the recovery of diverse native plant communities of coastal sites (Buisson et al. 2021).Our study is a step forward towards understanding potential interactions between plant invasions and climate change, revealing an unexpected complexity in the outcome of this interaction for belowground microbial communities.More field studies, including long-term surveys of restored ecosystems are needed and, ideally, should be integrated into the global framework for monitoring soil biodiversity and ecosystem function (Guerra et al. 2021).

Figure 1 .
Figure 1.Differences in read depth (a, c) and observed ASVs richness (b, d) of fungi and bacteria amongst treatments.Boxplot represents the median value and the interquartile range (IQR).Whiskers extend to ± 1.5 * IQR.Outliers (i.e.data beyond the end of the whiskers) are plotted individually.Significance letters are based on negative binomial GLM contrasts (p < 0.05).Note that the ASV richness analysis for bacteria was only marginally significant (p = 0.056) and significant differences amongst treatments disappeared when sequencing depth is accounted for.Treatment: Control; Invasive species (Inv); Invasive species and climate change scenario (Inv_ClimCh).

Figure 2 .
Figure 2. Differences in taxonomic composition (relative abundance of ASVs corresponding to different phyla [a, b] or classes [c, d]) amongst treatments.Differences in taxonomic composition amongst samples are shown in the Suppl.material 1: fig.S3).Taxonomic abundance was computed with function trans_abund() in microeco package.Treatment: Control; Invasive species (Inv); Invasive species and climate change scenario (Inv_ClimCh).

Figure 3 .
Figure 3. Boxplot of the differences in soil properties amongst treatments.Soil properties have been summarised into Principal Components after varimax rotation (Rotated Component, RC).Only components contributing to ≥ 0.7 cumulative variance after varimax rotation are retained for further analyses.Main variable loadings for each component can be found in Table 2. Boxplot represents the median value and the interquartile range (IQR).Whiskers extend to ± 1.5 * IQR.Outliers (i.e.data beyond the end of the whiskers) are plotted individually.Significance letters are based on TukeyHSD tests (p < 0.05).Treatment: Control; Invasive species (Inv); Invasive species and climate change scenario (Inv_ClimCh).

Figure 4 .
Figure 4. Ordination plot representing differences in community composition amongst treatments.Differences in community composition are computed based on Aitchison distance after Centre-Log Ratio (CLR) transformation of the data.Principal components are extracted with function rda() in library vegan and represented with function gg_ordiplot() in library ggordiplots.Treatment: Control; Invasive species (Inv); Invasive species and climate change scenario (Inv_ClimCh).

Figure 5 .
Figure 5. Differences in relative abundance of trophic groups amongst treatments.Only main trophic groups and those where significant differences were observed are shown.Boxplot represents the median value and the interquartile range (IQR).Whiskers extend to ± 1.5 * IQR.Outliers (i.e.data beyond the end of the whiskers) are plotted individually.Significance letters are based on TukeyHSD tests (p < 0.05).Functional groups were assigned with FUNGuild (Nguyen et al. 2016) and only assignments with "Highly Probable" or "Probable" confidence were considered.The variable "Dung saprotroph" was log-transformed for MANOVA analysis.Treatment: Control; Invasive species (Inv); Invasive species and climate change scenario (Inv_ClimCh).

Table 1 .
Summary statistics of the number of reads (read depth) and ASVs.

Table 2 .
ANOVA tests after multivariate analysis of variance (MANOVA, Pillai-Bartlett statistic = 1.208;F 12,32 = 4.06, p < 0.001) of differences in soil properties amongst treatments.Soil properties have been summarised into Principal Components after varimax rotation (Rotated Component, RC).Soil variables with large loading (≥ 0.70) are shown for each component.A full table of PCA variable loadings is provided in the Suppl.material1:tableS2.