Experimental warming reduces ecosystem resistance and resilience to severe flooding in a wetland.

Climate warming and extreme hydrological events are threatening the sustainability of wetlands across the globe. However, whether climate warming will amplify or diminish the impact of extreme flooding on wetland ecosystems is unknown. Here, we show that climate warming significantly reduced wetland resistance and resilience to a severe flooding event via a 6-year warming experiment. We first found that warming rapidly altered plant community structure by increasing the dominance of low-canopy species. Then, we showed that warming reduced the resistance and resilience of vegetation productivity to a 72-cm flooding event. Last, we detected slower postflooding carbon processes, such as gross ecosystem productivity, soil respiration, and soil methane emission, under the warming treatment. Our results demonstrate how severe flooding can destabilize wetland vegetation structure and ecosystem function under climate warming. These findings indicate an enhanced footprint of extreme hydrological events in wetland ecosystems in a warmer climate.

INTRODUCTION

Wetland ecosystems are changing rapidly in response to climate driving factors, especially climate warming (–) and extreme hydrological events (, ). As a typical extreme hydrological event, the abrupt occurrence of flooding has profound impacts on wetland structure and function (–). Species-level plant biomass or ecosystem-level productivity commonly decreases after flooding events in multiple ecosystems, including croplands (–), grasslands (, ), forests (), and wetlands (). However, the future changes of global wetlands are particularly uncertain because of the unknown role of extreme flooding in shaping ecosystem responses to climate warming. Resistance and resilience are two critical components of ecosystem stability (). Resistance describes the ability of the ecosystem to maintain its original level during climate extremes, and resilience measures the rate at which a process returns to its pre-extreme level (, ). Whether wetland resistance and resilience to severe flooding will be amplified or diminished by future climate warming remains unclear.

Because temperature regulates the rate of almost all biochemical processes, climate warming can increase vegetation productivity (, ) and shift plant species composition in many ecosystems (–). In grasslands and forests, experimental warming usually triggers rapid plant ecophysiological responses and slow changes in species composition (, ). Although there are only a few relevant long-term experiments (), the shifts in plant species composition can mediate and even reverse the short-term responses of ecosystem productivity to climate warming (, ). Climate warming can rapidly alter wetland vegetation productivity and species composition because of the rapid species-specific adaptations (, , ). Thus, the warming effect on resistance and resilience of wetland productivity to extreme flooding could depend on the changes in plant structural traits such as canopy height and stem density ().

Over the last three decades, more than 200 manipulative experiments have explored the ecosystem responses to climate warming (fig. S1) (). Although there are >30 published warming studies on global wetland ecosystems, no study has reported the impact of warming on resistance and resilience of wetland productivity to extreme flooding. Given that vegetation productivity is the primary process of the wetland carbon cycle, the uncertain role of severe flooding on vegetation productivity also limits our projections for carbon cycling in the global terrestrial-aquatic interface under future climate warming (, ).

Here, we report results from a 6-year wetland warming experiment in the Yellow River delta of China. The warming plots were heated continuously by infrared heaters from 1 November 2014, with the soil temperature increased by 2.4°C over the 6 years (Fig. 1A). We measured environmental factors (e.g., soil microclimate and soil salinity), plant species composition, and many ecosystem functions, such as net primary productivity and soil carbon effluxes (i.e., CO2 and CH4). Although annual total precipitation was similar among the 6 years (Fig. 1B), there was an extreme flooding event owing to a 347-mm rain event during August to September in 2016 (Fig. 1C). The severe flooding event allowed us to quantify flooding resistance and resilience of wetland vegetation based on the measurements from 2015 to 2020 (details in Materials and Methods) (, ). This study aimed to explore the role of extreme hydrological events in regulating the wetland ecosystem response to climate warming.

Fig. 1.

Temporal dynamics of soil temperature, climate, and the occurrence of flooding by year.

Temporal dynamics of soil temperature and annual precipitation in the experimental site from 2015 to 2020 (A) and the probability distributions of annual precipitation (B) and single maximum precipitation event (C) over 1961–2020. The vertical lines show the annual precipitation (B) and maximum precipitation event (C) in each year from 2015 to 2020.

RESULTS

Flooding resistance reduced by warming

The warming treatment significantly increased soil temperature by 2.4° ± 0.1°C across the 6 years (P < 0.001; Fig. 1A and tables S1 and S2). Phragmites australis and Suaeda glauca, two dominant species, contributed ~96% to total aboveground net primary productivity (ANPP) (fig. S2). In the ambient control plots, the aboveground biomass of P. australis had significant interannual variation, ranging from 657.0 ± 73.8 to 938.5 ± 60.2 g m−2 year−1 in 2017 (table S1). The analysis by repeated-measures analysis of variance (ANOVA) showed that warming decreased the aboveground biomass of P. australis by 62.5% (P < 0.001) across the 6 years, with an insignificant difference among years (P = 0.08; tables S1 and S2). The aboveground biomass of S. glauca ranged from 26.4 ± 6.0 g m−2 year−1 in 2019 to 488.5 ± 34.1 g m−2 year−1 in 2015 (table S1). The positive warming effect on S. glauca varied significantly among years (P < 0.001; tables S1 and S2), ranging from 34.3% in 2016 to 228.0% in 2018. In 2015, the aboveground biomass of S. glauca was enhanced by 667.1 ± 80.3 g m−2 year−1 under warming (P < 0.01; Fig. 2A). In addition, warming decreased the height of P. australis but increased the height of S. glauca (Fig. 2, A to F). Therefore, warming reduced the growth of P. australis but increased that of S. glauca in 2015 (both P < 0.001; Fig. 2, A and G).

Fig. 2.

Plant height and aboveground biomass of two dominant species in control and warming plots during the 6 years.

(A to F) Plant height (cm) and flooding depth (cm) in control and warming plots during the 6 years. (G to L) Aboveground biomass (g dry matter m–2) of two dominant species, i.e., P. australis (PA) and S. glauca (SG), in control (C) and warming (W) plots during the 6 years. Error bars represent SE across replicates (n = 4). **P < 0.01.

In 2016, an extreme flooding event occurred in the middle of the growing season, with a maximum flooding depth of 72.4 cm and event duration of 57 days (Fig. 2B). On the basis of the meteorological data over 1961–2020 in the Yellow River delta, we identified this extreme flooding as the “131-year flood” (figs. S3 and S4). As shown in Fig. 2B, the flooding depth was higher than the canopy height of S. glauca, leading to the mortality of all individuals of S. glauca in both control and warming plots. The change in ANPP under warming was mainly contributed from the aboveground biomass of P. australis (Fig. 2H). In 2016, warming decreased the aboveground biomass of P. australis and ANPP by 53.5 and 42.6%, respectively (both P < 0.01; table S1 and Figs. 2H and 3A). Compared with the average ANPP during the normal years (mean across all nonclimate event years), ANPP in 2016 (extreme year) declined by 20.3% in control plots, less than that in warming plots (table S1). As a result, the negative response of P. australis and the flooding-inhibited positive response of S. glauca (Fig. 2H) jointly led to a 44.7% reduction in resistance of ANPP to the extreme flooding event (P < 0.01; table S1 and Fig. 3B).

Fig. 3.

Warming effect on ecosystem resistance and resilience.

Warming effects on ANPP (g dry matter m−2 year−1) (A) and its resistance (B) and resilience (C) to flooding. Error bars represent SE across replicates (n = 4). **P < 0.01. (D) The relationship between carbon fluxes and ANPP across all plots and years. The solid line represents a statistically significant correlation (P < 0.05) of ANPP with gross ecosystem productivity (GEP; g C m−2 year−1), ecosystem respiration (Reco; g C m−2 year−1), net ecosystem CO2 exchange (NEE; g C m−2 year−1), soil respiration (Rsoil; g C m−2 year−1), soil heterotrophic respiration (Rh; g C m−2 year−1), and root respiration (Rroot; g C m−2 year−1). The squares represent CO2 uptake. The circles represent CO2 or CH4 release. The dashed line shows the insignificant relationship between aboveground respiration (Rabove; g C m−2 year−1) and ANPP (P > 0.05).

Decreased resilience to flooding under warming

There was no extreme rainfall (Fig. 1C) or flooding (Fig. 2, C to F) event from 2017 to 2020, with the maximum flooding depth ranging from 3.1 cm in 2017 to 29.8 cm in 2019 (Fig. 2, C to F). The warming-induced reductions in biomass of P. australis (from 472 g m−2 year−1 in 2017 to 826 g m−2 year−1 in 2018) were significantly greater than the warming-induced increases in biomass of S. glauca (from 44 g m−2 year−1 in 2019 to 240 g m−2 year−1 in 2020) (table S1). Thus, although warming significantly enhanced the biomass of S. glauca in 3 of the 4 years from 2017 to 2020, the magnitude of the biomass increase in S. glauca could not compensate for the decreased biomass of P. australis each year (Fig. 2, I to L).

After the extreme flooding year, ANPP recovered more slowly in the warming plots than in the control plots (Fig. 3A). After the flooding event in 2016, ANPP in the control plots recovered gradually from 62.9% in 2017 to 65.4% in 2020 of ANPP in 2015 (Fig. 3A). In warming plots, ANPP further declined to 24.3, 18.6, and 14.9% in 2017, 2018, and 2019, respectively, of ANPP in 2015 (Fig. 3A). Compared to the mean ANPP over the normal years (i.e., 2015 and 2017–2020), the control plots showed a rapid recovery, from 86.5% in 2019 to 95.7% in 2018 (table S3). The warming plots showed a slow recovery and reached 59.8% of mean ANPP over the normal years (Fig. 3A). As a result, warming reduced ANPP over 2017–2020 (Fig. 3A), leading to a significantly lower resilience (−40.9%; P < 0.01) of ANPP to the extreme flooding event in warming plots than in control plots (Fig. 3C).

Decelerated carbon cycling after extreme flooding

The regression analysis showed that ANPP played an important role in influencing carbon fluxes in the ecosystem (Fig. 3D). A further analysis based on a structural equation model (SEM) further confirmed the dominant role of ANPP in driving the response of ecosystem carbon cycling to climate warming (fig. S5). Across the 3 years from 2018 to 2020, warming significantly decreased NPP (net primary productivity), ANPP, and BNPP (belowground net primary productivity) by 50.1, 54.5, and 33.8%, respectively (all P < 0.01; table S4 and Fig. 4, A and B). The rates of most key carbon fluxes were decreased by climate warming over 2018–2020 (Fig. 4A). Warming not only reduced the gross ecosystem productivity (GEP) by plants (−16.2%; P < 0.001; Fig. 4C) but also respiratory CO2 release of soil microbes (−25.0%; P < 0.01), bulk soil (−21.3%; P < 0.01), and the entire ecosystem (Reco; −10.0%; P < 0.01). Net ecosystem CO2 uptake (−25.7%; P < 0.001; Fig. 4C) and soil CH4 emissions (−84.8%; P < 0.01; Fig. 4E) were both significantly lower in warming plots than in control plots.

Fig. 4.

Effects of warming on carbon processes after the extreme flooding year.

In (A), the red values represent the relative responses of carbon processes to warming over 2018–2020. (B to E) Means ± SE of ANPP (g dry matter m−2 year−1), BNPP (g dry matter m−2 year−1), NPP (g dry matter m−2 year−1), GEP (g C m−2 year−1), Reco (g C m−2 year−1), NEE (g C m−2 year−1), Rsoil (g C m−2 year−1), Rh (g C m−2 year−1), Rroot (g C m−2 year−1), Rabove (g C m−2 year−1), and soil CH4 emission (CH4; mg C m−2 year−1). Error bars represent SE across replicates (n = 4). **P < 0.01.

Warming decreased net ecosystem CO2 uptake by 32.0, 34.6, and 12.9% in 2018–2020, respectively (all P < 0.05; fig. S6A). GEP was decreased by 27.8 and 16.7% under warming in 2018 and 2019, respectively (all P < 0.05; fig. S6B). Warming reduced Reco by 21.8 and 14.7% in 2018 and 2019, respectively (both P < 0.05; fig. S6C). In 2020, warming had no significant effect on GEP or Reco (P = 0.349; fig. S6, B and C). In addition, we calculated the annual net global warming potential (GWP) of greenhouse gases (i.e., CO2, CH4, and N2O). Experimental warming made the annual net GWP from −4914.1 ± 37.03 to −3657.5 ± 106.01 g CO2 equivalent m−2 year−1 (table S5).

DISCUSSION

This study has demonstrated that climate warming can reduce the resistance and resilience of wetland ecosystems to extreme flooding (Fig. 3) by shifting the plant species composition (Fig. 2). In the warming plots, the dominant species changed from P. australis to S. glauca before the extreme flooding event (Fig. 2G). The increased soil salinity under warming can explain the rapid and contrasting responses between the two species (+49.0%; P < 0.01; fig. S7). For example, the warming-induced change in the aboveground biomass of S. glauca (ABSG) was positively affected by soil temperature but not influenced by soil salinity across all plots and years (fig. S8). As shown in the additional salinity tolerance experiment, the warming-induced increase in the aboveground biomass of S. glauca (ABSG) mainly resulted from the higher salinity tolerance of S. glauca than P. australis (fig. S9). However, the warming-induced shift of dominant species from P. australis to S. glauca in 2015 reduced the flooding resistance in 2016 because of the lower height of S. glauca (24.3 ± 2.3 and 33.3 ± 4.2 cm in control and warming plots, respectively; Fig. 2B) than the flooding depth (72 cm; Fig. 2B). These findings indicate a trade-off between plant salinity tolerance and productivity resistance to flooding under climate warming in coastal wetlands. The critical role of plant functional traits in shaping community structure under extreme climate events has also been detected in other ecosystems, such as the Amazonian rainforest () and temperate grasslands (, ).

The sustainability of the wetland functions, e.g., net carbon sink, strongly depends on carbon cycling through multiple processes within the ecosystem (, ). In this study, the reduced resilience of vegetation productivity to extreme flooding triggered a decelerated carbon cycle under climate warming (Fig. 4). The observed lower rates of C fluxes under warming were contrary to most model predictions in global wetlands (, ). In our field experiment, the negative warming effect on net ecosystem CO2 exchange (NEE) results from the greater reduction in GEP than aboveground and root respiration (Fig. 4C). We found that warming significantly reduced light-saturated photosynthesis but increased dark respiration of P. australis (fig. S10). These results suggest that increasing respiration is an essential physiological mechanism for plants to adapt to salinity stress (, ). We further found that the flooding depth in 2016 in our study is close to the projected mean maximum flood depth of global wetlands over 2091–2100 by global land-surface models (fig. S11). Thus, this study indicates that vegetation structure is critical in regulating the direction of wetland ecosystems in response to extreme flooding events.

We also conducted a gradient flooding experiment and found that a 60-cm-depth flooding event could significantly decrease NEE (−11.7%), Reco (−11.6%), and GEP (−11.7%) (P < 0.05; Fig. 5, A to C, and table S6). However, there was no significant response of ecosystem CO2 fluxes after any flooding event (Fig. 5, B, D, and F). The rapid recovery of ecosystem CO2 fluxes after flooding events supports the high flooding resilience of plant community and ecosystem CO2 sink in the control plots of the warming experiment. We further found that the wetland ecosystem under the control treatment was a sink for CO2 and N2O but a source for CH4 over 2018–2020 (table S5). Their net GWP was negative (table S5). However, experimental warming significantly reduced the net GWP of the three greenhouse gases (table S5). This finding implies a weakened role of coastal wetlands in mitigating climate change if more extreme flooding events occur in a warmer climate.

Fig. 5.

Ecosystem CO2 fluxes in a flooding gradient experiment.

Effect of flooding on NEE (A and B), GEP (C and D), and Reco (E and F) during three periods (i.e., preflooding event, midflooding event, and postflooding event).

Note that this study has not used a fully balanced statistical design (i.e., warmed × flooded) to explore the role of flooding in regulating wetland responses to climate warming. An ideal experimental design would include four treatments, including “unwarmed + unflooded,” “warmed + unflooded,” “unwarmed + flooded,” and “warmed + flooded.” Thus, some cautions should be taken to extrapolate our findings based on the warming experiment and an extreme flooding event. However, the findings in this study have implications for predicting the future dynamics of the wetland carbon cycle under climate change. First, vegetation structure is critical in predicting the responses of wetland ecosystems to extreme flooding events. As shown in fig. S12, we found that incorporating plant structural traits could improve a process-based ecosystem model in simulating GPP (gross primary productivity) dynamics in 2016 in our site. Second, soil salinity regulates the responses of vegetation structure and carbon sequestration to climate change in coastal wetlands (–). Third, the state-of-the-art global land-surface models have not fully coupled the hydrological cycle and vegetation dynamics, so they may overlook the negative impacts of extreme flooding events on ecosystem carbon processes (e.g., fig. S13).

In summary, our study provides experimental evidence for an enhanced footprint of extreme flooding events in a wetland ecosystem under climate warming. Many earth system models project that temperature increases will coincide with a higher frequency of severe flooding globally (fig. S11) (, –). Although those models are evolving to include terrestrial-aquatic interfaces and wetland ecosystems (), hydrological processes are still weakly coupled with vegetation dynamics and biogeochemistry (, ). This study also indicates that manipulative experiments are useful to improve the modeling of the wetland carbon cycle in response to climate extremes. Thus, our findings highlight the urgent need to incorporate process-based interactions between climate extremes and vegetation dynamics in future projections of global wetland changes.

MATERIALS AND METHODS

Site description

The experiment locates in a coastal wetland of the Yellow River delta (37°45′50″N, 118°59′24″E) in Dongying City, Shandong Province, China. The Yellow River delta has a warm temperate and continental monsoon climate with distinctive seasons. This site’s mean annual air temperature is 12.9°C, with the mean seasonal temperature ranging from 26.7°C in summer to −2.8°C in winter (). The average annual precipitation over 1961–2020 is 592 mm, with nearly 70% falling from June to September. The saline soil of this experiment site is classified as sandy clay loam to the Chinese classification (). The vegetation community is dominated by P. australis, S. glauca, Suaeda salsa, Apocynum venetum, Tripolium vulgare, Cynanchum chinense, Triarrhena sacchariflora, Sonchus arvensis, and Imperata cylindrical. More detailed information for each plant species is shown in table S7.

Experimental design

A manipulative warming experiment was established in the Yellow River delta on 1 November 2014. We used a randomized complete block design with two treatments (i.e., control and warming) and four replicates for each treatment. The size of each plot was 3 m by 4 m, and the distance between adjacent plots was about 3 m. All the warming plots were heated continuously by infrared heaters (Kalglo Electronics, Bethlehem, PA, USA) suspended approximately 1.75 m above the ground. We suspended one “dummy” heater at the same height in each control plot to mimic the shading effects of the infrared heater. All the heaters under the warming treatments have a mean radiation output of approximately 1600 W during the entire year. Before the treatment, i.e., in June 2014, we surveyed the background states of the plant community in each plot. We selected a 1 m by 1 m quadrat in each plot to measure the number of individuals, mean height, and canopy cover for each species. The 1 m by 1 m quadrat was divided into 100 grids equally. The abundance and canopy cover of each species were summed from all of the 100 grids. This method has been widely used to monitor plant cover in nonforest ecosystems (, , ). The data of species composition and plant traits in each plot are shown in table S8. The one-way ANOVA analyses showed an insignificant difference in the plant community (i.e., community height, total cover, species height, number, and cover) between control and warming plots before the treatment (table S9). The infrared heaters increased the canopy and soil temperature (fig. S14) but had an insignificant effect on air temperature and air humidity (fig. S15). The annual mean air temperatures in 2015, 2016, and 2018 were 13.0°, 13.1°, and 13.0°C, respectively, values that were slightly lower than the 6-year mean air temperature (i.e., 13.3°C). The warmest climate conditions occurred in 2017 (13.8°C), 2019 (13.5°C), and 2020 (13.6°C) (fig. S16). In control plots, soil temperature ranged from 13.3° ± 0.28°C in 2019 to 14.5° ± 0.08°C in 2017. The warming treatment increased soil temperature by 2.4° ± 0.09°C over the 6 years, ranging from +2.9° ± 0.16°C in 2015 to +1.8° ± 0.08°C in 2020 (table S1).

Species salinity tolerance experiment

We conducted an additional experiment to compare the salinity tolerance between the two dominant plant species, i.e., P. australis and S. glauca. First, we transplanted seedlings of the two species into separate pots (one seedling per pot; caliber: 21.0 cm; bottom diameter: 15.0 cm; height: 18.5 cm). The soil in the pots is low-salinity fertile, with a mean electronic conductivity of 0.87 ± 0.05 dS m−1. Then, we applied the salinity treatments of 0, 20, 40, 60, or 80 PSU (practical salinity units, ‰) by controlling the salinity of standing water (~5 cm). Each treatment was replicated four times, and the dynamics of plant growth were measured for 80 days. A nonlinear regression analysis was applied to fit the relationship between soil salinity and the survival time of the two species (fig. S9).

Flooding gradient experiment

In 2017, we conducted a flooding gradient experiment to examine the response of our ambient ecosystem to different flooding events. We used a randomized block design with five flooding depths (i.e., control, 10, 20, 40, and 60 cm). We replicated each flood treatment by nine times with 45 mesocosms (50 cm in diameter, 100 cm in height, and 40 cm inserted into the soil). We maintained the flooding depth throughout the experiment by adding water to compensate for the evaporation and infiltration. The flood treatment lasted from August 8 to 23. We measured GEP, ecosystem respiration (Reco), and their difference (i.e., NEE) every 2 to 3 days before, during, and after the experimental period.

Meteorological data

We installed a sensor of soil microclimate at a depth of 10 cm in the center of each plot. The sensor automatically measures soil temperature, soil moisture, and soil salinity (i.e., electrical conductivity) at a 2-hour interval (EM50 data logger, Decagon Devices, Pullman, WA, USA). Besides, the daily air temperature was measured with a temperature probe (HMP45C, Vaisala, Helsinki, Finland). Daily precipitation was measured by tipping bucket rain gauges (Texas Electronics, Dallas, TX, USA).

Plant community monitoring

We established one permanent 1 m by 1 m quadrat in each plot since April 2015. The height, cover, and abundance of each species in each quadrat were measured. The measurements were done biweekly in each quadrat from April to November each year.

Above- and belowground net primary productivity

We harvested the current-year aboveground biomass in early October to estimate ANPP. The current-year aboveground biomass was clipped at 3 cm above the ground from a 1 m by 1 m square. The dry mass of each species was determined by oven-drying at 70°C to constant weight.

We also used a nondestructive method to validate the biomass response in this experiment. The nondestructive method included four steps: (i) We recorded each species’ abundance, height, and cover when harvesting the biomass every year; (ii) a regression equation between plant traits and harvested biomass was developed for each species under both control and warming treatments; (iii) the abundance, height, and cover of each species in the permanent 1 m by 1 m quadrat were measured simultaneously during the biomass harvest; and (iv) we used the regression equations from the second step to estimate the biomass in each permanent 1 m by 1 m quadrat. The selected regression equations and the selective criterion are shown in fig. S17. The nondestructive ANPP was calculated as the sum of aboveground biomass in all species. As shown in figs. S18 and S19, the findings are consistent between the destructive and nondestructive methods.

We estimated BNPP using the root ingrowth core method (). The end of the plant-growing season usually occurs in October in the Yellow River delta. We first took out soil cores with a diameter of 10 cm and a depth of 0 to 40 cm from each plot in early November of the previous year. Then, we immediately filled those cores with sifted rootless soil from the same depth outside the plot. In late October, we used a smaller soil auger (8 cm in diameter) to collect root samples at the center of the original root ingrowth holes. The dry mass of the root was determined by oven-drying at 70°C to constant weight.

Measurements of ecosystem-level CO2 fluxes

Ecosystem CO2 fluxes were measured by an infrared gas analyzer (LI-6400, LI-COR, Lincoln, NE, USA) attached to a transparent chamber (0.5 m in diameter and 0.6 m in height). In each plot, one circle frame (0.5 m in diameter and 0.1 m in height) was permanently inserted into the soil at a depth of 7 cm in April 2017. During the measurements, one small electric fan was running continuously to mix the air inside the chamber. We used another high chamber (0.5 m in diameter and 1.2 m in height) for the measurements when the plant height was over 0.6 m. Nine consecutive CO2 concentrations were taken at 10-s intervals after the gas concentration was linearly increasing within the chamber. NEE was calculated as the rate of CO2 concentration change over time. Then, the chamber was ventilated and reseated on the frame and covered with an opaque cloth for the measurement of Reco. GEP was calculated from NEE and Reco. In this study, positive NEE (or GEP) values represent CO2 uptake. Ecosystem CO2 fluxes were usually measured three times per month on clear and sunny days from 9:00 to 11:00 a.m. from November 2017 to November 2020.

Measurements of soil respiration and methane emission

Soil CO2 flux was measured twice each month using an Li-8100 connected to a soil chamber (LI-COR Inc., Lincoln, NE, USA). A cylindrical polyvinyl chloride (PVC) collar (21.3 cm in diameter and 11.4 cm in height) was inserted 2 to 3 cm into the soil in each plot. Living plants inside the PVC collars were clipped at the soil surface at least 1 day before the measurements to eliminate the effect of aboveground plant respiration. The clipped plant material was left inside the collars. We inserted a deep soil collar (diameter: 21.3 cm; height: 35 cm) into the soil to exclude root respiration in each plot. The soil CO2 flux in the deep collar represented soil heterotrophic respiration (Rh). Soil CH4 flux was measured on the PVC collars with an LGR ultraportable greenhouse gas analyzer (UGGA; Los Gatos Research Inc., San Jose, CA, USA). All soil CO2 and CH4 fluxes were measured from 9:00 to 10:00 a.m. from November 2017 to November 2020. The root respiration (Rroot) was calculated as the difference between total soil respiration (Rsoil) and heterotrophic soil respiration (Rh). Aboveground plant respiration (Rabove) was calculated as the difference between Reco and Rsoil.

Calculation of ecosystem resistance and resilience to extreme flooding

Resistance and resilience are two components of ecosystem temporal stability (, ). They have been defined as the ability of an ecosystem to maintain its state and recover from disturbances (, ). In this study, resistance describes the ability of the ecosystem to maintain its original levels during the extreme flooding year. Resilience measures the rate of a process recovering to its pre-extreme flooding level. Similar to the definition in previous studies (, ), we defined resistance (Ω) and resilience (Δ) aswhere ‾ Yn is ANPP during the normal years (i.e., mean ANPP across all nonclimate event years). Ye and Ye+1 represent ANPP of the year during and after the climate event, respectively.

Calculation of GWP of greenhouse gases

We calculated the GWP of three greenhouse gases, including CO2, CH4, and N2O. Soil N2O flux was measured continuously for 3 years (November 2017 to November 2020) using an automated soil N2O measurement system (UGGA, Los Gatos Research Inc., Mountain View, CA, USA) with eight chambers. Eight soil collars (diameter: 19.8 cm; height: 11.4 cm) were inserted 2 cm into the soil in each plot. Soil N2O flux of individual chambers was measured once every 2 hours. All living plants inside the collars were carefully clipped from the soil surface during the entire study period. The GWP index was defined to measure the time-integrated global mean radiative forcing of a pulse emission of a specific compound relative to that of CO2 (). According to the Intergovernmental Panel on Climate Change (IPCC AR6), the average GWPs for CO2, CH4, and N2O are 1, 28, and 273, respectively, in the time horizon of 100 years. In this study, we used the following equations to calculate GWPs for all greenhouse gases and their net GWPwhere FCO, FCH, and FN were annual fluxes of CO2, CH4, and N2O between terrestrial ecosystems and the atmosphere, respectively.

Data analysis

Repeated-measures ANOVAs were used to examine warming, year, and their interactive effects on soil temperature, soil salinity, primary productivity, and all carbon fluxes. One-way ANOVA was used to examine warming effects on primary productivity, species biomass, and carbon fluxes. Linear regressions were performed to explore the relationships between carbon fluxes and ANPP. The statistical analyses were conducted using SPSS 17.0 (SPSS for Windows, Chicago, IL, USA). We applied an SEM to explore the mechanisms of the warming effect on carbon fluxes. In the SEM analysis, we first considered a full model that included all possible pathways. Then, we sequentially eliminated insignificant pathways until we attained the final model (χ2 = 54.3, df = 39, P = 0.053; fig. S3). This analysis was conducted using IBM SPSS Amos 21.0 (SPSS for Windows, Chicago, IL, USA). Multiple linear regression was performed to explore the relationships between aboveground plant biomass and species traits (). A significance level of P = 0.05 was used for all statistical tests.