Freshwater Releases Into Estuarine Wetlands Change the Determinants of Benthic Invertebrate Metacommunity Structure

In recent years, the relative importance of the processes driving metacommunity composition has aroused extensive attention and become a powerful approach to identify community patterns and their regulatory mechanisms. We investigated variations in the composition of benthic community in restored wetlands and natural wetlands in the Yellow River Delta (Shandong Province, China). First, spatial structures within each wetland were modeled with Moran eigenvector maps. Next, the variation in community structure among local environmental and spatial variables was partitioned using constrained ordination, and the “elements of metacommunity structure” analysis was used to determine the patterns of best fit for species distributions within metacommunities. Finally, the null model was used to analyze non-random patterns of species co-occurrence. The community structure of benthic invertebrates in restored wetlands and natural wetlands differed significantly. The benthic invertebrate metacommunity structure showed a nested distribution in restored wetlands and a quasi-Clementsian structure in natural wetlands. Pure environmental fractions and pure spatial fractions were critical in regulating benthic invertebrate metacommunities of restored wetlands. In natural wetlands, pure spatial fractions and the interaction between environmental and spatial factors (shared fractions) played a major role in the metacommunity. A species co-occurrence analysis showed that species co-occurred more frequently than expected by chance, demonstrating that biotic interactions were not the main driver of metacommunity structures in both wetland types. Accordingly, the benthic invertebrate metacommunity in estuarine wetlands following freshwater releases was mostly determined by environmental and spatial effects, which resulted in a metacommunity with nested distribution. These results are important for biodiversity protection and ecosystem management of estuarine wetlands in the Yellow River Delta.

In recent years, the relative importance of the processes driving metacommunity composition has aroused extensive attention and become a powerful approach to identify community patterns and their regulatory mechanisms. We investigated variations in the composition of benthic community in restored wetlands and natural wetlands in the Yellow River Delta (Shandong Province, China). First, spatial structures within each wetland were modeled with Moran eigenvector maps. Next, the variation in community structure among local environmental and spatial variables was partitioned using constrained ordination, and the "elements of metacommunity structure" analysis was used to determine the patterns of best fit for species distributions within metacommunities. Finally, the null model was used to analyze non-random patterns of species co-occurrence. The community structure of benthic invertebrates in restored wetlands and natural wetlands differed significantly. The benthic invertebrate metacommunity structure showed a nested distribution in restored wetlands and a quasi-Clementsian structure in natural wetlands. Pure environmental fractions and pure spatial fractions were critical in regulating benthic invertebrate metacommunities of restored wetlands. In natural wetlands, pure spatial fractions and the interaction between environmental and spatial factors (shared fractions) played a major role in the metacommunity. A species co-occurrence analysis showed that species co-occurred more frequently than expected by chance, demonstrating that biotic interactions were not the main driver of metacommunity structures in both wetland types. Accordingly, the benthic invertebrate metacommunity in estuarine wetlands following freshwater releases was mostly determined by environmental and spatial effects, which resulted in a metacommunity with nested distribution. These results are important for biodiversity protection and ecosystem management of estuarine wetlands in the Yellow River Delta.

INTRODUCTION
The mechanisms underlying both the patterns and maintenance of biodiversity within communities have been a core research topic in community ecology (Rosenzweig, 1995). Understanding the relative importance of processes driving the composition of community is the main way to reveal the drivers of community assembly (Bell, 2010;Caruso et al., 2011). According to niche and neutral theory, community assembly depends on abiotic factors, biotic interactions, priority effects (i.e., competitive dominance due to early colonization), and dispersal processes (Leibold et al., 2004). Discerning the relative importance of these factors has been a major challenge for understanding the assembly of ecological communities (Holyoak et al., 2005;Hildrew, 2009;Heino et al., 2015a). However, the relative importance of these factors varies with biota and spatiotemporal scales (Qian and Ricklefs, 2007;Qian, 2009;Logue et al., 2011;Heino et al., 2012). The degree to which the influence of these factors when the focus is on local-scale communities remains crucial . Therefore, knowledge on the determinants of local taxonomic diversity is the key to understanding ecological community composition (Daniel et al., 2019) and protecting biodiversity (Heino, 2013c).
In terrestrial, freshwater, and marine systems, the effects of spatial and environmental gradients on community composition have been widely studied (Soininen et al., 2007;Heino and Tolonen, 2017). Environmental gradients are closely related to species distribution in freshwater ecosystems (Heino and Soininen, 2010;Heino, 2011;Saito et al., 2015). For example, previous studies have shown that environmental factors, including lake size, depth, pH, nutrients, and biotic interactions, such as predation, are the main determinants of species distribution patterns in lakes . Due to dispersal barriers and isolation effects, the spatial positioning of lakes affects species colonization patterns. Species distributions are also affected by dispersal and other stochastic forces (Tonkin et al., 2015). Considerable variability amongst regions and benthic invertebrate groups has been reported, and environmental variables and species classification appear to be more important than spatial variables in driving the structure of ecological communities (Cottenie, 2005;Thornhill et al., 2017). Additionally, increasing evidence supports the importance of species sorting signals along geographical and environmental gradients in aquatic ecosystems (Heino, 2013a;Alahuhta et al., 2018) and emphasizes the influence of biotic interactions on biological community (García-Girón et al., 2020). Biotic interactions were the primary driving force of beta diversities of multiple organismal groups in ponds and, to a lesser extent, the abiotic environment (García-Girón et al., 2020). Hence, more biological groups and different aquatic ecosystems should be examined to gain a more general picture of the patterns of community organization at local scales (Cai et al., 2017).
The metacommunity concept has the potential to integrate local and regional dynamics within a general community ecology framework, providing a valuable analytical opportunity to identify community patterns and their regulatory mechanisms (Thompson et al., 2020). It refers to communities of interacting species connected by dispersal and responding to local environmental conditions (Leibold and Chase, 2017). Metacommunity theory integrates the abiotic environment, biotic interactions, and dispersal events as drivers of biological diversity across spatial scales. Assessing the relative importance of these drivers is a core objective of metacommunity ecology (Meynard et al., 2013;Brown et al., 2017). During the development of the metacommunity framework, six species distribution patterns (i.e., the nested subset, chessboard, Clementine gradient, Gleason gradient, uniform interval gradient, and random pattern) were identified, and each pattern was regulated by a specific mechanism (Leibold and Mikkelson, 2002;Leibold et al., 2004). Therefore, the metacommunity approach provides a research framework for identifying the mechanisms underlying species co-existence and community assembly at local scales (Leibold et al., 2004;Hill et al., 2017;Guo et al., 2019;Liu et al., 2019).
Different aquatic ecosystems exhibit fundamental differences in environmental heterogeneity, connectivity, and spatial extension, which provides challenges and opportunities for examining the organization of metacommunities at local scales (Brown et al., 2011;Heino, 2011;Lindstrom and Langenheder, 2012;Heino et al., 2015a). Estuarine wetlands which are transitional areas between land and sea are affected by the interactions among multiple systems including rivers, oceans, wetlands, and upland catchments (Yang et al., 2017a. Because of extensive interactions between freshwater and saltwater, estuarine wetlands form unique ecological patterns and exhibit biogeochemical cycles and biological communities different from lakes, rivers, and oceans. In estuarine wetlands, benthic invertebrates are directly involved in the material circulation and energy flow of ecosystems, and are sensitive to environmental changes and disturbances. Therefore, they are often used as biological indicator for bioassessments Wu et al., 2017Wu et al., , 2019Lu et al., 2019;Guan and Wu, 2021). In recent years, due to the multiple influences of natural factors and human activities, many estuarine wetlands have been degraded, their biodiversity has been diminished, and their ecosystem structure and function have been compromised. To prevent further degradation of the ecosystems, the management of the reserve has carried out ecological restoration projects for the degraded wetlands in the Yellow River Delta. However, restoration measures included the construction of levees for impoundments which blocked the connectivity between the restored wetlands and the sea. The diversity of benthic invertebrates might respond to the resulting changes in the physical and chemical properties of the water body, and the driving factors of community structure in restored wetland may differ from those in wetlands directly connected to the sea. Therefore, it is important to evaluate the benthic invertebrate metacommunity patterns and drivers after long-term freshwater releases, which can provide basic data and theoretical support for biodiversity conservation and restoration, as well as ecosystem management of the estuarine wetlands of the Yellow River Delta.
Most previous studies have focused primarily on community assembly mechanisms in freshwater ecosystems such as lakes, streams, and ponds, and few empirical studies have been performed in estuarine wetlands using similar approaches as in freshwater systems (Tonkin et al., 2018). In the present study, we assessed the responses of benthic invertebrate metacommunities to ecological restoration utilizing freshwater releases and compared them with nearby natural wetlands, which have not received freshwater releases (i.e., reference areas).
Our main goals were to: (1) assess the differences in benthic invertebrate community structure between restored and natural wetlands; and (2) identify the major factors driving benthic invertebrate metacommunity structure in the two types of wetlands. Our findings may help wetland managers to assess the ecological effects of freshwater releases, provide valuable insights for enhancing the restoration of estuarine wetland biodiversity, and improve the management of wetland ecosystems.

Study Site
The Yellow River Delta estuarine wetland is located within Shandong Province, China (37 • 34.768 ∼ 38 • 12.310 N, 118 • 32.981 ∼ 119 • 20.450 E), between Bohai and Laizhoubays (Figure 1). This area is a typical fan-shaped delta encompassing a binary phase structure of river sediments covering a marine layer. As this wetland is located at mid-latitudes, it has a warm-temperate semi-humid continental monsoonal climate. The study area is characterized by annual average temperatures between 11.7 and 12.6 • C, including an extreme maximum temperature of 41.9 • C, an extreme minimum temperature of -23.3 • C, a frost-free period of 211 days, and annual precipitation between 530 and 630 mm.
Since 2010, the Yellow River Delta Nature Reserve Administration has carried out a large-scale restoration project in degraded estuarine wetlands, reshaping the wetlands' hydrological situation. The restoration measures include the construction of tidal barriers and water diversion canals and the implementation of freshwater releases through the Qingshuigou flow path of the Yellow River from June to July every year, normally for 20-30 days (Li et al., 2016;Yang et al., 2017b). After 5 years of freshwater releases, the landscape of the restored area has changed dramatically (Dong et al., 2014;Yang et al., 2017a). This project has successfully expanded the area of reed wetlands, improved the water and salt conditions, and increased the biodiversity of the Yellow River Delta Wetland .
The barriers between the restored and natural wetland are earth embankments, which limits the extent of freshwater releases and creates different habitat features on either side, while preventing the entry of seawater during high tides and storms. Sampling sites in restored wetlands were located at least 100 m from the barriers to ensure no connectivity with natural wetlands (Figure 1).

Sample Processing and Analysis
In October 2017 and May 2018, benthic invertebrates were collected from 32 sampling sites (19 restored wetlands and 13 natural wetlands) in the Yellow River Delta Wetland. Benthic invertebrates were collected at each sampling site using a 1 mm mesh D-shaped sweep net with a diameter of 35 cm. This net was dragged horizontally from the inside of the water body to the shore and scraped the benthic surface . Our sampling plan captured different subhabitats. Four sweeps were performed in each site and mixed into a composite sample, stored in a numbered plastic bag, preserved in 95% alcohol, and returned to the laboratory for identification and classification. Benthic invertebrates were separated from associated materials and identified to the genus or species level (Morse et al., 1994).
A YSI multi probe water quality system (556MPS; Yellow Springs Instruments, United States) was used to measure the physicochemical properties of the water of the two types of wetlands. Measurements included pH, electrical conductivity (Cond), dissolved oxygen (DO), salinity, total dissolved solids (TDS), oxidation-reduction potential (ORP), chloride ions (Cl − ), and chlorophyll-a (Chl-a). In addition, total nitrogen (TN), total phosphorus (TP), total suspended solids (TSS), total organic carbon (TOC), total carbon (TC), inorganic carbon (IC), and dissolved organic carbon (DOC), as well as nitrate (NO 3 − ), ammonium (NH 4 + ), sulfate (SO 4 − ), carbonate (CO 3 − ), and bicarbonate ion (HCO 3 − ) concentrations were analyzed in water samples following the Surface Water Environmental Quality Standards of China (GB3838-2002). Some of the measured water quality variables showed significant differences between the two types of wetlands (Figure 2, p < 0.05).

Benthic Invertebrate Community Structures and Water Quality
The Shannon-Wiener diversity index (H ), the Simpson index (λ), and the Number of species (S) were calculated to compare benthic invertebrate community structures in both restored and natural wetlands. First, all environmental variables and biodiversity indices were tested for normal distribution using the Kolmogorov-Smirnov test prior to analyses. Then, independent sample t-tests were used to analyze differences in water quality variables and biodiversity indices between restored and natural wetlands. The biodiversity indices were calculated with the "Picante" package in R (R Core Team, 2020). The Kolmogorov-Smirnov test and the independent samples t-test were calculated by SPSS statistics version 21.
Non-metric multidimensional scaling (NMDS) was used to identify differences in the spatial distribution of benthic invertebrate metacommunities structure between different wetlands based on species abundance data of each sampling site. Differences were tested by analysis of similarities (ANOSIM). NMDS and ANOSIM were calculated using the "metaMDS" and "anosim" functions of the "vegan" package in R.

Statistical Analysis of Environmental and Spatial Data
To determine the key environmental variables affecting the benthic invertebrate community, detrended correspondence analysis (DCA) based on benthic invertebrate abundances was used to test for linearity, and axis lengths were all greater than 4. Therefore, canonical correspondence analysis (CCA) was selected for subsequent analysis. The data of the benthic invertebrate metacommunities was optimized by eliminating rare species with a frequency of occurrence less than or equal to 2. The optimized species data and environmental data (except pH) were log 10 (x+1) transformed prior to CCA. Then, a forward selection process was performed based on the adjusted R 2 to choose significant environmental variables  FIGURE 1 | Sampling sites in restored and natural wetlands in the Yellow River Delta.
FIGURE 2 | Significant differences in water quality variables between restored and natural wetlands. Blanchet et al., 2008). DCA and CCA were performed using the "vegan" package in R.
Moran's eigenvector maps (MEM) were used to model spatial structures among sites within the estuarine wetlands derived from the geographic coordinates of the study sites . The MEM analysis produces a set of orthogonal spatial variables, representing the spatial variation across a range of spatial scales Declerck et al., 2011). These variables can be used as explanatory variables in direct gradient analysis to describe spatial patterns Frontiers in Ecology and Evolution | www.frontiersin.org in communities (Borcard and Legendre, 2002;Dray et al., 2012;. In our study, MEM spatial eigen functions were computed using the "PCNM" function of the "PCNM" package in R .

Variation Partitioning
To explain the relative contribution of environmental filtering and the effects of the spatial factors, variation partitioning was applied on redundancy analysis models (RDA) for the significant environmental variables and spatial variables.  , 2006). We reported adjusted R 2 values for all analyses because they are unbiased estimates of the explained variation . Variation partitioning analysis was implemented using the "varpart" function of the "vegan" package in R. The significance of explained variables was tested using the function "anova" in the package "vegan" in R.

Elements of Metacommunity Structure Analysis (EMS)
The EMS can be used to determine the patterns of the benthic invertebrate metacommunities in restored and natural wetlands. Here, the analytical approach followed previous studies by Heino et al. (Heino et al., 2015a(Heino et al., ,b, 2017, and results were interpreted according to Presley et al. (2010). Prior to analysis, the reciprocal average (RA) method was used to rearrange the species presenceabsence matrix (sites by species). The RA method generates one or more RA axes by maximizing the maximum correspondence between species scores and site scores (Leibold and Mikkelson, 2002). Our analysis was performed only for the first axis (Gao et al., 2016). Based on the metrics: (1) coherence, (2) species turnover, and (3) boundary clumping, EMS was used to identify the structure of the metacommunity (e.g., random, checkerboard, nestedness, evenly spaced, Clementian, Gleasonian, and quasi structural patterns; Presley et al., 2010;Heino et al., 2015b).
(1) Coherence was evaluated by comparing the observed number of embedded absences (EAbs) in the ordinated matrix to a null distribution of embedded absences from simulated matrices (Leibold and Mikkelson, 2002;Presley et al., 2010). The results of coherence were characterized as: (i) non-significant coherence, i.e., with a random structure, (ii) significantly negative coherence, i.e., EAbs are significantly higher than expected by chance, which is a checkerboard pattern, reflecting competitive exclusion, or (iii) significantly positive coherence, i.e., EAbs are significantly lower than expected by chance, indicating that communities are structured along environmental gradients. Significantly positive coherence refers to nestedness, evenly spaced gradients, and Gleasonian structure or Clementsian structure, which needs to be further analyzed and confirmed by species turnover and boundary clumping (Leibold and Mikkelson, 2002).
(2) Species turnover was measured as the number of times one species replaces (Rep) another from site to site in an ordinated matrix (Presley et al., 2010). The results of species turnover were characterized as: (i) non-significant, i.e., quasi-structure, (ii) significantly negative, i.e., Rep is significantly lower than expected by chance, referring to a nested pattern, or (iii) significantly positive, i.e., Rep is significantly higher than expected by chance, referring to evenly spaced, Gleasonian or Clementsian structures (Presley et al., 2010). The evenly spaced Gleasonian and Clementsian metacommunity structures can be separated based on boundary clumping (Leibold and Mikkelson, 2002). Furthermore, non-significant negative species turnover could indicate quasi-nestedness, and non-significant positive species turnover could indicate quasi-Gleasonian, quasi-Clementsian, or quasi-evenly spaced gradients. Similarly, non-significant positive turnover can be separated by boundary clumping (for details, see Presley et al., 2010).
(3) Boundary clumping was assessed using Morisita's index (I) and a subsequent Chi-squared test comparing observed and expected distributions of range boundary locations. The results of boundary clumping were characterized as: (i) Gleasonian structure, i.e., I-values were close to 1, (ii) Clementsian structure, i.e., I-values were significantly greater than 1, or (iii) evenly spaced distributions, i.e., I-values were significantly greater than 1. A p-level of 0.05 was selected to test the statistical significance in all analyses. All EMS analyses were performed using the "metacom" package in R (Dallas, 2013).

Co-occurrence Analysis
The patterns of species co-occurrence in the benthic invertebrates metacommunity were analyzed using the C-score (Gotelli, 2000;Gotelli and Ulrich, 2012). The C-score (Checkerboard score, Stone and Roberts, 1990) is based on the average co-occurrence rate among all possible pairs of species in a presence-absence matrix. The checkerboard pattern of benthic invertebrates reflects whether competitive interactions among species occur (Gotelli and Mccabe, 2002). Therefore, indices of the C-score have a good power to examine species co-existence patterns at local scales (Gotelli, 2000;Gotelli and Ulrich, 2012). The calculation method of C-score was as: where R is the row totals of sites containing both species, R i and R j are the matrix row totals for species i and j, and S is the number of sites in which both species occur (Gotelli and Mccabe, 2002).
If the presence-absence data matrices have a significantly higher C-score than randomly generated matrices, then a substantial number of species pairs co-occur (simulated value) less often than by chance (observed value), suggesting that spatial distributions may be structured by interspecific competition. In contrast, low C-scores suggest that several species co-occur more frequently than expected by chance, and species can aggregate in the same metacommunity due to environmental filtering (Gan et al., 2019;Guo et al., 2019).
To compare the significance of results across studies, standardized effect sizes (SES) should be calculated for each matrix. The standardized effect size measures the number of standard deviations that the observed index is above or below the mean index of the simulated communities (Gotelli and Mccabe, 2002). Based on the C-scores from the observed and the simulated communities, the standard effect sizes were calculated as: Assuming the SES is normally distributed, the 95% confidence interval of the SES ranges between -2.0 and 2.0 (Gao et al., 2016). C-score analysis was performed using the "EcoSimR" package in R software.

Benthic Invertebrate Community Structure
In total, 96 taxa were recorded at the 32 sites, including 3 phyla, 7 classes, 20 orders, and 59 families. Among these, there were 6 polychaetes, 2 oligochaetes, 58 aquatic insects, 7 bivalves, 13 gastropods, 3 crustaceans, and 7 malacostracans. Aquatic insects were the main components of restored wetlands, accounting for 74.36% of the total number of taxa. However, in natural wetlands, the benthic invertebrate community was mainly composed of annelids, bivalves, malacostracans, and gastropods (Supplementary Table 1).
The results of the NMDS (stress = 0.08) analysis showed that the benthic invertebrate community structure of restored and natural wetlands were clearly divided into two groups (Figure 3). Furthermore, there were significant differences between these two groups (R = 0.68, p = 0.001, Figure 2). The Shannon-Wiener index (H ), the Simpson index (λ), and the Number of species (S) of benthic invertebrates in the restored wetlands of the Yellow River Delta were significantly higher than those of natural wetlands (p < 0.05). Diversity indices of benthic invertebrates differed significantly between the two types of wetlands (Figure 4).

Metacommunity Structure Analysis in Restored Wetlands
According to the results of the EMS analysis, benthic invertebrate metacommunities showed a significantly positive coherence (i.e., the Abs was significantly lower than expected by chance), indicating that communities were structured along environmental gradients. To distinguish the structures of benthic invertebrate metacommunities, species range turnover and boundary clumping needed to be further evaluated. The results of species turnover showed that the benthic invertebrate metacommunities reflected a nested distribution pattern (i.e., lower turnover than expected by chance, Table 1 and Supplementary Figure 1).
The environmental variables selected in the CCA models, i.e., the main determinants of community structure in restored wetlands, were Cond, Salinity, and ORP. We used the selected main environmental variables and the significant spatial MEM variables (MEM2, MEM4) for further analyses ( Table 2). According to variance partitioning results, environmental and spatial variables in combination explained 43% of the interpretation rates [(S), (E), and (S+E), Table 2]. The amount of variance explained by the pure spatial fraction [27.2%, (E)] and the pure environmental fraction [11%, (S)] were significant (p < 0.005, Table 2), indicating that environmental and spatial factors played an important role for benthic invertebrates ( Table 2). The interaction of environmental and spatial factors [4.8%, (E+S)] explained little variation in community composition.
Furthermore, most of the 19 sites in restored wetlands exhibited non-random patterns of species co-existence within sites (Figure 5). The SES < -2 indicated an aggregated distribution of benthic invertebrate assemblages. Finally, the species co-occurrence pattern analysis clearly demonstrated that biotic interactions were not the main driver of the metacommunity structure.

Metacommunity Structure Analysis in Natural Wetlands
Based on the results of EMS, benthic invertebrate metacommunities showed positive coherent metacommunity structures, non-significant positive species turnover, and clear boundary clumping. Hence, the benthic invertebrate metacommunities of natural wetlands corresponded best with the Quasi-Clementsian structure ( Table 1 and Supplementary Figure 1).
After the forward selection procedure, the CCA and Moran's eigenvector maps identified two statistically significant environmental variables (pH and IC) and one significant spatial variable (MEM1, Table 2) for variance partitioning analysis. According to the results of variance partitioning analysis, 50.6% of the total variation in the benthic invertebrate community composition was explained by the spatial and environmental factors in combination [(S), (E), and (S+E); Table 2]. The    amount of variance accounted for by the pure spatial fraction [1.77%, (S)] was significant, but that of the pure environmental fraction [0.73%, (E)] was not. In addition, the amount of variance explained by the shared fraction, which reflected the interaction of environmental and spatial factors [25.6%, (E+S) in Table 2], was higher than that of other fractions, indicating that the interaction of environmental and spatial variables was the main driver for benthic invertebrate community composition of natural wetlands. Furthermore, frequency histograms of standardized effect sizes showed an aggregated distribution among the 13 benthic invertebrate metacommunities in natural wetlands, indicating that biotic interactions were not the main driver of metacommunity structure (Figure 5).

DISCUSSION
We examined three sources of variation in community structure, i.e., spatial factors, environmental variables, and biotic interactions. Our results indicated that the benthic invertebrate community structure was significantly altered due to changes in water quality as a result of freshwater releases, and the underlying processes of community composition were also changed. The benthic invertebrate metacommunity structure showed a nested distribution in restored wetlands and a quasi-Clementsian structure in natural wetlands. Further analyses suggested that possible biotic interactions were not the key drivers of metacommunity structure at the local scale in estuarine wetlands. Environmental filtering and dispersal appeared to be key drivers of metacommunity structure in restored wetlands, and dispersal played a major role in regulating benthic invertebrate communities of natural wetlands. At the local scale of estuarine wetlands, however, possible biotic interactions seemed not to play a role in metacommunity structure.
The environmental heterogeneity and water quality changes caused by restoration-related freshwater releases led to changes in benthic invertebrate community composition. Our results showed the benthic invertebrate community structure of restored wetlands and natural wetlands were significantly different. Yang et al. (2019) reported that the composition of macrobenthos functional groups tended to become more diverse after freshwater releases, which is consistent with our findings. Since the implementation of the ecological restoration project of degraded wetlands in the Yellow River Delta, freshwater releases have alleviated the increasing water and salt stress of degraded wetlands. At the same time, freshwater releases had a positive impact on water quality, soil organic matter content, plant communities, and bird communities in degraded wetlands, as well as direct or indirect effects on the species richness and diversity of benthic invertebrate communities (Cui et al., 2009). In another study, we found that salinity is an important indicator to distinguish hydrologic characteristics between the two types of wetlands . Due to the construction of levees around the degraded wetlands for impoundment, the hydrological connectivity between the freshwater release areas and the ocean has been cut off, and the habitat environment has been greatly changed. Benthic invertebrates are sensitive to environmental changes and disturbances, ultimately leading to significant differences in community composition between the two types of wetlands. Given these potential effects on benthic invertebrate taxa and diversity in both wetland types, there is a need to identify the major drivers of community assembly in these estuarine wetlands.
We found a nested distribution in restored wetlands. It has been shown that nested distribution is not necessarily a result of biotic interactions (Heino, 2009) and may be caused by community size (Verbruggen et al., 2012), spatial isolation (Butaye et al., 2001), and environmental filtering (Heino et al., 2015b). As a result of habitat changes in restored wetlands, sensitive species may disappear along a gradient, and this environmental filtering can cause community nestedness. Indeed, Verbruggen et al. (2012) also confirmed that the varying tolerance levels of taxa for relevant environmental factors might be responsible for the observed nestedness. Therefore, the relative influence of environmental and spatial variables also depends on species traits. Some species exhibit wide environmental tolerances, while others have very specific habitat requirements (Rodil et al., 2017). The presence of suitable habitats and the ability of species to reach these habitats is what affects the extent of species distribution. Explaining differences in species life histories can provide a better understanding of species distribution and coexistence patterns (Corte et al., 2018). In addition, dispersal limitation causes niches to remain vacant, thereby obscuring expected species-environment associations under environmental filtering and potentially resulting in community nestedness (Lekberg and Koide, 2005). We found a quasi-Clementsian structure for the benthic invertebrate metacommunities in natural wetlands. Due to the inherently highly heterogeneous natural wetland systems, species responses to environmental gradients are typically more complex than a simple gain or loss of species along ecological gradients, making Clementsian gradients a frequent occurrence (Ward et al., 2002;Heino, 2013b). Moreover, natural estuarine wetlands are typically connected to the sea, and the length of the environmental gradient increases with spatial extent, resulting in stronger correlations between species composition and the environment among sites, ultimately leading to a Clementsian structure (Soininen, 2014;. Variance partitioning allowed us to evaluate the independent contributions of spatial and environmental variables to the observed benthic invertebrate distribution in two wetland types. In restored wetlands, our analyses showed that both environmental and spatial variables explained substantial portions of the variability in benthic invertebrate community composition, although environmental variables had a slightly stronger influence. The Yellow River Delta ecological restoration and protection project blocked the connectivity of the restored wetlands to the sea, and spatial factors, therefore, restricted the dispersal of benthic invertebrate metacommunities. Moreover, habitat changes resulted in higher sensitivity of benthic invertebrates to environmental heterogeneity. Thus environmental filtering played a leading role in determining species composition (Rádková et al., 2015). In natural wetlands, spatial variables played a major role in regulating benthic invertebrate metacommunities, and the interaction between environmental factors and spatial variables was also crucial in the metacommunity. It has been shown that less dispersive species are more strongly affected by spatial processes than highly dispersive species (Rao et al., 2020). Additionally, the mechanisms determining community structure changed according to the spatial scale considered (Rao et al., 2020). The results suggested that the spatial component plays a greater role in metacommunity organization in open estuarine waters (e.g., shallow beaches) than in less open, environmentally controlled aquatic systems (Rodil et al., 2017).
Biotic interactions must be assessed to accurately understand metacommunity organization at local scales (García-Girón et al., 2020). Our results showed that biotic interactions appeared to have little effect on the structure of benthic invertebrate metacommunities of restored and natural wetlands. This may be due to the fact that resource-based niche partitioning and interspecific competition play minor roles in influencing benthic invertebrate communities at small scales (Gao et al., 2016;Guo et al., 2019). Moreover, anthropogenic disturbances and management measures may also lead to low overall levels of biotic interactions (Ochoa-Hueso, 2016). However, in contrast to our results, Zhao et al. (2019) found that biotic interactions significantly influenced the patterns of biomass, species richness, and community composition of bacteria, diatoms, and shaker mosquitoes along a water depth gradient. This may have been due to the fact that this study only used biotic predictors as surrogates for local biotic constraints and that the role of spatial factors was not considered. However, the exact reasons for these patterns are still not completely understood. Therefore, additional measurements and variables (e.g., food web structure and energy flow, functional properties, and phylogeny) should be considered to analyze the relative contributions of biotic interactions properly.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.