Quanzhou Bay is an open bay in Fujian Province, southeastern China. It is fringed by a large area of wetlands and was designated as the Quanzhou Bay Estuary Wetland Natural Reserve in 2003. Since 1984, the invasion of S.alterniflora has been observed here, which caused many problems in the following decades, such as the degradation of native mangrove community, the decrease in benthic biomass and the loss of bird habitat (Wang et al., 2008; Cui et al., 2023; Yang et al., 2014). To prevent the further degradation of native community, the restoration by replacing the invasive S.alterniflora with mangrove plants has been carried out since 2022. S.atlerniflora in the study area was removed by mowing and plowing in later September 2022 and replaced by mangrove plants in early November 2022 (Figure 1), the details of mangrove plants are shown in Supplementary Table S1.

Field samplings were carried out in September 2022 (Sep 2022), October 2022 (Oct 2022), December 2022 (Dec 2022), March 2023 (Mar 2023), June 2023 (Jun 2023) and November 2023 (Nov 2023). The samplings in Oct 2022 and Dec 2022 were carried out about 10 d after the clear of S.alterniflora and the plants of mangrove seedlings. The sampled stations are shown in Figure 1, among with the control station was located in mudflat and was without any plants throughout the investigation. In order to avoid the residual signal of the S.alterniflora burned deeply in situ, surface soil (~5 cm) were collected with a sampling tube (2.9 cm inner diameter) to detect the content and the δ¹³C and δ¹⁵N of SOM (δ¹³C_SOM, δ¹⁵N_SOM). Samples of the neighbored seawater and plants (S.alterniflora and Kandelia candai) were also collected for the detecting of δ¹³C and δ¹⁵N. The above samples were all kept at −20 °C before analysis. Another three replicates of 30 g surface soil were collected and homogenized under anoxic conditions. Subsamples of 5 g soil were transferred to 50 mL tubes, cooled with dry ice and transported under an Ar atmosphere to the laboratory for molecular analyses.

Triplicate soil samples were collected in a 1 m × 1 m square at each station and mixed evenly before analysis. Suspend organic matters in water samples were collected by filtering 4 L of water through GF-75 filters (47mm, 0.3 μm pore size). The filters were previously combusted for 4 h at 450°C to remove organic matter. All water samples were then filtered through a 200 μm mesh sieve to remove large detritus. In the laboratory, the collected soil, plants and membrane samples were treated with HCl vapor (48 h) to remove inorganic carbon and then dried at 60°C. The dried samples were then fully pulverized, the contents of total organic carbon (TOC). total nitrogen (TN), and the δ¹³C and δ¹⁵N values of the samples were then measured using a Finnigan Delta V Advantage isotope ratio mass spectrometer interfaced with a Carlo Erba NC 2500 elemental analyzer, with an analytical precision <0.2‰.

The biomass and the composition of benthic microbes, and the abundance functional genes for microbial CO₂ fixation (cbbM and cbbL) were observed in Nov 2023. Generally, soil genomic DNA was extracted using E.Z.N.A. Soil DNA Kit (Omega Bio-tek, Inc., USA), the concentration and quality of the genomic DNA were checked by NanoDrop 2000 spectrophotometer (Thermo Scientific Inc., USA). The V3-4 hypervariable region of bacterial 16S rRNA gene were amplified with the universal primer 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACNNGGGTATCTAAT-3’). For each sample, 8-digit barcode sequence was added to the 5’ end of the forward and reverse primers (Allwegene Company, Beijing). The PCR was carried out on a ABI 9700 PCR instrument (Applied Biosystems, USA) using 25 μL reaction volumes, containing 12.5 μL 2×Taq PCR MasterMix (Vazyme Biotech Co.,Ltd, China), 3 °C μL (2ng μL^-1), 1 μL Forward Primer(5 μM)、1 μL Reverse Primer (5 μM), 2 μL template DNA, and 5.5 μL ddH₂O. Cycling parameters were 95 °C for 5 min, followed by 28 cycles of 95 °C for 45 s, 55 °C for 50 s and 72 °C for 45 s with a final extension at 72 °C for 10 min. The PCR products were purified using a Agencourt AMPure XP Kit (Beckman Coulter, Inc., USA). Sequencing libraries were generated using NEB Next Ultra II DNA Library Prep Kit (New England Biolabs, Inc., USA) following the manufacturer’s recommendations. The library quality was assessed by Nanodrop 2000 (ThermoFisher Scientific, Inc., USA), Agilent 2100 Bioanalyzer (Agilent Technologies, Inc., USA), and ABI StepOnePlus Real Time PCR System (Applied Biosystems, Inc., USA), successively. Deep sequencing was performed on Illumina Miseq/Nextseq 2000/Novaseq 6000 (Illumina, Inc., USA) platform at Beijing Allwegene Technology Co., Ltd. After the run, image analysis, base calling and error estimation were performed using Illumina Analysis Pipeline Version 2.6 (Illumina, Inc., USA).

Qualified sequences were clustered into operational taxonomic units (OTUs) at a similarity threshold of 97% use Uparse algorithm of Vsearch (v2.7.1) software. The BLAST tool was used to classify all OTU representative sequences into different taxonomic groups against Silva138 and GenBank non redundant nucleus database(nt) Database, and e-value threshold was set to1e-5. QIIME (v1.8.0) was used to generate rarefaction curves and to calculate the richness and diversity indices based on the OTU information. Based on the results of taxonomic annotation and relative abundance, R (v3.6.0) was used for bar-plot diagram analysis.

The expression of functional genes for microbial CO₂ fixation (cbbM, cbbL) and nitrogen cycle (nxrA, nxrB, AOA-amoA, AOB-amoA, nirS, nirK, nosZ, nrfA, narG, norB and nifH) was analyzed by quantitative real-time PCR (qRT-PCR). Generally, the total RNA samples were isolated by TRIzol reagent (TIANGEN BIOTECH, Beijing). Then the RNA purity and concentration was measured by using the NanoPhotometer spectrophotometer (IMPLEN, CA, USA). After detecting, cDNA was synthesized using 2 μg RNA using the PrimeScriptTM RT reagent Kit with gDNA Eraser (TaKaRa). The specific primer for function genes were designed using Primer 5.0 by Allwegene Technology (Allwegene Technology Co., Ltd. Beijing, China), the details are shown in Supplementary Table S2. The 16S rRNA was used as internal reference gene. qRT-PCR reaction was performed using SYBR® Premix Ex TaqTM II (Tli RNaseH Plus) and was conducted on ABI 7500 Real-time Detection System (Thermo Fisher Scientific, USA). The PCR reaction was carried out with the following reaction conditions: 95°C for 30 s; followed by 45 cycles of 95°C for 5 s, 60°C for 40 s. Samples for qRT-PCR were run in 3 biological replicates with 3 technical replicates and the data were represented as the mean ± SD (n = 3) for Student’s t-test analysis. The relative gene expression was calculated using the 2^-△△CT algorithm (Pfaffl, 2001).

Pearson’s correlations analysis and ANOVA test were conducted using the Statistical Package for Social Sciences program (version 19.0) and R (v3.6.0). The source analysis of organic matters based on δ¹³C, δ¹⁵N values were conducted by IsoSource (version 1.3), with the increment of 1% and the tolerance of 0.1% (Phillips and Gregg, 2003). The cluster analysis of OTUs and the diversity analysis of benthic community were completed by Mothur(v1.48.0) and R(v3.6.0).

The geo-chemical soil parameters are shown in Supplementary Table S3. The contents of the total organic carbon (TOC) and nitrogen (TN) in soil ranged from 0.17 ~ 1.89% (0.64 ± 0.27%, n = 36) and 0.02 ~ 0.15% (0.08 ± 0.03%, n = 36) during the investigation. Except for the maximum observed at A (1.89%) in Dec 2022, the contents of TOC at each stations decreased from Sep 2022 to Dec 2022 and recovered from Dec 2022 to Nov 2023, with less spatial variation (Figure 2A). Higher TN were also observed in Sep 2022, which decreased rapidly following the removal of S. alterniflora and slightly varied in the following months (Figure 2B). The C: N_SOM ranged from 5.3 ~ 38.8 (10.4 ± 5.3, n = 36), which increased rapidly in Oct 2022, and kept high in the following investigation. The carbon and nitrogen storage in the surface 5 cm soil were observed by multiply the contents of TOC and TN with the bulk density (Supplementary Table S3). The surface carbon storage in the study area ranged from 106.1 ~ 656.9 t ha^-1 (292.4 ± 120.1 t ha^-1, n = 36), with similar temporal variation with TOC (Figure 2C). The surface nitrogen storage ranged from 2.3 ~ 13.6 t ha^-1 (7.0 ± 2.8 t ha^-1, n = 36), which decreased rapidly after the removal of S.alterniflora, and slightly increased in the following mouths (Figure 2D). Generally, the spatial variation of SOM content was not obvious throughout the investigation.

The δ¹³C_SOM and δ¹⁵N_SOM in the study area ranged from −28.60‰ ~ −20.96‰ (−23.15 ± 1.64‰, n = 36) and 4.03‰ ~ 14.42‰ (6.41 ± 2.57‰, n = 36) respectively. δ¹³C_SOM varied slightly during the investigation (Figure 3), except for the minimum observed at A in Dec 2022 (−28.60‰). Whereas the δ¹⁵N_SOM in Jun 2023 was obviously higher (Figure 3) than that in other seasons.

The absolute abundance of 16S rRNA, cbbM and cbbL ranged from 1.3×10⁷ ~ 3.7×10⁷ Copy g^-1, 2.6×10⁸ ~ 6.2×10⁸ Copy g^-1 and 8.9×10⁶ ~ 1.1×10⁷ Copy g^-1 respectively. The abundance of 16S rRNA and cbbL has less spatial variation, while the abundance of cbbM at C, D and the control station are obviously higher than that at the other stations (Figure 4).

The composition of benthic microbes at phylum, class, order, family, genus and species level are shown in Figure 5. Most of the benthic microbes in the study area belonged to proteobacteria (Figure 5A), among which Gammaproteobacteria was the major one, accounting 22.7% ~ 29.7% (26.3 ± 2.8%, n=6) of the total biomass (Figure 5B). At the species level, most of the observed microbes were uncultured species, accounting more than 60% of the total biomass at all stations (Figure 5F). Cluster analysis of the observed OTUs indicates that the microbial composition at A and B are similar with each other, while those at D and E are more similar to that at the control station (Figure 5). However, the microbial diversity at the sampled stations has no obvious difference (Supplementary Table S4).

The relative expression of cbbL, cbbM, and functional genes for nitrogen cycle (nxrA, nxrB, AOA-amoA, AOB-amoA, nirS, nirK, nosZ, nrfA, narG, norB and nifH) was calculated by 2^-△△CT algorithm. The ratio of the 2^-△△CT observed at A ~ E to that at the control stations were used to indicate the spatial variation of the expression of the functional genes in the study area. The results are shown in Figure 6. Generally, the expression of cbbL at A ~ E were lower than that at the control station, while the expression of cbbM at A, B and E were higher (Figure 6A). The expression of functional genes for ammonia oxidation at A ~ E were lower than that at the control station (Figure 6B), while that for denitrification were obviously higher (Figure 6C). The maximum expression of nitrification and nitrogen fixation were observed at A, while that at the other stations were lower or similar to the control station (Figure 6D).

During our investigation, the carbon and nitrogen storage in the surface soil (5 cm) exhibited a rapid decrease following the removal of S.alterniflora, but gradual increased since the plant of mangroves (Figures 2C, D). At most of the sampled stations, the carbon storage observed in Nov 2023 were higher than that in Sep 2022. This indicate that mangrove plants may have greater contribution to SOM than S.alterniflora. Previous studies have reported that the production of C3 plants tend to exhibit greater isotopic fractionation between ¹²C and ¹³C (Reibach and Benedict, 1977; Roeske and O’Leary, 1984; Drake, 2014). Hence the photosynthetic products derived from C3 plants are usually characterized by lower δ¹³C than that derived from C4 plants. In this study, δ¹³C_SOM was found to be negatively correlated with TOC (Figure 7). This also suggests that mangrove plants (C3 plants) may have replaced S.alterniflora (C4 plants) to be the dominant source of SOM in the study area. To test this hypothesis, δ¹³C_SOM and δ¹⁵N_SOM were used as indicators for tracing the source of SOM (Phillips and Gregg, 2003).

Apart from S. alterniflora and mangrove plants, the production of benthic algae and the imports of external organic matter (by seawater or freshwater) may also have great contribution to SOM (Zhang et al., 2010; Alongi, 2014; Gao, 2019). In addition, the study area was located in a region with a dense human population, the imports of anthropogenic sewage may also be considered (Bianchi, 2011; Li et al., 2023). The contribution of the above sources to total SOM was observed in three steps. First, the sources at A ~ E were divided into S.alterniflora, mangrove plants, and “other source”. The control station without S.alterniflora and any mangrove plants was set to observe the δ¹³C and δ¹⁵N of “other source” (Figure 1). Based on multiterminal element mixing method (two-terminal elements in Sep 2022 and Oct 2022), the contributions of the above sources to total SOM at A ~ E were observed. Then, the composition of the “other source”, which was mainly composed by benthic algae, seawater, freshwater and sewage, was analyzed by a four-terminal element mixing method. Finally, by synthesizing the results of the above steps, the contributions of different sources to total SOM were observed. The δ¹³C and δ¹⁵N of all end-elements are shown in Table 1, the contribution of mangroves plants was observed by summing the contribution of all species at each station (Supplementary Table S1).

The change in the major source of the SOM in the study area was as follows (Figure 8): (1) in Sep 2022, most of the SOM was derived from external sources, accounting for >50% of total SOM. Benthic algae was another important source, with a higher contribution (26.1 ~ 32.2%) than the dominant S. alterniflora (5.0 ~ 22.9%). (2) In Oct 2022, after the clearing of S. alterniflora, a simultaneous decrease in the contributions of S.alterniflora and seawater was observed, from 3.4 ~ 21.3% to 0.2 ~ 9.9% and from 21.7 ~ 26.3% to 16.6 ~ 18.4%, respectively. A considerable increasing contribution of sewage was also observed, from 3.1 ~ 4.6% to 30.9 ~ 34.2%. (3) From Dec 2022 to Jun 2023, a constantly increasing contribution of mangroves was observed. In Dec 2022, approximately 81.2% of total SOM at A was derived from mangrove plants. In addition, the contents of TOC (r=0.765, p<0.05, n=15) and TN (r=519, p<0.05, n=15) derived from mangrove plants were positively correlated with the individual numbers of mangrove plants (Supplementary Table S1). (4) In Nov 2023, mangrove plants was still the major source of SOM at all stations. Another point to note is that the contribution of S.alterniflora increased, which was similar with the contribution of mangroves at E. The signal of S.alterniflora was observed throughout the investigation. We propose that the relatively weaker signal observed in Oct 2022 may be derived from the residues of S.alterniflora buried in situ, while the stronger signal starting in Dec 2022 may indicate the possible recurrence of S. alterniflora, which was indeed observed in Jun 2023 and Nov 2023.

The increasing content of SOM indicates that the restoration by replacing the invasive S.alterniflora with mangrove plants may be beneficial for the carbon storage function in the study area. Though the signal of S.alterniflora can be observed throughout the investigation, its contribution to total SOM was not as high as except, even in Sep 2022, while the strong signal of mangrove plants was observed at all stations since Dec 2022 (Figure 8). This indicates that the direct contribution of mangrove plants to total SOM may be stronger than S.alterniflora. The most possible reason for this is the different carbon allocation of mangrove plants and S.alterniflora. It has been reported that more than half of the carbon fixed by mangrove plants may be conserved in soil (Alongi, 2014; Kristensen et al., 2008). However, most of the carbon fixed by S.alterniflora may be converted to its vegetation biomass, rather than being conserved in soil (Hemminga et al., 1996; Castañeda-Moya et al., 2013; Zheng et al., 2023). Even so, the potential contribution of S.alterniflora can not be ignored. In this study, the contribution of seawater decreased following the removal of S.alterniflora (Figure 8), suggesting the strong accumulation of external organic matters by the roots of S.alterniflora (Castañeda-Moya et al., 2013; Zheng et al., 2023). Hence, we propose that the decrease in the carbon and nitrogen storage after the removal of S. alterniflora from Sep 2022 to Oct 2022 may be mainly caused by the loss of external organic matters derived from seawater (Figure 8). However, such loss can be offset with the growth of mangrove plants, as the accumulation of external organic matters by the complex roots of mature mangrove plants may be stronger than that of S.alterniflora (Alongi, 2014; Kristensen et al., 2008). Another reason for the increasing carbon storage in the study area may be the different C:N ratio of organic matters derived from different source. It has been reported that organic matter derived from S.alterniflora and seawater has lower C:N ratio and tannin content than that derived from mangrove plants (Feng et al., 2015; Wu, 2018; Gao, 2019). This can be further confirmed by the increasing contribution of mangrove to SOM, and the correlations between the C: N_SOM, TOC and δ¹³C_SOM (Figure 7). Due to the biological preference for N during the utilization of organic matter, SOM with higher C:N ratio is usually characterized by lower bioavailability and higher stability (Dauwe and Middelburg, 1998; Jílková et al., 2020). Thus, the replacement of S.alterniflora by mangroves may simultaneously increase the content and the stability of SOM, and hence promote the carbon storage function of the study area in long term (Figure 9).

Generally, our results reveal the benefits of the restoration for the carbon storage function of coastal areas. However, the persistent signal of S.alterniflora reveals its long-term effect, highlighting the need for the long-term monitoring of the study area. Our results also indicate that δ¹³C and δ¹⁵N may be effective in the early warning of S.alterniflora invasion and sewage discharge, and may be helpful in the future restoration and monitoring of coastal wetlands.

During our investigation, the negative correlation between the abundance of 16S rRNA and C: N_SOM was observed in Nov 2023 (r=−0.950, p<0.01, n=6), suggesting that C: N_SOM may be efficient in indicating the bioavailability of SOM. The difference in C: N_SOM between A ~ E and the control station increased with the contribution of mangroves (r=0.472, p<0.01, n=30), indicating the significant contribution of organic matters with higher C: N ratio and tannin contents by mangrove plants. Therefore, we propose that higher abundance of 16S rRNA may be observed at stations with lower contribution of mangroves and higher contribution of seawater, algae or S.alterniflora. However, there was no statistically significant correlation between 16S rRNA abundance and the contribution ratio of the present source of SOM in 4.1. This suggests the existence of other potential sources of bioavailable organic matter for the growth and the activity of the benthic microbes in the study area, which mainly composed by heterotrophic proteobacteria (Figure 5). Microbial assimilation of CO₂ may be the major one.

Microbial assimilation of CO₂ is a ubiquitous process in soils (Berg, 2011; Wood et al., 1941), which has been widely reported (Beulig et al., 2015; Yuan et al., 2012; Nowak et al., 2015; Ge et al., 2016; Spohn et al., 2019). RubisCo is the key carboxylating enzyme for the CO₂ fixation based on Calvin-Benson-Bassham cycles, and has been reported to be highly abundant in agricultural, forest, and wetlands (Nanba et al., 2004; Tolli and King, 2005; Nowak et al., 2015). Hence, the abundance of cbbL and cbbM, the functional genes for the form of RubisCO, was widely used to indicate the activity of benthic microbial CO₂ fixation (Nanba et al., 2004; Tolli and King, 2005; Nowak et al., 2015). During our investigation, the activity of benthic nitrogen cycle were also observed by detecting the relative expression of the associated functional genes (Figure 6). The obvious sewage contribution to SOM (Figure 8) and the higher δ¹⁵N_SOM observed in the study area suggests the potential active denitrification and ammonia oxidation in the study area (Craine et al., 2015; Denk et al., 2017), which have been widely reported in coastal areas under the effects of sewage (Reis et al., 2019; Wang et al., 2022). As most of the microbes associated with benthic nitrogen cycle are heterotrophic (Reis et al., 2019), its activity may largely reflect the bioavailability of organic matters in the study area. During our investigation, the abundance of cbbM was about 30 to 60 times higher than that of cbbL (Figure 4). When comparing the relative expression, the maximum expression of cbbL was observed at the control station, while that of cbbM was observed at A (Figure 6). Though there was no statistical correlation between 16S rRNA, cbbM and cbbL abundance, the significant correlations between the expression of cbbM and the activity of the functional genes for nitrogen cycle were observed (Figure 10). This suggests that microbial assimilation of CO₂ may facilitate the benthic nitrogen cycle processes, thus may serve as an important source of bioavailable carbon for the benthic microbes in the study area. Therefore, we propose that though the removal of S.alterniflora may reduce the bioavailability of SOM, the CO₂ fixed by benthic microbes may serve as an alternative source of bioavailable carbon and hence support the benthic community in the study area (Figure 9). However, the effects of microbial utilization of CO₂ on the isotopic composition of SOM remains unclear at present, which need to be considered in future investigation to understand the contribution of microbial CO₂ fixation to the carbon storage in the study area.