Plant biomass and soil organic carbon are main factors influencing dry-season ecosystem carbon rates in the coastal zone of the Yellow River Delta

Coastal wetlands are considered as a significant sink of global carbon due to their tremendous organic carbon storage. Coastal CO2 and CH4 flux rates play an important role in regulating atmospheric CO2 and CH4 concentrations. However, the relative contributions of vegetation, soil properties, and spatial structure on dry-season ecosystem carbon (C) rates (net ecosystem CO2 exchange, NEE; ecosystem respiration, ER; gross ecosystem productivity, GEP; and CH4) remain unclear at a regional scale. Here, we compared dry-season ecosystem C rates, plant, and soil properties across three vegetation types from 13 locations at a regional scale in the Yellow River Delta (YRD). The results showed that the Phragmites australis stand had the greatest NEE (-1365.4 μmol m-2 s-1), ER (660.2 μmol m-2 s-1), GEP (-2025.5 μmol m-2 s-1) and acted as a CH4 source (0.27 μmol m-2 s-1), whereas the Suaeda heteroptera and Tamarix chinensis stands uptook CH4 (-0.02 to -0.12 μmol m-2 s-1). Stepwise multiple regression analysis demonstrated that plant biomass was the main factor explaining all of the investigated carbon rates (GEP, ER, NEE, and CH4); while soil organic carbon was shown to be the most important for explaining the variability in the processes of carbon release to the atmosphere, i.e., ER and CH4. Variation partitioning results showed that vegetation and soil properties played equally important roles in shaping the pattern of C rates in the YRD. These results provide a better understanding of the link between ecosystem C rates and environmental drivers, and provide a framework to predict regional-scale ecosystem C fluxes under future climate change.

Introduction Carbon dioxide (CO 2 ) and methane (CH 4 ) are key greenhouse gases (GHGs) that make substantial contributions to global warming [1]. Numerous studies have estimated global wetland CO 2 and CH 4 fluxes, but with great uncertainties, mainly due to complicated environmental drivers [2][3][4]. Coastal wetlands have been recognized as the most vulnerable and sensitive ecosystems, because they act as the ecotone between terrestrial and aquatic ecosystems [5]. Coastal estuary wetlands store at least 430 Tg of carbon (C) with a C sequestration rate of 45 g C m -2 yr -1 , playing an important role in the global carbon cycle as natural carbon pools [6]. The coastal wetland is one of the most important wetland types for understanding C flux dynamics due to the high variations involved with water conditions, sedimentation characteristics, and vegetation types [7]. Coastal wetlands can act as greenhouse gas sinks via C burial, sediment deposition, and plant biomass accumulation, and as greenhouse gas sources through the release of CO 2 and CH 4 produced by the decomposition of organic matter [8], so they are of vital importance in governing the atmospheric concentrations of CO 2 and CH 4 [9] . However, due to the complicated interaction of environmental factors including vegetation and soil properties, how to disentangle the contributions of multiple drivers to CO 2 and CH 4 fluxes in estuary wetland remains unclear.
Vegetation exerts a major influence on C fluxes including net ecosystem exchange (NEE), ecosystem respiration (ER), gross ecosystem productivity (GEP), and CH 4 flux [10][11][12][13]. Niu et al. [14] found that a shift in the coverage of dominant plants could regulate the ecosystem C fluxes including NEE, ER, and GEP. A meta-analysis showed that NEE, ER, and GEP varied with different vegetation types in global coastal wetlands [12]. A previous study demonstrated that the plant biomass of Arctophila fulva was a strong predictor of C flux in an arctic tundra wetland [10]. Spartina alterniflora could alter the relationship between CH 4 and electron acceptors, resulting in an increase of CH 4 flux [11].
Soil properties have been found to be strongly associated with ecosystem C fluxes. Previous studies have reported that soil organic carbon (SOC), dissolved organic carbon (DOC), NH 4 + , pH, and salinity all exhibited dominant effects on ecosystem C fluxes [2,[15][16][17][18][19]. For example, Ardón et al. [20] identified salinity and hydrology as the most important determinants of GHG fluxes in salt marsh. However, other studies have considered DOC, SOC, and pH as the related factors to predict CH 4 flux in coastal wetlands [15,16]. Furthermore, the soil water level has also been considered as the main driving factor of regulating CO 2 and CH 4 fluxes in estuary and other coastal wetlands [18,19,21,22]. However, the underlying mechanism of regulating ecosystem C fluxes remains unclear.
Spatial structure, which represents underlying effects of the heterogeneity of environmental factors influences the pattern of C fluxes in a different way to biological and environmental factors acting on community and ecosystem [23]. Previous studies have found that soils with similar environmental characteristics have similar microbial communities, which is important for C release through C decomposition and CH 4 production [24,25]. It was found that spatial heterogeneity imposed strong influence on the variation of soil CO 2 efflux in tropical riparian ecosystems [4]. However, understandings of the role spatial structure plays in determining ecosystem C fluxes in coastal wetlands is still limited.
The Yellow River Delta (YRD) wetland, which is the largest wetland ecosystem in the warm temperate zone of China, has an obvious dry season (April to June) and rainy season (July to September) [26][27][28]. Previous studies have reported seasonal variations in ecosystem C fluxes (CO 2 and CH 4 ) based on continuous flux measurements [28][29][30] and found that dry-season had the second largest contribution to the CO 2 and CH 4 emissions in the YRD [29]. Therefore, in this study, we emphasized the relationship between dry-season ecosystem C rates (GEP, ER, NEE, and CH 4 ) with vegetation and soil properties and disentangled their contributions to the C rates. Specifically, we compared dry-season ecosystem C rates from 13 locations across three vegetation types at a regional scale in the YRD. The objectives of this study are (1) to test whether soil salinity was the primary determinant affecting ecosystem C rates in the YRD; and (2) to disentangle the contributions of vegetation, soil properties, and spatial structure on dryseason ecosystem C rates at a regional scale.

Site description
The coastal zone of the Yellow River Delta (Fig 1) has a temperate semi-humid continental monsoon climate [31]. The annual average temperature is 12.9˚C. The average annual precipitation is 530-630 mm and the rainfall mostly precipitates from July to September (rainy season). The dry season (April to June) accounts for about 30% of the annual precipitation [28]. The abundant vegetation of this research area includes Phragmites australis, Suaeda heteroptera Kitag, and Tamarix chinensis. Saline soil is the main soil type and the soil texture is sandy loam. Plant biomass and soil properties of each location are listed in Table 1. Locations of samples collected are listed in S1 Table. All necessary permits were obtained for the described field research. We are authorized by the Administrative Office of Yellow River Delta National Nature Reserve to carry out soil and plant collection.

Ecosystem carbon fluxes measurement
Previous studies have reported seasonal variation of ecosystem C fluxes (CO 2 and CH 4 ) based on continuous flux measurements in the Yellow River Delta [29,30], while this study emphasized the relationship between dry-season ecosystem C rates (GEP, ER, NEE, and CH 4 ) with plant and soil properties. Therefore, dry-season ecosystem C rates were determined in May 2018. Ecosystem C rates including GEP, ER, NEE, and CH 4 were measured by the closed chamber method. A 40 × 40 cm square stainless steel frame was inserted into the soil at each location. The CO 2 and CH 4 fluxes were measured using a greenhouse gas analyzer (UGGA, LGR, USA) attached to a transparent chamber (0.4 m × 0.4 m × 0.5 m). During the measurements, the chamber was placed on the frame and sealed with water [32], and two small fans ran continuously to mix the air inside the chamber. The relationship between CO 2 concentration and time was linear during the measurement interval. The CO 2 flux rates were then determined from the slope of the CO 2 concentration-time function. The coefficients of determination (r 2 ) of the linear function were greater than 0.95. Following the measurement of NEE, the chamber was vented for 1-2 min, replaced on the frame, and covered with an opaque cloth to measure ER. Negative NEE values represented net CO 2 uptake by the ecosystem, and positive NEE values represented net CO 2 loss. GEP was calculated as the difference between NEE and ER. The CO 2 and CH 4 flux measurements were performed between 9:00 am and 14:00 am. Air temperature was determined manually by inserting a digital thermometer onto the transparent chamber when ecosystem CO 2 fluxes were measured. Soil temperature was measured.by inserting a digital thermometer into soil at the depth of 0-5 cm.

Plant and soil sampling and element analyses
After the carbon flux measurements, the plant coverage was estimated using a 40 × 40 cm quadrat with a grid size of 4 × 4 cm in each plot, and plants were subsequently harvested for biomass determination. Five soil cores (3.5 cm diameter) were collected randomly from each plot at a 0-15 cm depth and mixed to one composite sample. The samples were passed through a 2 mm sieve and divided into two parts. One part of fresh soil was used for the determination of soil water content (SWC) and the analysis of ammonium (NH 4 + ), nitrate (NO 3 -), microbial biomass carbon (MBC), and dissolved organic carbon (DOC). The other part was air dried for the determination of soil pH, salinity, soil organic carbon (SOC), total phosphorus (TP), and available phosphorus (AP). Soil NH 4 + and NO 3 concentrations were determined by extraction with 2 M KCl solution followed by colorimetric analysis on a FIAstar 5000 Analyzer (FIAstar 5000 Analyzer, Foss Tecator, Hillerød, Denmark). Soil MBC was estimated by using a chloroform fumigation extraction method [33]. Soil DOC was extracted by adding 50 mL of 0.5 M potassium sulfate to subsamples of 12.5 g of homogenized soil, and by agitating on an orbital shaker at 120 rpm for 1 h. The filtrate was analyzed using a TOC analyzer (multi N/C 3100, Analytik Jena, Germany). The soil pH was determined in a 1:2.5 soil:water solution (w/v), and the soil conductivity was determined as an index of salinity. SOC and TN were analyzed using a C/N analyzer (multi-N/C 3100 Analytik Jena AG, Germany). TP was analyzed by applying the Murphy Riley method following perchloric acid digestion [34] and AP was determined by treatment with 0.5 mol L -1 NaHCO 3 followed by molybdenum blue colorimetry [35] using a spectrophotometer (UV2550, Shimadzu, Japan).

Statistical analyses
One-way ANOVA with the Duncan test was used to determine the differences of dry-season ecosystem C rates between vegetation types. Stepwise multiple regression was performed using step function based on the smallest AIC selection in in R v3.5.0 with the vegan package [36]. Unconstrained ordination (principal component analysis, PCA) was used to compare ecosystem carbon rates among plots (n = 39). Constrained ordination (redundancy analysis, RDA) was used to represent the relationships between environmental factors (vegetation, soil and space structure), plot patterns, and ecosystem carbon rates [37]. Vegetation type (S. heteroptera, P. australis, and T. chinensis) was considered as a qualitative factor and other environmental factors as quantitative factors.
In order to separate the effects of environmental factors on ecosystem carbon rates, the variation partitioning procedure was conducted. The environmental factors were divided into three groups: (1) vegetation factors including vegetation type, plant biomass, and plant cover; (2) soil factors including soil water content (SWC), soil pH, soil salinity, soil organic carbon (SOC), total phosphorus (TP), available phosphorus (AP), dissolved organic carbon (DOC), microbial biomass carbon (MBC), NH 4 + (ammonia), and NO 3 -(nitrate); (3) spatial structure (x, y, xy, x 2 , y 2 , x 2 y, xy 2 , x 3 , y 3 ), where nine terms of latitudinal (x) and longitudinal (y) coordinates were used to calculate a cubic trend surface. Spatial trend surface analysis is one of the quantitative ecological methods used to study the relation between spatial structure and community [23]. PCA, RDA, and variation partitioning analysis were performed in R v3.5.0 with the vegan package [36].

Variation in ecosystem carbon rates
The ecosystem carbon rates varied with different vegetation types (Fig 2). Specifically, stands of P. australis had the highest NEE, ER, and GEP with -1365.3, 660, and -2025.5 μmol m -2 s -1 , respectively, while S. heteroptera exhibited the lowest values with -272.3, 43.9, and -316.2 μmol m -2 s -1 , respectively. The difference of NEE between S. heteroptera and T. chinensis plots was not significant (P > 0.05, Fig 2A). Although the NEE of T. chinensis was significantly lower than that of P. australis, ER was almost the same as the latter (Fig 2B). Meanwhile, the GEP of P. australis was much greater than that of T. chinensis, which resulted in the highest value of NEE in P. australis (Fig 2A and 2C). With regard to CH 4 flux, the S. heteroptera and T. chinensis plots exhibited -0.019 and -0.12 μmol m -2 s -1 , respectively, and P. australis exhibited 0.27 μmol m -2 s -1 (Fig 2D). The first axis of PCA ordination explained 82.4% of the variation in the ecosystem carbon rates, mainly reflecting vegetation types (Fig 3A and 3B). P. australis and T. chinensis with higher NEE, ER, GEP, and CH 4 were plotted along the right side of the first axis. S. heteroptera with lower NEE, ER, GEP, and CH 4 were concentrated on the left side of the first axis. The second axis of PCA ordination explained 13.7% of the variation in the ecosystem carbon fluxes, mainly associated with soil properties and space variation. The positions of P. australis and T. chinensis were separated with the second axis (Fig 3A).

Relationship between ecosystem carbon rates and environmental factors
Ecosystem carbon rates across three vegetation types at the regional scale were distinguished by environmental factors with the RDA ordination (Fig 4A and 4B). The first axis described 75.5% of variation in the ecosystem carbon rates, mainly associated with SWC, pH, and salinity. The second axis explained 5.7% of the variation, primarily related to vegetation type, plant biomass, and plant coverage. Soil properties including MBC, DOC, and NO 3 did not show strong relationships, therefore were not considered as the major factors influencing ecosystem carbon rates (Fig 4B).
Stepwise multiple regression analysis demonstrated that 84% of the variation in NEE could be explained by DOC, plant coverage, and biomass together. SOC, plant coverage, and biomass together contributed for 90% of the spatial variation in ER. Plant coverage and biomass together explained 94% of the variation of GEP. Of the variation of CH 4 rate, 66% could be attributed to the combination of SOC, salinity and AP (Table 2).

Variation partitioning
Forward selection of the three groups of environmental factors with RDA suggested that the variation of ecosystem C fluxes was significantly associated with the vegetation (plant biomass) and soil properties (SOC and NH 4 + ). The variation partitioning results demonstrated that the total explained variation in ecosystem carbon fluxes was 82.8% and the undetermined variation was 17.2% (Fig 5). The vegetation, soil, and spatial factors explained a variation of 13.6%, 16.8%, and 11.3%, respectively, while the largest fraction in the determined variation was the combination of vegetation, soil, and space (21.6%). In addition, the overlap of vegetation and soil, vegetation and space, and soil and space explained a variation of 5.5%, 7.1%, and 6.9%, respectively (Fig 5).

Effects of vegetation on ecosystem C rates
In the exploration of the primary drivers regulating coastal wetland C rates (NEE, ER, GEP, and CH 4 ) and disentangling the relative contributions of multiple environmental factors (vegetation type, plant biomass, coverage, soil pH, salinity, SOC, SWC, TP, AP, DOC, MBC, NH 4 + and NO 3 -) and spatial structure on ecosystem C rates, in this study, plant biomass showed as  The mean value of NEE, ER, and GEP of the P. australis stand was the highest among the three species (Fig 2), suggesting that vegetation plays an important role in regulating ecosystem CO 2 exchange ( Table 2). The role of biotic control of plants on ecosystem CO 2 exchange has been evaluated in various Chinese coastal wetland ecosystems [6,30,38]. Higher ecosystem CO 2 exchange rates have been attributed to the higher plant biomass, as plant biomass can australis, and Tamarix chinensis), plant biomass and plant coverage. Soil factors include soil water content (SWC), soil pH, soil salinity, soil organic carbon (SOC), total phosphorus (TP), available phosphorus (AP), dissolved organic carbon (DOC), microbial biomass carbon (MBC), NH 4 + (ammonia), and NO 3 -(nitrate). Vegetation type was plotted as centroids (qualitative factor) and others were plotted as vectors (quantitative factors).
https://doi.org/10.1371/journal.pone.0210768.g004 Table 2. Results of stepwise multiple regression analysis. Independent variables: vegetation type, plant biomass, plant coverage, soil water content (SWC), soil pH, soil salinity, soil organic carbon (SOC), total phosphorus (TP), available phosphorus (AP), dissolved organic carbon (DOC), microbial biomass carbon (MBC), NH 4 + (ammonia), NO 3 -(nitrate), soil temperature at depth of 5 cm (Ts), and air temperature (Ta). The bold numbers represent significance (p < 0.05). regulate the supply of plant litter and hence organic matter decomposition, which in turn affects the rate of CO 2 production [39]. In this study, plant biomass and coverage were included in the regression models of NEE, ER, and GEP (Table 2). These results suggest that the variations in plant biomass could be one possible reason accounting for the difference in ecosystem CO 2 exchange among the three vegetation types in the YRD. It has been proposed that plant biomass could be used to indicate the levels of CH 4 flux because the litter input provides C resources for the growth of methanogenesis [10,40]. For instance, Andresen et al. [10] reported that CH 4 flux responded linearly to the biomass increase of Arctophila fulva. However, regression results suggested that plant biomass was not the main factor influencing CH 4 emission in this study ( Table 2). The difference of mean CH 4 rate between the three vegetation types was significant and P. australis had the highest CH 4 emissions with the others as a sink of CH 4 (Fig 2). Differences in the CH 4 rate among species were likely to be related to a combination of factors such as the growth form [41] and the depth preference among species, which may influence soil temperature [42], methanogenesis, and CH 4 transport [40]. The high CH 4 emission rate of P. australis might be attributed to several reasons. Firstly, large aerenchyma conduits and the presence of this species in relatively deeper water, resulting in the capability of P. australis roots to extend into CH 4 -rich soils without competition for methanogenic substrates and allowing CH 4 to be ventilated straight from the soil of these aquatic systems and into the atmosphere [43]. Secondly, high plant biomass and the perennial nature of P. australis could change soil micro-environment (e.g. higher soil moisture, anaerobic condition) and thus result in high CH 4 emission [40,43]. These findings suggest that although plant biomass is not the best predictable variable in this coastal wetland, plant species plays an important role in regulating the CH 4 rate on a community scale.

Effects of soil properties on ecosystem C rates
Microbiological processes and the roles of microorganisms could be another reason for the variation in ecosystem CO 2 and CH 4 rates. Their activities are controlled by several biological, chemical, and physical factors in soil [19,20,22,28,29,44]. Therefore, soil properties including soil organic carbon (SOC), total nitrogen (TN), inorganic nitrogen, total phosphorus (TP), available phosphorus (AP), salinity, and pH are closely related to CO 2 and CH 4 fluxes [45]. Previous studies have reported that the variation of ecosystem CO 2 exchange (NEE, ER and GEP) is associated with the shift in sediment pH and salinity [46]. This could be attributed to the high sensitivity of soil microbes to changes in pH and salinity [47]. Previous researches also demonstrated that ecosystem carbon fluxes were regulated by soil N and P concentration and nutrient stoichiometry [48,49]. The present study showed that NEE had a significant relationship with DOC, plant coverage, plant biomass, TP, AP, NO 3 -, MBC, and Ta; ER had a significant relationship with SOC, plant coverage, plant biomass, TP, and soil water content (SWC); and GEP was closely related with plant biomass, plant coverage, AP, and NO 3 -. It is worth noting that in contrast to previous studies, NEE, ER, and GEP were not significantly associated with soil pH and salinity in this coastal zone of the YRD.
Methanogenic bacteria are pH sensitive and most of them grew over a relatively narrow pH range of about 6-8, with an optimum pH of 7.7 for methane emission in coastal wetlands [50]. In this study, pH ranged from 8.0 to 9.0, which is not an optimum pH for methane emission. The CH 4 rate in coastal wetlands was enhanced at low salinity due to the intense oxidation or alleviation of competition by more efficient sulfate-and nitrate-reducing bacteria than methanogenic bacteria and a previous study has shown a strong negative relationship between CH 4 emission and salinity [45]. In this study, except for SOC and salinity, the CH 4 rate was significantly associated with AP, which was probably due to the intense anthropogenic nutrient inputs [51].

Future needed studies and implication for dry-season ecosystem C rates in coastal wetlands under climate change
It should be noted that only the dry-season ecosystem C rates (GEP, ER, NEE and CH 4 ) were measured in this study, which hampered calculating the yearly C fluxes. However, the conclusions drawn from the results of the relationship between dry-season ecosystem C rates and vegetation and soil properties were not compromised. Further studies should investigate the longterm ecosystem C rates and their relationships with vegetation and soil factors, which provide the foundation in which to explore the mechanism of the variation of ecosystem C fluxes.
In our study area, the vegetation was not homogeneous and the community types comprised of the native species P. australis, S. salsa, and T. chinensis. The spatial distribution pattern of vegetation was mostly identified as patches of these species or mud flats (bare soil) on the scale of meters to tens of meters due to the soil salt content gradient [6]. Wetland vegetation is the carrier of coastal C sink, which is of vital importance to mitigate climate change [5]. In addition, we did not sample the underground biomass of these stands to explain the variation of ecosystem C rates. These parameters warrant further investigation.
Previous study has demonstrated that the average annual precipitation has decreased with a rate of 4.5 mm yr -1 and the annual air temperature has increased by 1.7˚C over the past 55 years in the YRD [28]. This shift of precipitation and air temperature might induce drier and hotter seasons in the YRD in the future, which has a great possibility to increase ER, thus greater C loss from this coastal wetland. These findings have important implications for predicting ecosystem C rates to changes in precipitation and air temperature under climate change [52].

Conclusions
In conclusion, this study comprehensively investigated the regulation of ecosystem C flux rates under vegetation and soil property gradients in the coastal zone of the Yellow River Delta. This study showed that a combination of biotic and abiotic predictors, e.g., vegetation and soil, explained the majority of the variation in the investigated dry-season ecosystem C rates in the coastal zone of the Yellow River Delta. Plant biomass was found to be the main factor explaining all of the investigated carbon rates (GEP, ER, NEE, and CH 4 ), while soil organic carbon was shown to be the most important for explaining the variability in the processes of carbon release to the atmosphere, i.e., ER and CH 4 . Vegetation and soil properties played equally important roles in shaping the pattern of C rates through variation partition analysis. The results of this research provide a better understanding of the link between ecosystem C rates and environmental drivers, and provide a good basis to predict regional-scale ecosystem C fluxes under future climate change.
Supporting information S1