Study site and species
The study area is in a brackish marsh in the northern inlet of the Yangtze estuary, China (31°41′N, 121°41′E; Fig. 1b). This site is characterized by a subtropical monsoonal climate. The mean annual temperature (MAT) is 15.3°C, and the mean annual precipitation (MAP) is 1,022 mm. Tides are regularly semidiurnal, with an average tidal height ranging from 2 to 6 m. The brackish marshes at the study site exhibit clear shifts in the dominant plant species along an elevation gradient, with the low marsh dominated by Sc. mariqueter and the middle and high marshes dominated by P. australis35. Spartina alterniflora is native to North America and has now successfully invaded many of the world’s coastal ecosystems76. In the Yangtze estuary, S. alterniflora has altered the species composition of native marshes by invading the communities of Sc. mariqueter and P. australis, forming new plant communities composed exclusively of S. alterniflora in the low and middle marshes and exclusively P. australis in the high marsh, with striking zonation or a total conversion to S. alterniflora alone (Fig. 1b).
In the Yangtze estuary, S. alterniflora and P. australis have different phenologies. Although S. alterniflora and P. australis sprout from clonal ramets at the same time in March, P. australis green up (from April to May) and senescence approximately one month earlier than S. alterniflora (from December to January), respectively77,78. S. alterniflora and P. australis also have different morphologies. The ramets of S. alterniflora are denser (316 ± 17 ramets m− 2) but shorter (115 ± 8 cm) than those of P. australis (50 ± 6 ramets m− 2 and 298 ± 17 cm)23. These morphological differences between S. alterniflora and P. australis are particularly evident during their peak biomass season in October. Thus, P. australis patches in the S. alterniflora-invaded marshes can then be easily identified from the satellite images on Google Earth and high-resolution images taken with a drone (Supplementary Fig. 1 and see following methods).
Increasing density and size of P. australis patches over time.
The density and size of P. australis patches in an S. alterniflora-invaded marsh in the Yangtze estuary were estimated from 2014 to 2020 (i.e., before and after the reduction in N inputs that started in 2015). Before outlining P. australis patches, we determined the distribution area of S. alterniflora in the Yangtze estuary by using a pixel- and phenology-based mapping algorithm, time series Landsat 5/7/8 images, and the Google Earth Engine cloud computing platform. The pixel- and phenology-based algorithm was developed based on the differences in phenological characteristics among S. alterniflora and other native saltmarsh plants in two temporal windows (April-May and December-January) through the threshold values and frequencies of two vegetation indices (NDVI and enhanced vegetation index (EVI)) and one water index (land surface water index, LSWI) (for detailed information, please see78).
When evaluating high-resolution images (ca. 1 m) taken with a DJI Phantom 3 drone during low tide on clear, cloud-free days in December 2021, we found that P. australis appeared as blurry fungus-like spots when it flowered (Supplementary Fig. 1). These spots were outlined in an S. alterniflora-invaded marsh (up to ~ 14 km2 in 2020) from Google Earth during the flowering period of P. australis (from October to December). In 2020, we selected sixty areas (~ 330 × 420 m) in the invaded marsh, and then five hundred P. australis patches out of a total of 965 in these areas were validated with the high-resolution images taken with the drone in 2021 by determining the repeated presence of these patches at given locations between images of different resolutions (part of the validation is shown in Supplementary Fig. 1). The accuracy of the outlined P. australis patches reached 90.6% (± 1.4), with confounding caused mainly by omission errors. This high accuracy suggested that this approach was suitable for analysing the temporal dynamics of P. australis patches in the S. alterniflora-invaded marsh at an annual resolution.
We then outlined the P. australis patches within the distribution area of S. alterniflora on the historical satellite images in reverse chronological order (i.e., from 2020 to 2014), which made it less likely to miss P. australis patches when they were smaller in earlier years. The results for 2015 and 2017 were not obtained because of the lack of available satellite images during the P. australis flowering period. The satellite images with outlined P. australis patches were used to estimate the density and size of P. australis in the S. alterniflora-invaded marsh. Then, twenty 90,000 m2 areas (300 × 300 m) with more than 40% plant cover were randomly selected in the marsh at least 500 m apart to estimate the density and size of P. australis patches. Considering that the competitive ability of S. alterniflora and P. australis decreases and increases with marsh elevation56, the density and size of P. australis patches may increase significantly with marsh elevation. We further determined the spatial patterns in the density and size of P. australis patches across marsh elevation by randomly picking 20 grids (300 × 300 m) in this study area in only 2020 due to the lack of available marsh elevation data. Marsh elevation data in 2020 were obtained from the Enterprise Edition 91 Weitu Assistant.
We evaluated the temporal changes in the density and size of the P. australis patches in the invaded marsh using LME models, with year included as a fixed effect. Considering that both S. alterniflora and P. australis are perennial plants, we also included year as a random effect to account for the interannual impact. We further analysed the relationships between the density and size of P. australis patches and marsh elevation using simple linear regression. All analyses were conducted in R (v4.1.2) using the “nlme” package.
Decreasing growth performance of S. alterniflora over time.
To monitor the temporal effects of decreasing N availability on plant growth, we calculated the maximum NDVI of both plants over time at the landscape scale and measured plant biomass, stem density, plant height, and leaf N availability over time under ambient N conditions in the SINE experiment. We calculated the maximum NDVI from 2014 to 2020 following the methods of78,79. Briefly, Landsat ETM+/OLI images were taken at low tide on clear, cloud-free days in the fast-growing phase (June–August) of S. alterniflora and P. australis and then were used to calculate the NDVI of 15 randomly picked sites in the monospecific stands of S. alterniflora and the large P. australis patches (diameter, > 40 m) that were present throughout the period from 2014 to 2020.
In the SINE experiment, twenty 4 × 4 m plots were established in 2017. Each set of 10 plots, spaced 20 − 30 m apart, was established at similar elevations in monospecific stands of S. alterniflora and P. australis and assigned alternately to two treatments (n = 5 replicates): ambient N or N enrichment. To simulate the eutrophic conditions in the Yangtze estuary, the N enrichment plots were enriched using a dissolved inorganic N solution (NH4+-N:NO3−-N = 4:1 by weight) at a rate of 150 g N m− 2 yr− 1 monthly from March through December starting in 2017. Nitrogen was applied during neap tides to reduce nutrient removal by the tides and to ensure that plants had enough time to absorb the nutrients. For more details, please see23.
We deployed three small quadrats (20 × 20 cm for S. alterniflora and 50 × 50 cm for P. australis, different quadrat sizes due to large differences in their density and height (Fig. 2a-c)) in each plot to quantify plant growth performance from 2017 to 2020, as described in previous studies80,81. In the peak biomass season (October), fresh and healthy leaves of both S. alterniflora and P. australis were collected, oven-dried and ground to measure carbon (C) and N concentrations and the N:C ratio using a nutrient analyser (vario MACRO cube; Elementar, Hanau, Germany). We harvested aboveground biomass, counted stem density, and measured plant height for both S. alterniflora and P. australis at ground level in each quadrat. We then collected one soil block (10 × 10 cm) from each quadrat by carefully excavating to a depth of 30 cm to determine root biomass. The soil blocks were gently washed to remove soils and adhering materials, and belowground tissues were then sorted into roots and rhizomes or other plant material.
To determine the edaphic factors, we collected three soil cores (diameter, 4 cm; depth, 30 cm) from each quadrat. The soil samples were then pooled and mixed thoroughly, cleaned to remove stones and plant materials by hand, and weighed before and after being air-dried to determine the water content. The samples were then ground and sieved (mesh size < 0.2 mm) to measure the dissolved inorganic nitrogen (DIN) concentration and salinity. Soil DIN was extracted with 0.5 mol/L K2SO4 (1:5, w:v) and analysed colorimetrically for NH4+-N and NOx−-N on a microplate reader (Synergy 2; BioTek, Winooski, Vermont, USA). Soil salinity (1:5, w:v) was measured with a conductivity metre (SevenExcellence S479-uMix; Mettler-Toledo, Greifensee, Switzerland). We also compiled MAT from 2017 to 2020 from the Climatic Data Center of the National Meteorological Information Center in China (https://data.cma.cn/).
We evaluated the temporal changes in the maximum NDVI of each species using LME models, with year included as a fixed effect. In this model, year was also treated as a random effect to account for the interannual impact. We compared temporal patterns of growth performance (aboveground biomass, belowground biomass, root:shoot ratio, stem density, plant height, and leaf N availability) between S. alterniflora and P. australis using LME models, in which species, year, and their interactions were included as fixed effects, and year and plot ID were included as random effects. Multivariate regression analyses were performed separately for S. alterniflora and P. australis with LME models to assess relationships between growth performance and environmental variables (MAT, soil salinity, DIN, and water content).
The competitive advantage of S. alterniflora decreases with decreasing N availability.
Beginning in March 2021, we conducted a field experiment to compare the competitive advantage of S. alterniflora on co-planted P. australis under ambient N and N-enriched conditions. Eight plots (2 × 1 m, spaced > 2 m apart) were established in S. alterniflora-dominated marshes outside of the SINE experiment at similar elevations. We transplanted 8 S. alterniflora turfs with similar plant heights (~ 10 cm) from N-fertilized and unfertilized plots in the SINE experiment, which has run for 4 years since 2017, using PVC pipes (diameter, 10 cm; depth, 30 cm) and then deployed the turf randomly into each 2 m2 plot, with the top edges of the PVC pipes flush with the marsh substrate. The pipes were predrilled with 42 holes (diameter, 1 cm) into their walls to allow the soil pore water to flow in and out. The stem densities of these turfs were also counted separately before the experiment.
Four of the 8 plots were randomly assigned to each of two treatments: ambient N and N enrichment. N enrichment plots were fertilized with N monthly from March through December 2021 at a rate of 150 g N m− 2 yr− 1 following the method used in the SINE experiment to simulate the severe level of eutrophication in the Yangtze estuary that occurred before 2017. The ambient N plots without fertilization were used to represent the decreased N availability in the Yangtze estuary after 2017.
Considering that the competitive ability of S. alterniflora may decrease with age39, we determined the confounding effect of population age on the effects of N levels on the competitive ability of S. alterniflora in this transplant experiment. The S. alterniflora populations from both the N-fertilized and unfertilized plots in the SINE experiment were approximately 10 years old in 2020, as determined by historical images on Google Earth. Using this identification method, we also collected S. alterniflora turf of different ages (2 and 5 years) from S. alterniflora populations near the SINE experiment (100 ~ 300 m apart). The interval of population ages (2, 5, and 10 years) was similar to the period of fertilization (4 years) experienced by S. alterniflora in the N-fertilized plots. Finally, 32 S. alterniflora turfs (4 populations of different ages × 2 N levels × 4 replicates) in total were transplanted into these 2 m2 plots.
As N fertilization had no significant effect on the aboveground biomass of P. australis in the SINE experiment23, three disconnected clonal ramets of P. australis of similar height (~ 10 cm) were transplanted carefully from the unfertilized plots in the SINE experiment into each of the 32 turfs without harming the roots of S. alterniflora. Any ramets of P. australis that died within 1 month after being transplanted were replaced with healthy ramets. Considering that the aboveground parts of S. alterniflora occurring outside the PVC pipes in the marsh may affect the growth of plants inside the PVC pipes through light competition, we cleared them all to ground level within 1 m away during the growing season (from March to October). In October, the plants in each PVC pipe were harvested at ground level and sorted by species (S. alterniflora or P. australis). Their aboveground biomass was weighed separately after oven drying at 70°C to a constant mass. We analysed the effects of N enrichment, species, and population age on aboveground biomass and the effect of N enrichment and population age on P. australis dominance (represented by the proportion of aboveground biomass in the mixed stand) using LME models with the initial density of S. alterniflora in the transplanted turf as a random effect.
Reduced N availability shifts competitive outcomes from exclusion to coexistence
As the competitive advantage of S. alterniflora was primarily affected by N levels (Fig. 3), the competitive advantage of S. alterniflora before and after the reduction in N inputs could then be simulated in N-enriched (150 g N m− 2 yr− 1, as in the SINE experiment) and ambient N conditions, respectively. To assess the mechanisms by which N levels affect the competitive advantage of S. alterniflora over P. australis, we conducted another N-enrichment field experiment using three plant groups: monospecific stands of P. australis (blue quadrat); monospecific stands of S. alterniflora (orange quadrat); and mixed stands with both P. australis and S alterniflora (black quadrat) (see Fig. 4a). The three plant population quadrats were grown under N-enriched conditions (solid quadrat) or lower N ambient conditions (open quadrat), and the various quadrats were arranged along radial transects in a large, approximately circular patch of P. australis (diameter, ~ 40 m) in an S. alterniflora-dominated marsh in the Yangtze estuary in 2021 (n = 5 replicates of 3 group arrangements and 2 N conditions, for a total of 30 quadrats). The transects were established inside, outside, and at the margin of the P. australis patch (i.e., monospecific stands of P. australis or S. alterniflora and a mixed stand of both species). Ten sets of 3 quadrats (1 × 1 m) were established along each transect, with adjacent quadrats spaced 10 m apart. In each quadrat of the mixed stands, S. alterniflora and P. australis had similar stem densities (Fig. 4d). There were 5 replicates of each treatment, and all plots along each transect were assigned alternately to one of two treatments: ambient N or N enrichment.
In 2021, N was applied to the plots in the enrichment treatment following the method of the SINE experiment23,65. Plant survival was recorded separately for each species in a randomly positioned quadrat (50 × 50 cm) established in each plot during the growing season (from March through October) in the same year. Aboveground biomass and stem density were monitored in both the mono- and mixed stands during the peak biomass season (October) before (2020) and after N enrichment (2021). To quantify aboveground biomass, the plants in each quadrat were harvested at ground level and sorted into S. alterniflora and P. australis. Stem density was quantified, and aboveground biomass was determined after oven drying at 70°C to a constant mass. We then evaluated niche and fitness differences by using the aboveground biomass of S. alterniflora and P. australis in both the monospecific and mixed species stands in October 2020 and 2021 (for detailed calculation, please see the following subsection). A series of ANOVAs were performed with LME models to estimate the effects of N enrichment on niche and fitness differences between the two species separately, in which case plot ID was included as a random effect; to estimate the effects of N enrichment, species, and sampling time on plant survival in the mixed stands during the growing season, in which case sampling time and plot ID were included as random effects; and to estimate the effects of N enrichment and species on stem density and aboveground biomass during the peak biomass season, in which case plot ID was included as a random effect.
Measurement of the mechanisms of coexistence
The growth rates of S. alterniflora and P. australis before (2020) and after (2021) N enrichment in the monospecific (µi,0) and mixed stands (µi,j) were calculated as follows38:
$${\mu }_{i,x}=\frac{1}{T}·\frac{{AGB}_{2021}}{{AGB}_{2020}}$$
where x = 0 or j, T is 365 d, and AGB2020 and AGB2021 are the aboveground biomass of the focal species (i.e., S. alterniflora and P. australis) in the monospecific and mixed stands in 2020 and 2021, respectively.
We estimated niche differences (NDs) and relative fitness differences (RFDs) based on species sensitivity to competition as previously described38,82. We calculated the sensitivity of each species to competition for each replicate, which was defined as the amount by which the per-individual growth rate for each species was reduced when grown in the mixed species treatment (µi,j) relative to the growth rate that each species achieved in the respective monospecific stands (µi,0):
$${S}_{i,j}=\frac{{\mu }_{i, 0}-{\mu }_{i, j}}{{\mu }_{i,0}}$$
When Si,j < 1, i can outcompete j, and when Si,j > 1, i cannot invade j. The ND was then calculated as one minus the geometric mean of the species’ sensitivities:
$$ND=1-\sqrt{{s}_{i,j}{s}_{j,i}}$$
An RFD causes species sensitivities to diverge. For two species indexed such that Si,j ≥ Sj,i, the RFD is given as follows:
$$RFD=\sqrt{{S}_{i,j}/{S}_{j.i}}$$