The S. alterniflora invasion zone was located along the Jiangsu coast on the western edge of the Yellow Sea. The 954-km coastline lies at the center of the East Asian-Australasian Flyway and encompasses four major migratory bird staging areas, providing critical breeding, stopover, and wintering habitats for tens of millions of birds (Ma et al., 2023; Melville et al., 2016; Wang et al., 2022; Zhao et al., 2015). It also contains the world’s largest intertidal wetland (Barter and Riegen, 2004). Jiangsu has a temperate monsoon climate, with annual rainfall ranging from 800 to 1,200 mm (Qin et al., 2022) and average temperatures between 13 and 16 °C (Qin et al., 2022). Coastal soils are mostly saline-alkaline and alluvial, commonly occurring in tidal flats (Liu et al., 2015). Physical S. alterniflora removal began in August 2023 and was nearly complete by February 2024. Technicians first cut aboveground biomass using brush cutters and tractors equipped with mowing blades. Excavators and rotary tillers then ploughed and tilled the root zone to a depth of approximately 2 m (Supplementary Figure S3). Removal covered continuous patches averaging 5–10 ha per site. All treatments finished at least two months before subsequent soil, macrobenthos, and bird surveys (Qi et al., 2024). Within S. alterniflora patches, dense rhizome–root mats and accumulated litter increased sediment consolidation and surface shear strength. These conditions created firm, load-bearing surfaces that allowed safe foot access during surveys (Supplementary Figure S3). We defined two conditions: RA (removal area) and ARA (absence of removal) (Qi et al., 2024). Sites were selected from Lianyungang to Yangtze River estuary. Soil, macrobenthos, and bird surveys were conducted in ARA (April–July 2023) and repeated in RA (April–July 2024) to assess post-removal changes (Figure 1). All field surveys were conducted within the same seasonal survey window, and sampling months were kept consistent across years to ensure temporal alignment among soil, macrobenthic, and bird measurements and to minimize seasonal variability in bird activity and associated ecological processes (Bibby et al., 1998).

We collected soil samples in both RA and ARA at the start, middle, and end points of each 1-km transect, matching bird and macrobenthos survey locations. Transects were established perpendicular to the shoreline, each 1 km in length, with sampling points fixed at 0 m, 500 m, and 1,000 m. At each point, three replicate soil cores were collected (n = 9 per transect). Sampling was conducted in April, May, and July of each year to capture seasonal variability. In RA, S. alterniflora was removed by cutting aboveground biomass and tilling the soil to a depth of approximately 2 m.

At each site, three replicate soil samples (0–20 cm) were collected and dried using a Pilot10-15S vacuum freeze dryer (BoyiKang, Beijing, China). The pH was measured in a 1:5 soil–water mixture using a PHSJ-6L pH meter (INESA, Shanghai, China). Samples were obtained via the ring knife method, and total soluble salts (TS) were extracted using a 1:5 leaching method; the leachate was dried and weighed. For each replicate, sediment depth was standardized to 20 cm using a ring-knife sampler (100 cm³, model TZ-100, Tianjin, China). Total organic carbon (OC) was measured by potassium dichromate oxidation and ferrous sulfate titration (Lyu et al., 2023b). TN by the Kjeldahl method (Xu et al., 2016). AP by the ammonium molybdate colorimetric method (Moonrungsee et al., 2015), and AK via ICP-AES (Krogstad and Zivanovic, 2022).

Macrobenthic organism surveys were conducted along 1-km bird line transects at each site (Bibby et al., 1998; Seys et al., 2002) (Figure 1). Sampling was conducted during the breeding season to coincide with peak biological activity. At each transect, three sampling points were selected—at the start, middle, and end—to ensure spatial representation. Macrobenthic organisms were collected using a 20 × 30 cm quadrat (volume: 12,000 cm³) (Mereta et al., 2013). Sediments were sieved through a 0.5-mm mesh. At each of the three sampling points per transect (0 m, 500 m, 1,000 m), three replicate quadrats were taken, resulting in nine samples per transect. Sediment depth was standardized to 20 cm. All specimens were preserved in 4% buffered formalin solution, with exposure time limited to 48 h before transfer to 70% ethanol for long-term storage (Clifford, 1991; Mereta et al., 2013). In the laboratory, organisms were identified under a microscope, and their species richness was recorded (Clifford, 1991; Mereta et al., 2013). Macrobenthic organisms were identified to the lowest possible taxonomic level using standard taxonomic keys and regional monographs, primarily following Checklist of Marine Biota of China Seas (Liu, 2008), supplemented by widely used identification guides for polychaetes and mollusks (Fauchald, 1977; Gosling, 2008).

Bird surveys were conducted during the local peak breeding season. We established ten transects, each 1 km in length. Each transect was surveyed four times: twice during high tide and twice during low tide (Bibby et al., 1998). Transects crossed S. alterniflora stands and adjacent intertidal flats (Bibby et al., 1998; Jackson et al., 2021; Cheng et al., 2025). Observers maintained a distance of 40–120 m from the birds. To avoid disturbance, observers did not walk on dikes directly adjacent to birds but used transects farther away (Supplementary Figure S3). Observers used binoculars (Swarovski EL 10×42) to count individuals and spotting scopes (Swarovski ATX 25–60×85) to confirm species identity (Supplementary Figure S3).

Each survey lasted 45–120 minutes depending on bird abundance, but no survey exceeded 2 hours to avoid double-counting. One round of surveys covered all transects across the study sites. Four rounds were conducted from April to July. In each round, every transect was surveyed once, and no transect was surveyed more than once per day. Observers used handheld GPS units (Garmin GPSMAP 64st) to record transect locations. They maintained a steady pace of 0.5–1 km/h to ensure complete coverage and minimize missed detections. Bird surveys were conducted at comparable tidal states across sites, targeting similar stages of high or low tide. The reported morning (7:00-11:30 AM) and afternoon (3:00-6:30 PM) periods represent typical time windows during which these tidal conditions occurred, rather than fixed survey times (Sebastián-Gonzáez et al., 2010; Speckman et al., 2003). These times ensured detection of species that were more active during specific tidal phases (Sebastián-Gonzáez et al., 2010; Speckman et al., 2003).

We obtained the tidal flat area for both ARA and RA from the Tidal Flats Map of China (TFMC) (Chen et al., 2025). The map was generated using 30-m Sentinel-2 imagery from 2023 (ARA) and 2024 (RA). Minimum water extent was extracted using total water deficit, and maximum water extent was mapped using the modified normalized difference water index (Chen et al., 2025). For each site, we clipped TFMC tidal-flat polygons using GPS-defined transect buffers. We manually removed aquaculture ponds and artificial structures (e.g., dikes and seawalls) through visual interpretation of high-resolution remote sensing imagery to avoid misclassification as tidal flats. We excluded vegetated zones dominated by S.alterniflora using NDVI thresholding (NDVI > 0.3) combined with species-specific vegetation products. The difference between these two extents defined the tidal flat area. This mapping process was implemented in Google Earth Engine. Vegetated and aquaculture zones were excluded. TFMC product showed high classification accuracy, with an F1 score greater than 0.97.

We used the Wilcoxon rank-sum test to compare soil properties, macrobenthic communities, and bird communities between ARA and RA, due to violations of normality and homogeneity (confirmed by Shapiro-Wilk and Levene’s tests). Soil variables included pH, TS, OC, TN, AP, and AK; community metrics included macrobenthic and bird richness, macrobenthic abundance, and bird individual numbers. All analyses were conducted in R using the wilcox.test() function, with significance set at p < 0.05 and effect sizes reported.

Principal component analysis (PCA) was used to assess variation in soil properties within S. alterniflora-dominated habitats. Soil variables included pH, TS, OC, TN, AP, and AK. Euclidean biplots based on the first two principal components were generated using the FactoMineR package (Lê et al., 2008). Bird and macrobenthic community structures were assessed using non-metric multidimensional scaling (NMDS) (Agarwal et al., 2007), based on square-root-transformed species abundance and Bray–Curtis dissimilarity (Clarke and Warwick, 2001). NMDS was conducted separately for each group using the metaMDS() function in the vegan package (Oksanen, 2022), with 999 permutations and two dimensions (k = 2). Final stress values were < 0.3, indicating acceptable ordination fit. Group-level differences (ARA vs. RA) were tested using PERMANOVA (adonis()), and multivariate dispersion was assessed using betadisper(). Stress values and ordination plots were used to evaluate model fit and visualize community structure.

Phylogenetic trees were downloaded from BirdTree.org (https://birdtree.org/) using the “Hackett All Species: a set of 10,000 trees with 9,993 OTUs each” option (Jetz et al., 2012), and were subsequently pruned to include only the bird species present in our dataset. From these, 1,000 pseudoposterior phylogenetic trees were randomly selected, and a maximum clade credibility tree was constructed using TreeAnnotator in BEAST, with node heights set to the mean of ancestral estimates (Drummond and Rambaut, 2007; Ricklefs and Jønsson, 2014). Phylogenetic similarity at the species level was calculated using the variance–covariance matrix derived from the pruned tree. Nested random effects were included in the model by nesting individual identifiers within species names.

We used a Bayesian Phylogenetic Structural Equation Model (PSEM), implemented via the brms package, to test causal relationships among S. alterniflora types (ARA and RA), environmental variables (e.g., tidal flat area and soil physicochemical properties represented by PC1 and PC2), macrobenthic richness and abundance, and bird richness and individual number (Bürkner, 2017; Lefcheck, 2016). Due to collinearity in tidal flat area between ARA and RA, separate models were constructed for each habitat type. To reduce redundancy among soil variables while capturing key environmental gradients, the first two principal components were extracted (Waldemar, 1993). PC1 and PC2 served as proxies for soil properties, with PC1 reflecting variation in soil pH, OC, and TN, while PC2 reflected variation in TS and AP (Grace, 2006; Li et al., 2022). A phylogenetic covariance matrix derived from the pruned tree was included as a random effect to account for shared evolutionary history among species. The PSEM was specified a priori to include both direct and indirect pathways linking soil properties, macrobenthic communities, and bird communities. Specifically, bird richness and bird individual number were modeled as functions of PC1 and PC2 (direct soil effects) as well as macrobenthic richness and macrobenthic abundance (indirect effects). Path significance was assessed using posterior estimates and their 90% credible intervals (Li et al., 2022; Vehtari et al., 2018). Model specification followed recent methodological developments in Bayesian phylogenetic structural equation modeling, which demonstrate that complex causal structures with latent variables can be estimated under moderate sample sizes when model dimensionality is constrained (von Hardenberg and Gonzalez-Voyer, 2025). Accordingly, dimensionality reduction and habitat-specific model fitting were applied to limit the number of free parameters relative to the available sample size. The model estimated interactions among soil properties (PC1 and PC2), macrobenthic metrics (richness and abundance), and bird metrics (richness and individual number). Four Markov chains were run for 1,000 iterations, with the first 300 discarded as burn-in. All explanatory variables were standardized (mean = 0, variance = 1), and no strong multicollinearity was detected (VIF < 2) (Dormann et al., 2013; O’brien, 2007). Default priors were used. Convergence was assessed via the Rubin–Gelman statistic (R̂ < 1.1), and model performance was evaluated using leave-one-out cross-validation via the loo package (Vehtari et al., 2018). This procedure provided a formal assessment of model adequacy at the specified level of model complexity. Explanatory variables were considered significant if their 90% credible intervals excluded zero, with effect size interpreted as an indicator of relative importance.

All six soil variables differed significantly between ARA and RA (all p < 0.001) (Figures 2a-f; Supplementary Figures S2a-f). PH, OC, TS, and AP were significantly higher in RA (p < 0.001) (Figures 2a-d; Supplementary Figures S2a-d), with nearly a threefold increase. In contrast, TN and AK were significantly lower in RA (p < 0.001) (Figures 2e, f; Supplementary Figures S2e, f). PCA showed a significant difference between ARA and RA (p < 0.001) (Supplementary Figure S1). The first two principal components explained 83.20% of the total variance (PC1 = 50.70%, PC2 = 32.50%). RA samples clustered on the right side of PC1. These were associated with higher pH, OC, and AP. ARA samples clustered on the left and were linked to higher TN and AK Supplementary Figure S1).

We observed 48 macrobenthic species from 29 orders and 31 families, with a total of 913 individuals across all sites (Supplementary Table S1). Macrobenthic richness and abundance were significantly higher in ARA than those in RA (p < 0.001) (Figures 3a, b, d, e). NMDS ordination showed no significant difference in macrobenthic community composition between ARA and RA (p = 0.83) (Figure 3c).The polychaete Perinereis aibuhitensis (n = 408) dominated the assemblage, accounting for nearly half of all individuals. High abundances were also observed for the bivalve Glauconome primeana (n = 204) and the crab Sesarma plicata (n = 60), both largely restricted to ARA sites. Species occurrences differed between treatments: Perinereis aibuhitensis (n = 408), Glauconome primeana (n = 204), and Helice sheni (n = 39) were dominant in ARA, whereas Assiminea latericea (n = 34), Cyclina sinensis (n = 5), and Acetes chinensis (n = 5) were more frequently observed in RA.

We observed 50 bird species from 7 families, with a total of 79,320 individuals across all sites (Supplementary Table S2). Bird richness (p < 0.001) (Figures 4a–d) and abundance (p < 0.05) differed significantly between ARA and RA (Figures 4b–e). NMDS ordination showed no significant difference in bird community composition between ARA and RA (p = 0.19) (Figure 4c). The bird community was numerically dominated by the Pied Avocet (Recurvirostra avosetta) (n = 41,368), followed by the Far Eastern Curlew (Numenius madagascariensis) (n = 12,762) and the Little Egret (Egretta garzetta) (n = 8,649), which were recorded in both ARA and RA with comparable abundances. Most species were observed under both ARA and RA conditions. A subset of species, including the Pacific Golden-Plover (Pluvialis fulva) (n = 22), Long-billed Plover (Charadrius placidus) (n = 235), and Ruddy Turnstone (Arenaria interpres) (n = 6), was restricted to RA, while others, such as the Eurasian Spoonbill (Platalea leucorodia) (n = 10), Black-crowned Night Heron (Nycticorax nycticorax) (n = 45), and Chinese Egret (Egretta eulophotes) (n = 1), occurred only in ARA.

SEM for ARA revealed distinct patterns (Figure 5a; Supplementary Table S3). Tidal flat area had a significant negative effect on PCA2 (95% CI [−0.42, −0.19]) (Figure 5e). PCA1 had a significant positive effect on macrobenthic richness (95% CI [0.21, 0.44]) and macrobenthic abundance (95% CI [0.19, 0.42]) (Figures 5b, c). PCA2 also had a significant positive effect on macrobenthic abundance (95% CI [0.02, 0.27]) (Figure 5d). However, both macrobenthic richness and abundance had positive but non-significant effects on bird individual number and richness.

SEM for RA showed a different structure (Figure 6a; Supplementary Table S3). Tidal flat area had a significant positive effect on PCA1 (95% CI [0.09, 0.35]) (Figure 6e). PCA2 had significant negative effects on both macrobenthic richness (95% CI [−0.35, −0.11]) and abundance (95% CI [−0.36, −0.12]) (Figures 6b, c). Macrobenthic abundance had a significant positive effect on bird individual number and richness (95% CI [0.09, 0.29]) (Figure 6d). Model performance diagnostics indicated adequate explanatory power and stable predictive performance for both ARA and RA SEMs, with moderate conditional R² values and predominantly acceptable Pareto-k diagnostics from leave-one-out cross-validation (Supplementary Tables S4, S5).

Physical S. alterniflora removal significantly altered soil structure, with a significant pH increase after removal. This finding supports previous reports showing that S. alterniflora roots acidify the soil by secreting organic acids (Li et al., 2020). Root-mediated acidification has been shown to buffer sediment pH during invasion (Li et al., 2024), and its cessation following removal can therefore relax biological control over soil acidity. The rise in pH may result from the cessation of acid root exudation and increased intrusion of alkaline seawater during tidal flushing. Reduced microbial activity and base cation accumulation in the absence of vegetation may also contribute. TS increased significantly, consistent with Yang and Guo (2018), who found that S. alterniflora promotes desalinization through salt exclusion and evapotranspiration. This process implies that intact vegetation regulates salt distribution within surface soils (Manousaki and Kalogerakis, 2011), whereas removal disrupts this regulation and facilitates salt accumulation through evaporation and reduced ion uptake. After removal, vegetation loss likely enhanced surface evaporation and disrupted ion uptake, promoting salt accumulation in surface soils. Additionally, TS elevation may reflect upward movement of saline groundwater under diminished root barriers.

Additionally, we observed a significant increase in OC following S. alterniflora removal. This contrasts with previous studies indicating that S. alterniflora contributes to soil carbon inputs through litterfall and rhizodeposition (Yang and Guo, 2018). Such plant-derived inputs play a central role in maintaining soil carbon pools, and their removal can increase carbon vulnerability under disturbed conditions (Jackson et al., 2017). The increase in OC may be associated with short-term retention of belowground organic residues and altered sediment conditions following physical disturbance, rather than sustained carbon accumulation. Additionally, post-removal changes in microbial activity and sediment aeration may regulate the subsequent redistribution and turnover of OC. These findings are consistent with evidence showing that vegetation removal promotes nitrogen loss via leaching and ammonia volatilization (Cheng et al., 2025; Weidlich et al., 2020). TN declined, likely due to reduced root uptake, disrupted nitrification-denitrification, and microbial competition under stress. In contrast, AP increased significantly post-removal, possibly driven by phosphorus release from decomposing biomass, reduced microbial immobilization, and desorption under alkaline conditions, where higher pH weakens phosphorus fixation. Alkaline shifts have been shown to reduce phosphorus sorption capacity, thereby increasing short-term bioavailability following physical disturbance. AK declined sharply, consistent with previous findings under physical disturbance (Lv et al., 2023). Elevated salinity and increased Na⁺ competition can constrain potassium retention at exchange sites, particularly in mechanically disturbed soils. In high-salinity soils, potassium mobility can also be suppressed by elevated ionic strength.

Our results showed that macrobenthic species richness declined significantly after physically removing S. alterniflora. This finding is consistent with Lu et al. (2022), who observed minimal improvements in macrobenthic diversity following removal interventions. Previous studies indicate that mechanical removal can homogenize sediment structure and reduce fine-scale habitat complexity (Lu et al., 2022), thereby constraining recolonization by taxa with specific substrate or oxygen requirements. In our study, physical removal flattened sediment layers and increased surface exposure, reducing vertical gradients in grain size and redox conditions that support burrowing and oxygen-dependent taxa (Feng et al., 2014; Neira et al., 2006). These effects reflect the combined legacy of long-term S. alterniflora dominance and additional disturbance caused by post-removal tilling (Lyu et al., 2023b; Neira et al., 2006). Prolonged compaction limits vertical permeability and oxygen diffusion, forming stable anoxic layers that reduce habitat suitability for burrowing and oxygen-dependent invertebrates. Post-removal mechanical tilling likely disrupted the substrate, delaying recolonization by sensitive taxa such as Assiminea latericea and Bullacta caurina. Additionally, sediment exposure may suppress microalgae growth, reducing food availability for macrobenthic herbivores.

We also found that macrobenthic abundance decreased after removal. Although NMDS indicated no significant difference in macrobenthic community structure between ARA and RA, the decreases in richness and abundance suggest that macrobenthic communities are still in an early recovery stage following disturbance. This pattern reflects disturbance-driven population suppression rather than rapid species turnover (Feng et al., 2014; Lu et al., 2022). Species persistence across treatments indicates limited turnover, while reduced abundances reflect constraints on population performance under post-removal environmental conditions (Feng et al., 2014; Lu et al., 2022). Comparable patterns have been reported in disturbed intertidal systems where physical alteration reduces resource availability and physiological tolerance but does not exclude taxa from the regional species pool (Feng et al., 2014). However, our results suggest that bottom-up effects on birds remained weak in ARA, which may reflect limited energy transfer despite enhanced macrobenthic productivity. This result contrasts with Lyu et al. (2023b), who reported increases in macrobenthic density after S. alterniflora eradication. One possible explanation is that disturbance removed surface detritus and disrupted food availability. These effects may impact detritivores such as mollusks and polychaetes (Lyu et al., 2023b; Neira et al., 2006). Changes in salinity and pH may also cause stress and reduce abundance (Feng et al., 2014). Our findings support earlier research showing that short-term mechanical treatments can lower macrofaunal biomass and delay recovery (Cui et al., 2016). Site-specific comparisons showed that most locations exhibited reduced richness and abundance after removal.

We provide evidence that physical S. alterniflora removal was associated with changes in bird communities, mediated by shifts in both abiotic drivers and biotic trophic pathways. Structural equation models revealed distinct mechanisms in ARA and RA, suggesting a temporal shift in ecological control. In ARA, OC, TN, and TS had significant positive effects on macrobenthic richness and abundance, consistent with previous studies showing that organic carbon and nitrogen enrichment promote macrobenthic productivity, while moderate salinity enhances faunal tolerance (Beukema et al., 1993; Patten et al., 2017). The positive effect of pH on macrobenthic abundance also supports earlier findings suggesting that acid secretion by S. alterniflora creates favorable microhabitats for invertebrates (Feugere et al., 2021). However, macrobenthic variables had no significant effect on bird communities in ARA, suggesting that abiotic factors were the dominant drivers. This pattern may reflect early invasion dynamics, where habitat structure and openness, rather than trophic interactions, influenced bird use (Jackson et al., 2021). Nutrient enrichment likely benefited lower-level consumers, but limited energy transfer weakened bottom-up effects on birds during early S. alterniflora invasion.

In contrast, we found that in RA, AK exerted a significant negative effect on macrobenthic richness and abundance. This aligns with literature showing that excess potassium may cause ionic imbalance, hinder microbial stability, or limit detritivore recovery after soil disturbance (Mateos-Naranjo et al., 2012). More importantly, macrobenthic abundance increased bird abundance, supporting a bottom-up regulatory mechanism. This finding echoes Ma et al. (2023), who emphasized the importance of trophic connectivity in wetland ecosystems following macrobenthic change (Beukema et al., 1993; Ma et al., 2023; Patten and O’Casey, 2007). The observed trophic pathway indicates that bird recolonization was driven by macrobenthic prey availability rather than direct nutrient conditions, reflecting a full cascade from soil chemistry to higher trophic levels.

At the same time, the increase in bird individual numbers following removal should not be interpreted as evidence of immediate improvement in habitat quality or food availability (Mott et al., 2023). Instead, this pattern likely reflects enhanced habitat openness and accessibility after vegetation removal, which rapidly exposes tidal flats and reduces structural barriers to movement (Dai et al., 2025; Mott et al., 2023). Under such conditions, higher bird abundance may result from aggregation or redistribution responses, even when macrobenthic prey communities remain in an early recovery stage. We predicted that nutrient shifts caused by vegetation removal would reorganize trophic structure. Earlier studies indicated that dense S. alterniflora stands limited bird activity and created ecological traps (Lyu et al., 2023b; Patten and O’Casey, 2007; Tang et al., 2021). The legacy of habitat degradation likely delayed macrobenthic recovery and trophic restructuring. Our results support this, showing a shift in bird responses from abiotic control in ARA to biotic mediation in RA. This suggests that trophic recovery lagged behind soil reorganization, with macrobenthic restructuring ultimately driving avian diversity. The re-establishment of trophic linkages reflects an ecosystem in transition, where macrobenthic organisms mediate energy flow from soils to birds. However, this study compared absence of removal areas and removal areas across different years and did not include a simultaneous untreated control. As a result, the observed patterns should be interpreted as associations consistent with removal effects rather than definitive causal relationships. Although surveys were conducted within the same breeding-season window, interannual variability in climatic and hydrological conditions, including differences in precipitation, tidal inundation, and sediment dynamics, cannot be fully excluded. Future studies incorporating long-term monitoring and concurrent control sites would allow a more rigorous separation of removal effects from year-to-year environmental variability.