Legacy effects of invasive plant species on soil bacterial community assembly, β-diversity, and ecological interactions

Abstract

Invasions of exotic plant species pose a serious threat to local biodiversity and ecosystem functioning, with their effects on soil persisting even after removal. In a mesocosm experiment, we investigated the impact of two alien species, Conyza bonariensis (annual) and Solanum elaeagnifolium (perennial) on soil bacterial community after one year of growth (conditioning sampling), and their legacy effects on the bacterial community developed during the subsequent growth of a native grass species, Cichorium intybus (legacy sampling). We assessed the effects of these species by analysing soil enzymatic activity, bacterial community biomass and structure, β-diversity and the co-occurrence patterns of microbial members. Plant identity did not affect enzymatic activity, bacterial biomass and community composition. The communities across all treatments were dominated by the phylum Firmicutes particularly the Bacillus genus. The heterogeneity in the composition of bacterial communities between treatments (β-diversity) was higher at conditioning compared to legacy sampling while the niche width of the bacterial members expanded after C. intybus growth. β-diversity in soils with S. elaeagnifolium legacy was mainly driven by stochastic processes such as ecological or genetic drift while in soils with C. bonarienzis legacy, deterministic processes like environmental filtering played a dominant role. Regulation of microbial co-occurrence patterns was nearly equally influenced by stochastic and deterministic processes. However, the legacy effects of the invaders significantly impacted the robustness of bacterial networks to further disturbance, with the networks in C. bonarienzis exhibiting enhanced robustness. Our results suggest divergent management strategies for these two species: precautionary containment for S. elaeagnifolium vs. direct intervention for C. bonariensis.

Graphical abstract

Keywords

rank abundance models / S. elaeagnifolium / C. bonarienzis / soil enzyme activity / PLFAs

Highlight

	● No invasive-specific effects on bacterial community composition and biomass.
	● The identity of invaders’ legacy impacted the robustness of bacterial networks.
	● Legacy of C. bonarienzis exerts strong filtering effect on soil bacterial β-diversity.
	● S. elaeagnifolium legacy promotes the stochastic regulation of bacterial β-diversity.

1 Introduction

Biological invasions of exotic species rank among the most serious threats to local biodiversity and ecosystem functioning (Mooney and Hobbs, 2000). They can threaten the diversity or abundance of native species, destabilize invaded ecosystems, undermine economic activities dependent on these ecosystems (e.g., agricultural), and even pose risks to human health (Galil, 2007). The global economic impact of biological invasions is equivalent to that of natural catastrophes. From 1980 to 2019, financial losses due to invasive alien species amounted to $1208 billion (US), compared to nearly $1914 billion in losses caused by storms, $1139 billion by earthquakes and $1120 billion by floods (Turbelin et al., 2023). For instance, Solanum elaeagnifolium invasion affects negatively crops causing losses up to 75% in yield. These losses include decreased forage quality on grazing lands, decreased cropping land and amenity values of public space, increased water loss, increased water conveyancing costs and increased forest restoration costs (Uludag et al., 2016). Other examples of plant invaders with strong negative impact on human health were those of the common ragweed—Ambrosia artemisifolia L. (Knolmajer et al., 2024) and Cortaderia selloana (Liendo et al., 2023). Both species pose significant risks to human health due to the allergenic properties of their pollen. Climate warming is projected to expand the distribution of A. artemisifolia in north Europe causing allergy in about 44 million more people in the future (Knolmajer et al., 2024).

As plant invasions have become a global phenomenon, research has increasingly focused on their implications, as evidenced by the growing published work (Simberloff et al., 2013; Kumar and Singh, 2020; Tamburello and Litt, 2023). It has been suggested that invasive plant species introduce mechanisms that are novel to indigenous communities, potentially affecting the structure and functions of soil microbial community (Duda et al., 2003) as well as the plant-microbe interactions (Vitousek and Walker, 1989). One such mechanism involves changes in the production of litter and/or root exudates. Exotic plants, being different in the quality or quantity of above- and below-ground litter compared to native plants, can alter resource availability for belowground communities. Additionally, the release of novel chemicals with antimicrobial properties and the introduction of novel nutrient acquisition strategies, such as nitrogen fixation, can change critical biochemical processes (Hierro and Callaway, 2003; Callaway and Ridenour, 2004). Changes in root architecture or function can further induce differences in the local soil environment (van der Heijden et al., 1998). These impacts may arise from any one of these mechanisms or from multiple mechanisms acting in concert. For instance, changes in litter quality, root exudates, and nutrient acquisition strategies drive shifts in biochemical cycles—such as nitrogen cycling or nitrification activity—mediated by soil microbes. Frequently, these shifts give invaders a competitive advantage over native species.

The changes in abiotic and biotic soil properties caused by exotic plants are often referred to as “legacy effects”, a term coined by Connell (1978). These effects can persist in the soil for some time after the exotic plant species has disappeared, depending on the nature of the induced changes in soil conditions (van der Putten, 2003). Invasive species can exploit these legacy effects through self-induced modifications of soil conditions, whether biological, chemical or physical (Corbin and D’Antonio, 2012). These legacy effects can also potentially have long-term consequences for plant community diversity and productivity (Grman and Suding, 2010; Kulmatiski et al., 2011). Although numerous studies have examined how legacy effects influence the productivity and diversity of plant communities, relatively little is known about their impact on soil microbial communities. A meta-analysis by Torres et al. (2021) revealed that invasive plant legacies can affect bacterial diversity, but the influence on nitrifying bacteria and consequently on soil nitrification can vary from positive to negative or even remain neutral. According to Afzal et al. (2023), the legacy effects of invasive plants on soil nitrification are context-dependent and can differ significantly depending on the location, the invasive species involved, and other environmental factors. Consequently, a gap remains in our understanding on how invasion legacy effects shape soil microbial communities and their role in plant–soil feedback interactions.

During the phase that invasive plants establish in new areas (the so-called “conditioning” phase), soil microbes play a highly significant role. As shown in our previous work, the invasion with Solanum elaeagnifolium (a perennial plant) and Conyza bonariensis (an annual plant) differently affected the biomass of Gram-positive and Gram-negative soil bacteria in relation to unplanted soils (Kapagianni et al., 2021), indicating shifts in carbon availability. In this study, we extend our investigation by examining the legacy effect of these two invasive plants on the bacterial community in soils subsequently planted with a native plant species (Cichorium intybus). Specifically, after a year of invader growth, Cichorium intybus was grown in “conditioned” soils, and changes in attributes of the soil bacterial communities (such as diversity, community composition, and bacterial interaction networks) were assessed. This approach goes beyond merely quantifying aspects of bacterial community organization such as composition, and further explores the mechanistic underpinnings of community assembly, such as stochastic (ecological drift genetic drift, limited or unlimited dispersal) and deterministic (niche-based) processes.

Deterministic and stochastic processes, though conceptually at opposite ends of a continuum, occur simultaneously, with one typically dominating under specific circumstances (Stamou and Papatheodorou, 2023). Because community assembly is dynamic, the primary question is not which mechanism operates, but rather when and why one mechanism takes precedence over the other (Langenheder and Lindström, 2019). Deterministic processes tend to prevail in extreme (limited availability of resources or extreme abiotic conditions) and/or more heterogeneous environments, where the environment acts as a filter on microbial communities. Stochastic processes often dominate in more homogeneous environments. For instance, stochasticity in soil microbial communities predominated in deeper soil layers due to the less heterogeneous environment compared to topsoil (Luan et al., 2020). Similarly, nutrient-rich environments, such as those receiving N additions, may also favor stochasticity (Zhou et al., 2022). Invasive plants may alter resource availability through litter quality or novel compounds’ release, and through their root systems could modify soil heterogeneity. All these factors affect biotic interactions among microbial members, shifting community assembly forces between stochastic and deterministic processes (Li et al., 2022). In this study, we additionally aimed to identify potential relations between community assembly attributes and the functionality of the bacterial communities, using as proxy of functionality the activity of six enzymes involved in C, N and P cycles. Although, the study of Luan et al. (2020) provided evidence that bacterial community assembly may influence carbon metabolism in deeper soil layers, the underlying mechanisms behind this positive coupling still need further investigation. Specifically, we hypothesize that invasive plants increase β-diversity in bacterial communities via plant-specific modifications of soil heterogeneity, while subsequent native growth homogenizes conditions, thereby reducing β-diversity and shifting the balance between stochastic and deterministic processes. Furthermore, we predict that the annual C. bonariensis will impart a more disruptive legacy on soil bacterial networks than the perennial S. elaeagnifolium, making the former communities more vulnerable to future disturbances.

2 Materials and methods

2.1 Experimental design

In a mesocosm experiment, we investigated the impact of two alien species, Conyza bonariensis (annual) and Solanum elaeagnifolium (perennial) on soil bacterial community, as well as their legacy effects on the bacterial community developed during the subsequent growth of a native grass species, Cichorium intybus.

Soil used for mesocosms’ construction was collected from the topsoil (20 cm) from an area disturbed by grazing of Mount Chortiatis, Northern Greece. Soil was transferred in plastic bags, sieved (1-cm mesh) to remove coarse fragments, homogenised by hand mixing and placed in mesocosms (pots) in a glasshouse. The soil had the following properties: 51% sand, 24% clay, 25% silt, pH 5.5, organic carbon 2.96% dry soil and total nitrogen 1.7 mg g⁻¹ (Kapagianni et al., 2021).

The invasive plant species (S. elaeagnifolium and C. bonariensis) were grown individually in pots (constructed of PVC; 14 cm diameter, 15.5 cm height). Seeds of C. bonariensis and S. elaeagnifolium were collected from ex-arable fields. Five pots were sown with seeds of S. elaeagnifolium, five with seeds of C. bonariensis and five pots contained only soil, serving as controls. One plant was allowed to grow in each pot for 12 months, remaining until the natural withering and death of C. bonariensis. Any additional emerging weed seedlings were manually removed. Both the planted and control pots were watered regularly, with soil moisture maintained at 15% (v:v) by weight throughout the experiment (Bezemer et al., 2006), ensuring no leaching or drainage occurred from the pots. At the end of the 12-month period, a destructive sampling was conducted, referred to hereafter as “conditioning sampling”. At this sampling, C. bonariensis had died and was withered while S. elaeagnifolium produced fruits (berries), growth had ceased, and plant dormancy had begun. The samples’ codes at conditioning sampling were “CONT”: control, “CO”: C. bonariensis and “SOL”: S. elaeagnifolium.

After the conditioning sampling, the plants were removed from the pots. The soil in the pots was first homogenized according to the conditioning type and sifted to remove any remaining plant material (e.g., roots). The conditioned soil was then placed back into the same pots. All pots containing soil conditioned by C. bonariensis plants, S. elaeagnifolium plants, and control soil were sown with seeds of C. intybus (five pots per conditioning type). Similarly to the conditioning phase, one plant of C. intybus was allowed to grow in each pot, watered, and left to remain in the pots until it completed the vegetative stage, six months after establishment. The pots were watered regularly to keep soil moisture at 15% (v:v). After six months, we conducted the second destructive sampling, referred to hereafter as “legacy sampling”. The sample codes at legacy sampling were “PRCONT”: control legacy, “PRCO”: C. bonariensis legacy and “PRSOL”: S. elaeagnifolium legacy. The number of replicates in each treatment at both experimental phases was five. During the experiment, pots were frequently randomized within the glasshouse to account for any environmental variation. The glasshouse received natural light and the average air temperature and relative humidity ranged within 0–40 ^oC and 20–70% respectively, during the experiment.

At both samplings, soil samples were collected from the plant rhizosphere, the narrow area of soil that is influenced by plant roots. Plants were uprooted and loosely adhering soil was removed from the roots by gentle shaking for the collection of rhizosphere soil in sterile bags. We divided rhizosphere samples into two subsamples; one was stored in a refrigerator at 4 ^oC until the determination of the soil biological and chemical parameters, while the other, used for molecular determination was placed in −80 ^oC untill analysis. For DNA extraction, three out of the five samples from each treatment were used.

2.2 Bacterial biomass

Extraction and analysis of PLFAs from soil samples were performed within a week postsampling. The overall procedure involved the extraction of lipids, separation of phospholipids by column chromatography, methylation of esterified fatty acids in the phospholipid fraction and chromatographic separation and identification of the main components on a Trace GC Ultra gas chromatograph (Thermo Finnigan, San Jose, CA, USA) coupled with a Trace ISQ mass spectrometry detector, a split-splitless injector and an Xcalibur MS platform (Kapagianni et al., 2021). Specific PLFAs were used as indicators of bacterial biomass. These included i-15:0, a15:0, 15:0, i16:0, i17:0, 16:0, 17:0 (indicators of Gram-positive bacteria; McKinley et al., 2005; Rillig et al., 2006), 16:1w7c, 16:1w9c, 16:1w9t, 2OH16:0, 18:1w9t and cyc17:0 (Gram-negative bacteria; Zhang et al., 2010), and 10Me16:0, 15 10Me17:0, 10Me18:0 (actinomycetes; White et al., 1996). All of these were of bacterial origin and were chosen to represent bacterial biomass.

2.3 Enzyme analyses

The activity of six enzymes mediating the N, P and C cycles was assessed on <2 mm field moist samples. These enzymes are b-glucosidase (BG), peroxidase (POD) and polyphenol oxidase (PPO) involved in C-cycle, N-acetylglucosaminidase (NAG) and urease involved in N-cycle, and acid phosphomonoesterase (AP) that mediates P cycle. We determined enzyme activities (except urease activity) according to the procedures of Allison and Jastrow (2006) as modified for 96-well microplates (Kapagianni et al., 2021). Urease activity was determined according to the method of Sinsabaugh et al. (2000).

2.4 Microbial community composition

The DNA extraction was performed using the NucleoSpin Soil extraction kit (Macherey-Nagel, Düren, Germany) according to the manufacturer’s instructions using 250 mg of soil. Quantification of the extracted DNA was performed using a Quawell Q3000 microvolume spectrophotometer (Quawell Technology, San Jose, USA). Library preparation was performed using the 16S Barcoding Kit 1–12 (Oxford Nano-pore Technologies, Oxford, UK) following manufacturer’s instructions. The full-length 16S rRNA bacterial gene was amplified through PCR using the kit’s barcoded primers (27F: 5′-AGAGTTTGATCMTGGCTCAG-3′ and 1492R:5′-CGGTTACCTTGTTACGACTT-3′) allowing multiplexing. Amplification was conducted using 25 μL of repliQa HiFiToughMix mastermix (Quantabio, Beverly, MA, USA), 10 μL input DNA (10 ng of genomic DNA in 10 μL nuclease-free water), 1 μL of 16S barcode and 4 μL of nuclease-free water. The thermocycler was programmed with the following conditions: initial denaturation at 95 °C for 1 min followed by 25 cycles of denaturation (95 °C, 20 s), annealing (55 °C, 30 s), and extension (65 °C, 2 min), with a final extension at 65 °C for 5 min and termination of the reaction at 4 °C. The barcoded PCR products were visualized using gel electrophoresis to check the size and purity of the fragments. PCR products were then cleaned and purified with Agencourt AMPureXP magnetic beads (Beckman Coulter, CA, USA). The concentration of purified DNA amplicons was determined using the microvolume spectrophotometer. All samples were pooled equitably (100 ng of DNA) into a single library along with the reagents provided by the 16S Barcoding Kit. The library was then loaded and sequenced on the SpotON FlowCell (MinION R9.4.1) which was controlled through the MinKNOW software v.3.6.5.

Nanopore sequencing data were base-called, demultiplexed, and trimmed to remove barcode and adapter sequences using Guppy v3.2.10 (Oxford Nanopore Technologies). Any residual adapter sequences were subsequently detected and removed via Porechop v0.2.4 (Wick RR, see the website of github.com), and chimeric reads were discarded through the “discard_middle” option. The dataset then underwent further refinement with Chopper v0.2.0 from the NanoPack2 toolkit (De Coster and Rademakers, 2023), which involved: (a) cropping the first 25 nucleotides (–headcrop 25) based on initial FastQC v0.12.1 (Andrews, 2010) results indicating lower quality in the first 20 bases, (b) filtering reads to lengths between 1000 and 1600 nucleotides (–minlength 1000 –maxlength 1600), and (c) retaining only those sequences with a mean quality score ≥ 7 (–quality 7). Quality assessments of read lengths and Q scores were performed with MinIONQC v1.4.2 (Lanfear et al., 2019). Taxonomic classification was carried out using Kraken2 v2.1.2 (Wood et al., 2019) with the SILVA SSU v138.1 database (Quast et al., 2013) at a 90% confidence threshold. Estimated abundances were then calculated using Bracken v2.8 (Lu et al., 2017), and operational taxonomic units (OTUs) read counts were normalized to the median sequencing depth per sample. Only OTUs that constituted at least 0.01% of the total reads were selected for downstream analyses.

2.5 Statistical analyses

To detect changes between treatments in soil nutrient content, bacterial biomass, potential activities of six enzymes and bacterial niche width, we applied a one-way analysis of variance. In case of deviation from normality and the lack of homogeneity of variances, we applied a non-parametric test (Kruskal-Wallis). When the ANOVA test revealed significant differences, we further applied a post-hoc test. For parametric data we applied the Tukey test, while for non-parametric data the Bonferroni test. To reveal differences in relation to the overall soil enzymatic profile, a non-metric multidimensional scaling (NMDS) analysis was applied in R using the “vegan” package. Additionally, we applied a principal coordinates analysis (PCoA) based on Bray-Curtis index on bacterial community data using the “microeco” package. For the estimation of niche width, we used the Levis index of niche overlap in the “spaa” package in R (version 4.4.2).

2.5.1 β-diversity

To explore the biotic heterogeneity of the bacterial communities among treatments, the β-diversity (βSOR, sensu Sørensen) was estimated using the “beta.pair” function of the script “betapart” (Baselga and Orme, 2012) and run in R (version 4.2.0). To reveal mechanisms underpinning the variation in diversity among treatments, total β-diversity (β_SOR) was partitioned according to Baselga and Orme (2012) into the gradient component (β_grad) and the balance component (β_bal) produced respectively by stochastic and deterministic processes. The gradient component (β_grad) accounts for variation in the abundances due to the gain or loss of taxa, while the balance component (β_bal) indicates replacement of some taxa by others due to environmental sorting or spatial and historical constraints. To determine which component was more important (Xiong et al., 2021), the β-ratio index was calculated as the proportion of the gradient component to the total β-diversity. For β-ratio < 0.5, β-diversity is determined mainly by determinism, while for β-ratio > 0.5 stochasticity predominates (Xiong et al., 2021).

2.5.2 Rank abundance models

To explore whether the organization of the bacterial communities was governed primarily by stochastic or deterministic processes, five rank abundance models based on niche perspectives (pre-emption, broken stick, log-normal, Zipf and Zipf-Mandelbrot), as well as the Zero-Sum Multinomial (ZSM) model based on neutral perspectives (null model) were fitted on data from different treatments. The five former models were fitted using the function “radfit” (Gausian family), part of the Vegan package in R (Oksanen et al., 2017), while the fitting of the ZSM model was performed using the Tetame2 software (Jabot et al., 2008). The two parameters of the ZSM model, theta and dispersal rate were estimated according to the formula of Etienne. The goodness-of-fit of the models was compared using the Akaike information criterion (AIC), calculated with the equation AIC = −2 × loglikelihood + 2 × the number of parameters in the fitted model. The lower the AIC value the better the fit of the model.

2.5.3 Network analysis

The subjects of the network analysis were matrices depicting bacterial taxa co-occurrence patterns across different treatments. Taxa were represented by nodes, and the strength of co-occurrence among them was indicated by ties. Prior to network construction, the data underwent a filtering process to ensure that individual taxa were present at least in two out of three replicates. The sum of the filtered taxa was retained, maintaining taxon proportions, and the data were normalized by computing the relative abundance of each taxon. Co-occurrence patterns were determined using the CoNet software (Faust and Raes, 2016), which is an add-on to Cytoscape 3.9.1. software (Shannon et al., 2003). The methodology employed a set of metrics (Pérez-Valera et al., 2017) to measure the strength of connections, aiming to reduce the impact of relying on a single measure by requiring at least two supporting measures to confirm co-occurrence between two taxa. For the depiction of the co-occurrence patterns, we used Cytoscape 3.9.1 (Shannon et al., 2003) and UCINET 6 (Borgatti et al., 2002).

Among other network metrics estimated (Stamou et al., 2024), we assessed modularity (Q) and the Small-Worldness Index. For modularity, we used the Girvan-Newman optimization algorithm. Q values > 0.4 indicate a modular architecture (Luo et al., 2019). A network is characterized as Small-World network when it displays a similar shortest path and a much higher clustering coefficient compared to a random network. Consequently, when the ratio S (Index of Small-Worldness) = (Cl_real/Cl_null)/(L_real/L_null) exceeds 1, the network is classified as a Small-World network (Humphries and Gurney, 2008).

The network’s robustness against non-targeted disturbances (Critical function) was assessed using random graph theory. (Liu et al., 2017). Specifically, we estimated the p_s^r index, which posits that the average sum of the squares of the degrees surpasses twice the average degree of the network (K²/K > 2). The p_s^r index quantifies the critical fraction of node removal that triggers the network’s disintegration:

p_s^r=1-\frac{1}{K_o-1}

where K_o = K²/K, K is the average degree of the network, and K² is the average of the square of degrees. The larger the positive value of p_s^r, the more robust the network is, while the negative values pertain to disintegrated networks.

Another index designed to account for disturbances targeted at abundant taxa relies on examining the eigenvalues of the Laplacian and/or adjacency matrices (Liu et al., 2017). In this study, the natural connectivity of the networks was calculated as follows:

\overline{\lambda }=\mathrm{l}\mathrm{n}\left(\frac{1}{N}\right)\sum _{i=1}^N{e}^{{\lambda }_i}

where λ_i is the ith eigenvalue of the adjacency matrix. The larger the value of \overline{\lambda }, the more robust the network.

To explore the interplay between stochasticity and determinism in the network connections, the effective information index (EI) was computed following the method outlined by Klein et al. (2022). Initially, the tie strengths of each variable i were individually averaged, and the Shannon entropy representing uncertainty related to these averages was determined, referred to as the determinism component, \mathit{H}(\overline{W_i}). Subsequently, the Shannon entropy was computed separately for each variable, and the resulting individual Shannon indices were averaged to determine the stochasticity component, \overline{H(W_i)}. The effective information index (EI) was then derived as the difference between (\overline{W_i}) and \overline{H(W_i)}:

EI=H(\overline{W_i})-\overline{H(W_i)}

where H is the Shannon diversity and W_i is the strength of the node i. Finally, EI was network-size-normalized as follows:

Effectivness=\frac{EI}{{\mathrm{l}\mathrm{o}\mathrm{g}}_2\left(N\right)}

where N is the network size. Positive values indicate the prevalence of stochasticity, whereas negative values indicate the prevalence of determinism.

3 Results

3.1 Soil nutrients

Soil nutrients were recorded only at the end of the conditioning phase. During this phase, soil organic C, organic N, and the availability of K, P, NH₄ and NO₃ were estimated. Significant differences between control and conditioned pots were recorded for the concentrations of P, K and NO₃ (Fig. S1). However, soil K was significantly differentiated between Conyza and Solanum soils, while NO₃ was different only between control and solanum soils. No differences were recorded in the soil pools of C and N, nor in the NH₄.

3.2 Bacterial biomass and enzymes

Although there was a tendency for a decrease in bacterial biomass at legacy sampling, no significant difference was observed between the two samplings (Fig.1). Analysis of variance applied to AP, NAG and BG values revealed a significant effect of treatments on their values (Fig. S2); however, these differences were not consistent between control and conditioned pots or between the conditioning and the legacy samplings AP values in CONT and PRCONT pots were lower than in the other treatments in both samplings. For NAG values there was a tendency for a decline from conditioning to legacy sampling, but this was not statistically significant. PPO and POD are enzymes mainly involved in the degradation of more recalcitrant compounds. Significant differences in POD activity were recorded between the conditioning and the legacy sampling for both control and treated soils (Fig. S3a, S3b) whereas for PPO the differences were found only in the treated pots (Fig. S4a, S4b). Higher PPO and POD were recorded at the legacy sampling. Urease activity displayed significant differences between control and Conyza soils only at the conditioning sampling (Fig. S3c).

To assess the clustering of samples of different treatments in relation to their enzymatic profiles, an NMDS analysis based on Bray–Curtis distance was applied to the data (Fig.2). A clear distinction was observed between the two samplings as well as between the treated and the control soils. The distinction was greater at the conditioning sampling, while the identity of invader seemed to exert a marginal effect on soil functioning in both samplings.

3.3 Bacterial community composition

In total, 124 bacterial taxa were identified. The dominant phylum was Firmicutes exhibiting a relative abundance of 80-90% across all samples (Fig.3). The relative abundance of Proteobacteria was on average 10%, while that of Bacteroidota, Acidobacteriota and Chloroflexi was less than 10%, and that of Actinobacteriota was < 1%. The dominant genus was Bacillus, followed by Pseudomonas and Chitinophaga. Bacillus represented 82% to 95% of the total genera (Fig.4). To detect significant changes in community composition between treatments we applied one-way PERMANOVA based on Bray–Curtis distance. There were no significant differences between treatments in the same sampling point as well as between the conditioning and the legacy sampling (F=2.095, p=0.07).

3.4 Niche width

According to Fig.5, at the conditioning sampling the mean niche width of the bacterial community in CO and CONT soils was similar, while that of SOL soils was differentiated from the control. The SOL community exhibited the lowest value (a community consisting mainly of specialists). After the chicorum’s growth there was no difference in the mean niche width of the communities originated from the three treatments, although these exhibited higher values compared to those recorded at conditioning. To detect relationships between enzymes and mean niche width we applied Spearman correlation. The Spearmans’ rs (for non-normal data) between mean niche width and mean PPO and POD were 1 and 0.96 respectively, while no significant correlations between the other enzymes and niche width were recorded.

3.5 β-diversity of bacterial communities

According to the PCoA analysis, the biotic heterogeneity (β-diversity) was higher in conditioned pots (including the control) and became much lower in all pots after the C. intybus cultivation (Fig.6). The presence of the plant eliminated the heterogeneity induced by the legacy of the alien species or by the absence of plants.

Further, we detected differences concerning the mechanisms regulating the changes in β-diversity. The β-ratio (β_grad/β_glb) that describes the balance between stochastic and deterministic regulation of β-diversity, was <0.5 in CO pots in both samplings, indicating the predominance of deterministic mechanisms (Fig.7). In contrast, in SOL pots immediately after conditioning, the β-ratio was slightly above 0.5 showing a balance between determinism and stochasticity. However, after C. intybus cultivation (PRSOL), the regulation became mainly stochastic. Control pots exhibited the opposite trend: after a year of no-cultivation (CONT), β-diversity was regulated by stochastic processes (β-ratio > 0.5) while after cultivation (PRCONT), regulation was mainly deterministic. In summary, after conditioning, the regulation in β-diversity in CO pots was dominated by deterministic processes, in SOL pots there was a balance between determinism and stochasticity while in CONT pots stochasticity predominated. After C. intybus cultivation, regulation in PRCONT and PRCO was mainly deterministic, while in PRSOL it was primarily stochastic.

3.6 Rank abundance models

Five different models were applied to the rank abundance data of the bacterial taxa. Four of these are based on deterministic specifications, and one is based on neutral ones. According to the AIC values (Tab.1), two models showed the best fit: one determinist (Zipf-Mandelbrot) and one neutral (ZSM; zero sum model). The g value of the Zipf-Mandelbrot model describes the sharing of abundance between taxa; low values indicate higher evenness. g values represent ecosystem predictability, as they are related to the amount of information required to describe the system. In the same model, parameter b indicates environmental heterogeneity; higher b values are indicative of many environmental alternatives (high environmental heterogeneity). Analysing the g and b values of the model (Tab.2), revealed that except for a few taxa (two or three), the rest of the taxa attained similar abundances in a relatively homogeneous environment. This trend was more pronounced in SOL pots where b values were near zero in both samplings. An exception was found in the g (−9.96) and b (1.92) values in CO pots at conditioning which showed that in these pots a few taxa such as Bacillus, Pseudomonas, and Sphigomonas exhibited high abundances in a relatively heterogeneous environment. This suggests that these dominant taxa exploited most of the available resources. Conversely, after C. intybus cultivation, the distribution of abundance between bacterial taxa became more similar across all treated soils while the environment became more homogeneous in PRCONT and PRCO pots compared to the conditioning phase. High evenness and low heterogeneity resulting from completed interspecific competition could be indicative of a community consisting mainly of bacteria displaying specificity.

By examining the parameters of the neutral model (m=dispersion rate and I/J= rate of replenishment where J is the number of existing taxa and I is the number of new individuals, we concluded that the dispersion rate and the rate of replenishment were very low in all treated pots (Tab.2). Also, the limitations to dispersion did not differ substantially between samplings. During the conditioning phase, the limitation in dispersal in CONT pots was likely imposed by specific soil characteristics, as these pots remained unplanted for a year. In contrast for CO and SOL pots, the limitations were initially due to the legacy of alien plant species, and subsequently, the limitations were due to the combined legacy of both alien and the native plants.

3.7 Co-occurrence patterns

Bacterial networks consisted of fewer subgroups at the conditioning than at the legacy samplings in all soils (Tab.3, Fig. S4). This difference was more pronounced in SOL pots, where the network consisted of a single giant network at the conditioning sampling. At this stage, all networks exhibited SW properties (SW>1) and non-modular architecture (modularity < 0.4). The growth of C. intybus resulted in a slight increase in network fragmentation. At legacy sampling, the number of nodes increased, while the average neighborhood size decreased. Networks became modular, consisting of three modules each, and they continued to exhibit SW properties. The highest SW was observed in PRSOL and PRCO pots. In PRCO and CONT networks, the critical function showed a slight response to Chicorum’s growth (with small difference in values between the two samplings), while the value in SOL pots declined remarkably at the legacy sampling. This was also reflected in the natural connectivity of these soils. In contrast, in PRCO, both measurements of robustness were almost identical across the samplings, indicating no effect of Chicorum’s growth on these parameters. In the CONT network, robustness increased after plant growth when the removal concerned only specific bacterial taxa (natural connectivity), while it was unaffected when the removal followed a random process (critical function). The effectiveness parameter measures the contribution of stochastic and determinist processes to network architecture. In almost all networks, the values were near zero, indicating a balance between these two types of processes. A tendency towards stochasticity was observed in the case of SOL pots (0.16).

4 Discussion

4.1 Soil chemical environment and bacterial biomass

In a mesocosm experiment, we analyzed how two invasive species—Conyza bonariensis (annual) and Solanum elaeagnifolium (perennial)—influenced the structure and function of the soil bacterial community. We then examined how the legacy imprinted on these bacterial communities was reshaped by the growth of a native species, the grass Cichorium intybus.

The growth of invaders for one year modified the soil chemical environment by altering the availability of P, K and {\text{NO}}_{\text{3}}^{-} among CO, SOL and CONT pots. Invasive plants significantly decreased the amounts of P, K and {\text{NO}}_{\text{3}}^{-} compared to CONT pots. Differential effects of S. elaeagnifolim on the soil chemical environment have been recorded by Karmezi et al. (2023), who reported an increase in Corg and a decline in Norg in long-term S. elaeagnifolium-invaded sites, specifically in “poor” Quercus stands and “rich” Pinus sites, respectively. In our study, no effect of treatments on organic C and N was recorded, as above-ground plant residues were not incorporated into soil. Karmezi et al. (2023) also reported rapid N mineralization in S. elaeagnifolium soils, supporting the hypothesis that invasive plants facilitate their own growth by maintaining fast nutrient cycles (Kulmatiski et al., 2006). Indeed, the rapid growth of both invasive species likely resulted in a reduced amount of nutrient content in the soil. Although the nutrient status of the soil is generally related to its capacity to sustain microbial biomass, this was not the case in our experiment. Notably, invasion can have positive, negative, or neutral effects on soil microbial properties, but such studies typically compare invaded sites to those with native species. In our experiment, however, the control pots remained unplanted for a year. Consequently, the tendency for higher bacterial biomass in the invasive-treated pots, although non-significant, compared to controls, was likely attributable to root exudates released by the aliens; substances which served as food resources for microbes. The lack of the invader’s species-specific effects on bacterial biomass recorded in this study aligns with Ryan et al. (2023), who found no influence of cultivating six different grassland species on bacterial abundance. According to Torres et al. (2021), allelopathic substances released by invasive species had no effect on bacterial biomass. Indeed, although S. elaeagnifolium is known for allelopathic constituents (Sammani et al., 2013; Balah, 2015; Tsaballa et al., 2015) that exhibit antibacterial and antifungal properties (Bouslamti et al., 2022), the effect of S. elaeagnifolium on bacterial and even fungal biomass (unpublished data) was not negative compared to the controls. The lack of plant species-specific effects and their legacy was also recorded in the soil enzymatic profile, as revealed by the NMDS analysis. However, the effect of plant presence on enzymatic profile was significant in both samplings. Nunes et al. (2025), who examined the legacy effect of two types of vegetation on eucalyptus growth in a microcosm experiment, recorded an effect of plant presence on the potential activities of enzymes involved in soil carbon and sulfur cycles, while the activities of P-mineralization enzymes remained unaffected, a fact that was related to soil nutrient content.

4.2 Bacterial communities’ structure

Bacterial communities, regardless of treatment, were dominated by Firmicutes, with an abundance of 80−90%. The Bacillus genus alone accounted for over 80% of the relative abundance. Members of the Bacillus group are metabolically diverse and capable of forming endospores, which provide high resistance to radiation, desiccation, UV light, heat, and chemicals (Mandic-Mulec et al., 2016). The soil used in this experiment originated from a degraded, dry Mediterranean area subjected to long-term grazing pressure (Kapagianni et al., 2021). Under such adverse conditions, a large proportion of the microbial community likely consisted of Bacillus species that survived as endospores. When this soil was used to establish mesocosms and subjected to regular watering for 1.5 year, the endospores were activated, leading to their dominance in the community. The resulting priority effect exerted by Bacillus spp., which can produce antimicrobial agents and other bioactive compounds, likely shaped the broader composition of the microbial community. Interactions between Bacillus spp. and other microbes, such as arbuscular mycorrhizal fungi or Pseudomonas spp., are well-documented in the literature and depend on factors like soil water and nutrient status (Nikolaidou et al., 2021; Papatheodorou et al., 2021; Du et al., 2022).

The PERMANOVA analysis revealed no significant differences in the composition of bacterial communities among treated soils or between conditioning and legacy samplings, an unexpected finding. In many studies, the impact of plant invasion on soil microbial communities is inconsistent and depends on various factors. For instance, the effect of Phragmites invasion on plant and microbial community composition became stronger as salinity increased (Farrer et al., 2021). In a meta-analysis by Custer and van Diepen (2020), plant invasions were found to have highly heterogenous and limited impacts on microbial α-diversity. Further, Carey et al. (2015) examined the effect of different treatments, among them the effect of exotic plant invasion, and found that bacterial and archaeal communities showed limited variation in composition and diversity across treatments (plant invasion, clipping and N-fertilization). Even in plots subjected to three interacting factors, microbial communities resembled those in restored native grassland. The authors concluded that in Mediterranean grasslands, historical exposure to large seasonal and inter-annual variations in key soil properties, along with prior site history, may select for functionally plastic or largely dormant microbial communities. Among other reasons that could result in similar community composition no matter the invader or the invader’s legacy, in our study a possible explanation could be related with the intrinsic soil characteristics like soil pH which was around 5.5 (Kapagianni et al., 2021). pH has been frequently identified as the key modulating factor for variation in community composition (Rousk et al., 2010). It seems that pH overmarked changes induced by the two invasive species through their root exudates. Also, the overdominance of Bacillus could be related to acidic conditions since specific strains of Bacillus like B. coagulans, B. subtilis and few others show preference for acid soils due to their mechanisms that enable them to regulate their intacellular pH.

Both the deterministic Zipf-Mandelbrot and the neutral ZSM rank abundance models provided a good fit, with no significant differences observed between samplings. Specifically, the beta and gamma parameters of the deterministic model reflected communities with nearly equal participation of most taxa living in homogeneous environments. In such case, we anticipate taxa with narrow niches. Indeed, the mean niche width of the communities was narrow (ranging from one to two), which contrasts with findings by Li et al. (2022) where niche width for soil bacteria was significantly broader. This also contrasts with findings by Xie et al. (2021), where bacterial communities exhibited wide niche breadth due to high metabolic flexibility. Among these specialized communities, the one associated with S. elaeagnifolium at the conditioning sampling demonstrated the highest specificity. At legacy sampling, community specificity decreased, likely due to the introduction of new substances through chicory root exudates. Because community composition did not change across treatments and samplings, resource availability and distribution appear to be critical factors influencing community specificity. This is supported by the linear relationship between mean niche width and the mean activity of PPO and POD at both samplings. Based on these findings, the degree of community specificity seems closely tied to C cycling. Many species of Actinobacteria, Proteobacteria, and Firmicutes are known to be important chitinolytic bacteria in soil (Beier and Bertilsson, 2013). Since niche width was estimated using a niche overlap index, an increase in mean community niche width suggests a rise in taxa co-occurrences, which facilitates the decomposition of recalcitrant organic materials (PPO, POD). This is consistent with the observed increase in POD and PPO activity in the legacy sampling compared to conditioning sampling, indicating that more recalcitrant compounds accumulated in the soil following chicory growth. This accumulation may explain the tendency for reduced bacterial biomass in the legacy sampling.

Both deterministic and stochastic processes influenced the ranking of taxon abundances. Stochastic processes are typically linked with dispersal, ecological drift, and/or genetic drift. In this study, the parameters m (dispersal rate) and I/J (the ratio of immigrants to newborns) exhibited very low values, indicating minimal introduction of new taxa. Therefore, the observed stochasticity in ranking was likely driven by ecological drift (i.e., random fluctuations in the abundance of taxa across soil microsites due to birth and death rates). Alternatively, it could reflect limited dispersal, which may increase local diversification through genetic drift, resulting in transient taxa (Goss-Souza et al., 2020). Although Xu et al. (2022) suggested that generalist communities with wider niche breadths are less influenced by environmental factors and thus exhibit higher stochasticity, while specialized communities are predominantly shaped by deterministic processes, our results deviated from this pattern. Despite the communities being specialized, rank abundance patterns were influenced almost equally by deterministic and stochastic mechanisms. One plausible explanation for this discrepancy is that the indigenous species pool was limited, facilitating stochastic colonization of new microsites. Additionally, the overdominance of Bacillus species likely induced a priority effect, amplifying the role of stochasticity.

Stochastic and deterministic processes can exert varying influences on different aspects of community structure and organization (Stamou and Papatheodorou, 2023). When examining the two components of beta-diversity, the findings diverged from those derived from the abundance models. At the conditioning stage, the ratio β_grad/β_bal revealed that pot-to-pot variation in beta-diversity in CONT soils was predominantly regulated by stochasticity (β-ratio>0.5). Stochasticity often dominates under extreme (He et al., 2022) and/or highly homogeneous conditions. In this case, the unplanted soils, which lacked organic inputs for a year, created extreme conditions for bacterial communities and increased substrate homogeneity due to the absence of vegetation. Conversely, plant growth enhanced soil heterogeneity, both structurally through root system development, and chemically through rhizodeposition. The strong filtering effects of plant roots were evidenced by the predominance of the deterministic regulation of β-diversity in CO soils, and the equal contribution of deterministic and stochastic forces in SOL soils. This finding indicates that invasion shifted the balance between stochasticity and determinism. Notably, the result contrasts with Li et al. (2022), who reported that invaded soils are more stochastic. The discrepancy likely stems from differences in experimental controls: Li et al. (2022) used a vegetated, uninvaded control, whereas our experiment employed unplanted soils as the control.

Biotic heterogeneity, as measured by β-diversity, increased from CONT to SOL at conditioning but converged across treatments after chicory growth. C. intybus played a decisive role in homozenizing β-diversity, aligning with our initial hypothesis. During the legacy sampling, six months after chicory growth, deterministic processes continued to dominate β-diversity regulation in PRCO soils, whereas stochasticity became the main driver in PRSOL soils. This outcome underscores the legacy effects tied to the identity of the invasive species. Specifically, C. bonariensis exerted stronger filtering effects, while S. elaeagnifolium created more stochastic conditions for the bacterial community. Thus, S. elaeagnifolium, beyond its high invasive potential, with its legacy promoted stochastic regulation of β-diversity, fostering conditions conductive to ecological and genetic drift within the bacterial community.

4.3 Co-occurrence patterns

The bacterial networks describing microbial interactions featured relatively few nodes, reflecting the particularly narrow niche widths of the bacterial members. Although the difference in network connectivity between conditioned and legacy pots was slight, the organizational architecture changed significantly. Specifically, networks transitioned from a non-modular to modular structure following the growth of chicory. This suggests that the growth of C. intybus increased habitat heterogeneity, as the number of modules can represent the degree of niche differentiation (Wan et al., 2020). Given that members of each module tend to respond similarly to abiotic factors or act synergistically to access resources (Papatheodorou et al., 2021), we infer that in this more heterogenous environment, module members occupied specific soil microenvironments. The transition from a non-modular to a modular organization was driven by an increase in negative associations, signaling network disruption. Notably, all networks exhibited small-world (SW) properties. The greatest increase in SW after chicory’s growth was observed in PRSOL pots, making this network the least robust to further disturbances, as indicated by critical function and natural connectivity values. By contrast, CO networks demonstrated robustness against both targeted and random disturbances in both samplings. These findings are contrary to our initial hypothesis. The resilience of CO networks aligns with Li et al. (2022), who examined bacterial co-occurrence patterns in soils with and without invasion and found that networks in invaded soils were more robust and contained more keystone species.

Based on effectiveness values, the network architectures across all treatments were governed by a combination of deterministic and stochastic processes. In our networks, deterministic processes mainly operated at the module level, whereas stochastic processes ̶ such as ecological and genetic drift, and limited or unlimited dispersal acted at the level of the entire network. Following the concept that stochasticity increases with rising nutrient availability, while deterministic processes are more associated with low-nutrient conditions (Zhou et al., 2014; Feng et al., 2018) we concluded that both nutrient-rich and nutrient-limited patches co-existed in our mesocosms.

5 Conclusions

Although mesocosm experiments on plant invaders and their legacy effects on soil microbial communities cannot fully replicate field conditions—due to the absence of real-world complexity, including environmental variability, natural plant-root interactions, aboveground litter, unrestricted rooting depth, and limitations in microbial dispersal—they can effectively simulate the early stages of alien plant encroachment onto disturbed soil. As such, they serve as a valuable tool for studying and managing biological invasions. According to our outcomes, invaded soils exhibited higher bacterial biomass compared to unplanted soils, with no species-specific effects or difference in bacterial community composition. However, the legacy of the invasive species did not influence the balance of deterministic and stochastic processes shaping the architecture of soil bacterial networks or to community structure, as indicated by the rank abundance models. Conversely, the invaders’ legacy significantly impacted the robustness of bacterial networks to further disturbance. C. bonariensis’s legacy enhanced network robustness, whereas S. elaeagnifolium created a more disruptive effect, increasing vulnerability under both targeted and random disturbances. Additionally, the type of legacy influenced the components of β-diversity: regulation was primarily deterministic in soils with C. bonariensis and stochastic in soils with S. elaeagnifolium. The stronger filtering of C. bonariensis likely created a more extreme environment, whereas S. elaeagnifolium exerted greater variable selection on bacterial communities, promoting stochastic processes such as ecological and/or genetic drift. Consequently, we suggest that in soils with S. elaeagnifolium legacy, genetic drift may predominate due to limited dispersal rates. We suggested that the inherent characteristics of S. elaeagnifolium necessitate mainly precautionary actions to prevent invasion compared to C. bonariensis, whose invasion may be more straightforward to control. Although our results could inform management actions, they should be interpreted with caution, as this study focused solely on the short-term effects of invasion. Factors such as transient responses, delayed competitive interactions, and seasonal variations were not considered. Additionally, the experiment was conducted in a degraded Mediterranean soil with specific intrinsic properties that likely influenced the outcomes. To derive broader trends regarding the invasion effects of S. elaeagnifolium and C. bonariensis, research should include experiments across soils with diverse characteristics, presenting an important challenge for further investigation.

References

[1] Afzal, M.R., Naz, M., Ashraf, W., Du, D.L., 2023. The legacy of plant invasion: impacts on soil nitrification and management implications. Plants12, 2980.

[2] Allison, S.D., Jastrow, J.D., 2006. Activities of extracellular enzymes in physically isolated fractions of restored grassland soils. Soil Biology and Biochemistry38, 3245–3256.

[3] Andrews S. 2010. FastQC: A quality control tool for high throughput sequence data. available at the website of bioinformatics.babraham.ac.uk.

[4] Balah, M.A.A., 2015. Herbicidal activity of constituents isolated from Solanum elaeagnifolium (Solanaceae). Journal of Crop Protection4, 487–496.

[5] Baselga, A., Orme, C.D.L., 2012. betapart: an R package for the study of beta diversity. Methods in Ecology and Evolution3, 808–812.

[6] Beier, S., Bertilsson, S., 2013. Bacterial chitin degradation-mechanisms and ecophysiological strategies. Frontiers in Microbiology4, 149.

[7] Bezemer, T.M., Lawson, C.S., Hedlund, K., Edwards, .AR., Brook, A.J., Igual, J.M., Mortimer, S.R., van der Putten, W.H., 2006. Plant species and functional group effects on abiotic and microbial soil properties and plant–soil feedback responses in two grasslands. Journal of Ecology 94, 893–904.

[8] Borgatti, S.P., Everett, M.G., Freeman, L.C., 2002. UCINET for Windows: Software for Social Network Analysis. Harvard, MA, USA: Analytic Technologies.

[9] Bouslamti, M., Metouekel, A., Chelouati, T., El Moussaoui, A., Barnossi, A.E., Chebaibi, M., Nafidi, H.A., Salamatullah, A.M., Alzahrani, A., Aboul-Soud, M.A.M., Bourhia, M., Lyoussi, B., Benjelloun, A.S., 2022. Solanum elaeagnifolium var. obtusifolium (Dunal) Dunal: antioxidant, antibacterial, and antifungal activities of polyphenol-rich extracts chemically characterized by use of in vitro and in silico approaches. Molecules27, 8688.

[10] Callaway, R.M., Ridenour, W.M., 2004. Novel weapons: invasive success and the evolution of increased competitive ability. Frontiers in Ecology and the Environment2, 436–443.

[11] Carey, C.J., Beman, J.M., Eviner, V.T., Malmstrom, C.M., Hart, S.C., 2015. Soil microbial community structure is unaltered by plant invasion, vegetation clipping, and nitrogen fertilization in experimental semi-arid grasslands. Frontiers in Microbiology6, 466.

[12] Connell, J.H., 1978. Diversity in tropical rain forests and coral reefs. Science199, 1302–1310.

[13] Corbin, J.D., D’Antonio, C.M., 2012. Gone but Not Forgotten? Invasive plants’ legacies on community and ecosystem properties. Invasive Plant Science and Management5, 117–124.

[14] Custer, G.F., van Diepen, L.T.A., 2020. Plant invasion has limited impact on soil microbial α-diversity: a meta-analysis. Diversity12, 112.

[15] De Coster, W., Rademakers, R., 2023. NanoPack2: population-scale evaluation of long-read sequencing data. Bioinformatics39, btad311.

[16] Du, E.W., Chen, Y.P., Li, Y.H., Sun, Z.X., Gui, F.R., 2022. Rhizospheric bacillus-facilitated effects on the growth and competitive ability of the invasive plant Ageratina adenophora. Frontiers in Plant Science13, 882255.

[17] Duda, J.J., Freeman, D.C., Emlen, J.M., Belnap, J., Kitchen, S.G., Zak, J.C., Sobek, E., Tracy, M., Montante, J., 2003. Differences in native soil ecology associated with invasion of the exotic annual chenopod, Halogeton glomeratus. Biology and Fertility of Soils38, 72–77.

[18] Farrer, E.C., Birnbaum, C., Waryszak, P., Halbrook, S.R., Brady, M.V., Bumby, C.R., Candaele, H., Kulick, N.K., Lee, S.F.H., Schroeder, C.S., Smith, M.K.H., Wilber, W., 2021. Plant and microbial impacts of an invasive species vary across an environmental gradient. Journal of Ecology109, 2163–2176.

[19] Faust, K., Raes, J., 2016. CoNet app: inference of biological association networks using Cytoscape. F1000Research5, 1519.

[20] Feng, Y.Z., Chen, R.R., Stegen, J.C., Guo, Z.Y., Zhang, J.W., Li, Z.P., Lin, X.G., 2018. Two key features influencing community assembly processes at regional scale: initial state and degree of change in environmental conditions. Molecular Ecology27, 5238–5251.

[21] Galil, B.S., 2007. Loss or gain? Invasive aliens and biodiversity in the Mediterranean Sea. Marine Pollution Bulletin55, 314–322.

[22] Goss-Souza, D., Mendes, L.W., Rodrigues, J.L.M., Tsai, S.M., 2020. Ecological processes shaping bulk soil and rhizosphere microbiome assembly in a long-term amazon forest-to-agriculture conversion. Microbial Ecology79, 110–122.

[23] Grman, E., Suding, K.N., 2010. Within-Year soil legacies contribute to strong priority effects of exotics on native California grassland communities. Restoration Ecology18, 664–670.

[24] He, G.X., Peng, T.S., Guo, Y., Wen, S.Z., Ji, L., Luo, Z., 2022. Forest succession improves the complexity of soil microbial interaction and ecological stochasticity of community assembly: evidence from Phoebe bournei-dominated forests in subtropical regions. Frontiers in Microbiology13, 1021258.

[25] Hierro, J.L., Callaway, R.M., 2003. Allelopathy and exotic plant invasion. Plant and Soil256, 29–39.

[26] Humphries, M.D., Gurney, K., 2008. Network ‘small-world-ness’: a quantitative method for determining canonical network equivalence. PLoS One3, e0002051.

[27] Jabot, F., Etienne, R.S., Chave, J., 2008. Reconciling neutral assemblage models and environmental filtering: theory and an empirical test. Oikos 117, 1308–1320.

[28] Kapagianni, P.D., Topalis, I., Gwynn-Jones, D., Menkissoglu-Spiroudi, U., Stamou, G.P., Papatheodorou, E.M., 2021. Effects of plant invaders on rhizosphere microbial attributes depend on plant identity and growth stage. Soil Research59, 225–238.

[29] Karmezi, M., Krigas, N., Papatheodorou, E.M., Argyropoulou, M.D., 2023. The invasion of alien populations of solanum elaeagnifolium in two mediterranean habitats modifies the soil communities in different ways. Plants12, 2193.

[30] Klein, B., Swain, A., Byrum, T., Scarpino, S.V., Fagan, W.F., 2022. Exploring noise, degeneracy and determinism in biological networks with the einet package. Methods in Ecology and Evolution13, 799–804.

[31] Knolmajer, B., Jócsák, I., Taller, J., Keszthelyi, S., Kazinczi, G., 2024. Common ragweed—Ambrosia artemisiifolia L. : a review with special regards to the latest results in biology and ecology. Agronomy14, 497.

[32] Kulmatiski, A., Beard, K.H., Stark, J.M., 2006. Soil history as a primary control on plant invasion in abandoned agricultural fields. Journal of Applied Ecology43, 868–876.

[33] Kulmatiski, A., Heavilin, J., Beard, K.H., 2011. Testing predictions of a three-species plant–soil feedback model. Journal of Ecology99, 542–550.

[34] Kumar, R.P., Singh, J.S., 2020. Invasive alien plant species: their impact on environment, ecosystem services and human health. Ecological Indicators111, 106020.

[35] Lanfear, R., Schalamun, M., Kainer, D., Wang, W., Schwessinger, Β., 2019. MinIONQC: fast and simple quality control for MinION sequencing data. Bioinformatics35, 523–525.

[36] Langenheder, S., Lindström, E.S., 2019. Factors influencing aquatic and terrestrial bacterial community assembly. Environmental Microbiology Reports11, 306–315.

[37] Li, C.C., Bo, H.Z., Song, B.Z., Chen, X.C., Cao, Q.Q., Yang, R.R., Ji, S.P., Wang, L.F., Liu, J., 2022. Reshaping of the soil microbiome by the expansion of invasive plants: shifts in structure, diversity, co-occurrence, niche breadth, and assembly processes. Plant Soil477, 629–646.

[38] Liendo, D., Campos, J.A., Gandarillas, A., 2023. Cortaderia selloana, an example of aggressive invaders that affect human health, yet to be included in binding international invasive catalogues. NeoBiota89, 229–237.

[39] Liu, J., Zhou, M.X., Wang, S., Liu, P.H., 2017. A comparative study of network robustness measures. Frontiers of Computer Science11, 568–584.

[40] Lu, J., Breitwieser, F.P., Thielen, P., Salzberg, S.L., 2017. Bracken: estimating species abundance in metagenomics data. PeerJ Computer Science3, e104.

[41] Luan, L., Liang, C., Chen, L.J., Wang, H.T., Xu, Q.S., Jiang, Y.J., Sun, B., 2020. Coupling bacterial community assembly to microbial metabolism across soil profiles. mSystems5, e00298–20.

[42] Luo, X., Han, S., Fu, X., Li, X., Wang, L., Peng, S., Chen, W., Huang, Q., 2019. The microbial network in naturally fertile paddy soil possibly facilitates functional recruitment in the rice mature stage. Applied Soil Ecology 135, 174–181.

[43] Mandic-Mulec, I., Stefanic, P., van Elsas, J.D., 2016. Ecology of Bacillaceae. In: Driks, A., Eichenberger, P., eds. The Bacterial Spore: From Molecules to Systems. Washington: ASM Press59–85.

[44] McKinley, V.L., Peacock, A.D., White, D.C., 2005. Microbial community PLFA and PHB responses to ecosystem restoration in tallgrass prairie soils. Soil Biology and Biochemistry37, 1946–1958.

[45] Mooney, H.A., Hobbs, R.J., 2000. Invasive Species in a Changing World. 2nd ed. Washington, DC: Island Press, 457.

[46] Nikolaidou, C., Monokrousos, N., Kapagianni, P.D., Orfanoudakis, M., Dermitzoglou, T., Papatheodorou, E.M., 2021. The effect of Rhizophagus irregularis, Bacillus subtilis and water regime on the Plant–Microbial soil system: the case of Lactuca sativa. Agronomy11, 2183.

[47] Nunes, E.A., Cassiano, G.H., da Silveira, A.P.D., De Andrade, S.A.L., 2025. Soil legacies left by a 20-year eucalypt plantation and a secondary vegetation covers on young eucalypt plants and plant-soil feedback. Biology and Fertility of Soils61, 187–210.

[48] Oksanen, J., Simpson, G.L., Blanchet, F.G., Kindt, R., Legendre, P., Minchin, P.R., O'Hara, R.B., Solymos, P., Stevens, M.H.H., Szoecs, E., Wagner, H., Barbour, M., Bedward, M., Bolker, B., Borcard D., Carvalho, G., Chirico, M., De Caceres, M., Durand, S., Evangelista, H.B.A., FitzJohn, R., Friendly, M., Furneaux, B., Hannigan, G., Hill, M.O., Lahti, L., McGlinn, D., Ouellette, M.H., Cunha, E.R., Smith, T., Stier, A., Ter Braak, C.J.F., Weedon, J., Borman, T., 2017. Vegan: community ecology package. R package version 2.4-3 [Online]. available at the website of cran.r-project.org/web/packages/vegan/vegan.pdf.

[49] Papatheodorou, E.M., Monokrousos, N., Angelina, E., Stamou, G.P., 2021. Robustness of rhizosphere microbial communities of L. sativa originated from soils of different legacy after inoculation with Plant Growth Promoting Rhizobacteria. Applied Soil Ecology167, 104028.

[50] Pérez-Valera, E., Goberna, M., Faust, K., Raes, J., García, C., Verdú, M., 2017. Fire modifies the phylogenetic structure of soil bacterial cooccurrence networks. Environmental Microbiology19, 317–327.

[51] Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., Peplies, J., Glöckner, F.O., 2013. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Research41, D590–D596.

[52] Rillig, M.C., Mummey, D.L., Ramsey, P.W., Klironomos, J.N., Gannon, J.E., 2006. Phylogeny of arbuscular mycorrhizal fungi predicts community composition of symbiosis-associated bacteria. FEMS Microbiology Ecology57, 389–395.

[53] Rousk, J., Bååth, E., Brookes, P.C., Lauber, C.L., Lozupone, C., Caporaso, J.G., Knight, R., Fierer, N., 2010. Soil bacterial and fungal communities across a pH gradient in an arable soil. The ISME Journal4, 1340–1351.

[54] Ryan, K.B., De Menezes, A., Finn, J.A., Brennan, F.P., 2023. Plant species and soil depth differentially affect microbial diversity and function in grasslands. Journal of Sustainable Agriculture and Environment2, 397–411.

[55] Sammani, A., Shammaa, E., Chehna, F., 2013. Qualitative and quantitative steroidal alkaloids of Solanum species distributed widely in Syria by TLC and HPLC. International Journal of Pharmaceutical Sciences Review and Research23, 23–27.

[56] Shannon, P., Markiel, A., Ozier, O., Baliga, N.S., Wang, J.T., Ramage, D., Amin, N., Schwikowski, B., Ideker, T., 2003. Cytoscape: a software environment for integrated models of Biomolecular Interaction Networks. Genome Research13, 2498–2504.

[57] Simberloff, D., Martin, J.L., Genovesi, P., Maris, V., Wardle, D.A., Aronson, J., Courchamp, F., Galil, B., García-Berthou, E., Pascal, M., Pyšek, P., Sousa, R., Tabacchi, E., Vilà, M., 2013. Impacts of biological invasions: what’s what and the way forward. Trends in Ecology & Evolution28, 58–66.

[58] Sinsabaugh, R.L., Reynolds, H., Long, T.M., 2000. Rapid assay for amidohydrolase (urease) activity in environmental samples. Soil Biology and Biochemistry32, 2095–2097.

[59] Stamou, G.P., Panagos, P., Papatheodorou, E.M., 2024. Connections between soil microbes, land use and European climate: insights for management practices. Journal of Environmental Management360, 121180.

[60] Stamou, G.P., Papatheodorou, E.M., 2023. Deterministic versus stochastic control in β-diversity, abundance and co-occurrence patterns of a soil nematode assemblage living in a Mediterranean soil. Applied Soil Ecology188, 104879.

[61] Tamburello, N., Litt, M.A., 2023. Multiple impacts of invasive species on species at risk: a case study in British Columbia, Canada. FACETS8, 1–13.

[62] Torres, N., Herrera, I., Fajardo, L., Bustamante, R.O., 2021. Meta-analysis of the impact of plant invasions on soil microbial communities. BMC Ecology and Evolution21, 172.

[63] Tsaballa, A., Nikolaidis, A., Trikka, F., Ignea, C., Kampranis, S.C., Makris, A.M., Argiriou, A., 2015. Use of the de novo transcriptome analysis of silver-leaf nightshade (Solanum elaeagnifolium) to identify gene expression changes associated with wounding and terpene biosynthesis. BMC Genomics16, 504.

[64] Turbelin, A.J., Cuthbert, R.N., Essl, F., Haubrock, P.J., Ricciardi, A., Courchamp, F., 2023. Biological invasions are as costly as natural hazards. Perspectives in Ecology and Conservation21, 143–150.

[65] Uludag, A., Gbehounou, G., Kashefi, J., Bouhache, M., Bon, M.C., Bell, C., Lagopodi, A.L., 2016. Review of the current situation for Solanum elaeagnifolium in the Mediterranean Basin. EPPO Bulletin46, 139–147.

[66] van der Heijden, M.G.A., Klironomos, J.N., Ursic, M., Moutoglis, P., Streitwolf-Engel, R., Boller, T., Wiemken, A., Sanders, I.R., 1998. Mycorrhizal fungal diversity determines plant biodiversity, ecosystem variability and productivity. Nature396, 69–72.

[67] Van Der Putten, W.H., 2003. Plant defense belowground and spatiotemporal processes in natural vegetation. Ecology84, 2269–2280.

[68] Vitousek, P.M., Walker, L.R., 1989. Biological iInvasion by Myrica faya in Hawai'i: plant demography, nitrogen fixation, ecosystem effects. Ecological Monographs59, 247–265.

[69] Wan, X.L., Gao, Q., Zhao, J.S., Feng, J.J., van Nostrand, J.D., Yang, Y.F., Zhou, J.Z., 2020. Biogeographic patterns of microbial association networks in paddy soil within Eastern China. Soil Biology and Biochemistry142, 107696.

[70] White, D., Stair, J., Ringelberg, D., 1996. Quantitative comparisons of in situ microbial biodiversity by signature biomarker analysis. Journal of Industrial Microbiology & Biotechnology17, 185–196.

[71] Wood, D.E., Lu, J., Langmead, B., 2019. Improved metagenomic analysis with Kraken 2. Genome Biology20, 257.

[72] Xie, J., Wang, X.Q., Xu, J.W., Xie, H.W., Cai, Y.H., Liu, Y.Z., Ding, X., 2021. Strategies and structure feature of the aboveground and belowground microbial community respond to drought in wild rice (Oryza longistaminata). Rice14, 79.

[73] Xiong, D., Wei, C.Z., Wang, X.G., Lü, X.T., Fang, S., Li, Y.B., Wang, X.B., Liang, W.J., Han, X.G., Bezemer, T.M., Li, Q., 2021. Spatial patterns and ecological drivers of soil nematode β-diversity in natural grasslands vary among vegetation types and trophic position. Journal of Animal Ecology90, 1367–1378.

[74] Xu, Q.C., Vandenkoornhuyse, P., Li, L., Guo, J.J., Zhu, C., Guo, S.W., Ling, N., Shen, Q.R., 2022. Microbial generalists and specialists differently contribute to the community diversity in farmland soils. Journal of Advanced Research40, 17–27.

[75] Zhang, C.P., Xu, J., Liu, X.G., Dong, F.S., Kong, Z.Q., Sheng, Y., Zheng, Y.Q., 2010. Impact of imazethapyr on the microbial community structure in agricultural soils. Chemosphere81, 800–806.

[76] Zhou, J.Z., Deng, Y., Zhang, P., Xue, K., Liang, Y.T., Van Nostrand, J.D., Yang, Y.F., He, Z.L., Wu, L.Y., Stahl, D.A., Hazen, T.C., Tiedje, J.M., Arkin, A.P., 2014. Stochasticity, succession, and environmental perturbations in a fluidic ecosystem. Proceedings of the National Academy of Sciences of the United States of America111, E836–E845.

[77] Zhou, Z.H., Zheng, M.H., Xia, J.Y., Wang, C.K., 2022. Nitrogen addition promotes soil microbial beta diversity and the stochastic assembly. Science of the Total Environment806, 150569.

RIGHTS & PERMISSIONS

The Author(s) 2025. This article is published with open access at link.springer.com and journal.hep.com.cn