Ecological carrying capacity and carbon sequestration potential of bivalve shellfish in marine ranching: A case study in Bohai Bay, China

Introduction Shellfish play an important role in ecological restoration and as carbon (C) sinks, but studies on their ecological carrying capacity (ECC) and C sequestration potential are sparse. Methods In this study, we selected a 57-hectare artificial oyster reef in a typical marine ranching in Bohai Bay, China, to evaluate the ECC and their C sequestration potential of bivalve shellfish, and projecting their impact on functional groups in the system, with an Ecopath with Ecosim (EwE) food web model. We conducted four biological surveys to obtain the biomass measurements, with one conducted in each of the summer, autumn, and winter of 2019 and one in the spring of 2020; and the functional groups included in the surveys comprised fish, cephalopods, crustaceans, snails, bivalve shellfish, annelids, other macrobenthos, meiobenthos, starfish, sea cucumbers, zooplankton, phytoplankton, and detritus. Results and Discussion The EwE model prediction results showed that the ECC of bivalve shellfish was established to be 282.66 t/km2, far more than the existing quantity of 187.76 t/km2. Therefore, at present, the ecosystem of the study marine ranching is not yet mature. Moreover, our ecological network analysis parameters indicated that the marine ranching ecosystem will be mature and stable when the bivalve shellfish population reaches its ECC. However, the increase in bivalve shellfish biomass will result in a decrease in the population sizes of species competing for food resources with bivalve shellfish, mainly gobiid fish such as Tridentiger bifasciatus, Tridentiger trigonocephalus, Tridentiger barbatus. Simultaneously, when the bivalve shellfish reach their ECC, 29.23 t of CO2 can be sequestrated by bivalve shellfish, comprising 14.32 t being removed from the ecosystem as prey and 14.91 t being stored on the seafloor through biodeposition. Conclusion Therefore, the research demonstrated that, within the scope of ECC, the increasing bivalve shellfish can improve the C sequestration capacity of the marine ranch ecosystem, and effective management of bivalve shellfish in marine ranching can improve the economic benefits and C sink service functions of marine ranching.


Introduction
Increasing global greenhouse gas concentrations result in a series of environmental problems, such as global warming, ocean acidification, and extreme weather events (Pachauri et al., 2014;Yu et al., 2021). Coastal wetlands display a powerful carbon sink function, playing an important role in reducing carbon dioxide (CO 2 ) concentration and mitigating global climate change (Duarte et al., 2013). The carbon sequestered by coastal wetlands and seagrass beds is called coastal blue carbon (Zhao and Hu, 2019). During the coastal wetlands, marine ranching is a special offshore ecosystem (Dutta et al., 2019) that can sequester significant carbon (C) and mitigate the greenhouse effect (Shen et al., 2020). Bivalve shellfish is a key component in marine ranching, which represents important marine biological C sink due to high biodeposition, calcification rates (Pan et al., 2019a;Pan et al., 2019b). For example, oyster reefs exhibiting a strong C sequestration capability and serve as remarkable C reservoir (Grabowski and Peterson, 2007). However, the high population density of bivalve shellfish has negative-impacts on marine ecosystems, such as depressing the growth of phytoplankton and deteriorating bottom environments (Velasco et al., 2009). Therefore, it is critical to understand not only the C sequestration potential but also the ecological carrying capacity (ECC) of bivalve shellfish.
Carrying capacity is defined as the maximum biomass of a specific population that the ecosystem can support under specific environmental conditions, such as food and habitat availability (Filgueira et al., 2021). Inglis et al. (2000) classified the carrying capacity of shellfish into four categories: physical, productive, ecological, and social carrying capacity. ECC is the maximum biomass of a target species that does not cause unsustainable impacts on the species or ecosystem (McKindsey et al., 2006). Thus, to assess the ECC of bivalve shellfish, their impact on the ecosystem structure, function, and other functional groups needs to be examined at a holistic level. Food web models provide the means to conduct such research at the ecosystem level, with the Ecopath with Ecosim (EwE) model currently being one of the most popular models (Feng et al., 2022a). This model has been widely applied, for example, to calculate an ECC of 65 t/km 2 for Perna canaliculus in the Tasman and Golden Bays (Jiang and Gibbs, 2005), 297 t/km 2 for oysters farmed in Narragansett Bay (Byron et al., 2011), and 976 t/ km 2 for Crassostrea gigas in Sanggou Bay (Gao et al., 2020). Feng et al. (2022a) also used the EwE model to calculate the ECC of Mercenaria mercenaria soft tissue in monoculture and mixed culture ponds as 41.42 g/m 2 and 470.00 g/m 2 , respectively. It assesses the ECC of shellfish based on the principle of trophic dynamics and from the perspective of material and capacity balance, which fully considers the competitive relationship between species and trophic groups (Srithong et al., 2021).
The Bohai Bay, located on the east coast of China, has fertile water quality and abundant germplasm resources. Numerous oyster reefs and shell ridges once existed along its northwestern reaches (Yue et al., 2012), but currently, only patches of oyster reefs remain, distributed in the estuaries and shallows. This decline resulted from the rapid increase in destructive fishing activities and human interference in the coastal areas of the Bohai Bay, coinciding with advancing urbanization (Jin and Deng, 2000). To alleviate this ecosystem decline, the construction of artificial oyster reefs via marine ranching have been implemented at the confluence of the Liuan River and Bohai Bay since 2011. The marine ranching study located in Bohai Bay was recognized as one of the first National Marine Ranching Demonstration Areas in China in 2015. We hypothesized that a high bivalve shellfish biomass may have an adverse effect on the stability of system maturity and the abundance of other species in the region. As a result, the EwE model was implemented in this study to evaluate the ECC of bivalve shellfish in marine ranching areas and to estimate their C sequestration potential. The research results may offer a scientific reference for future assessments of the ecological service function of marine ranching and the sustainable development of the coastal region.

Study area
The study area was located in Bohai Bay, China (39°11'49.9" N, 118°59'25.3" E) (Figure 1), in which 93,0000 m 3 of various artificial reefs had been deployed by the end of 2018, including two artificial oyster reef areas of 200 and 57 ha, respectively. This study focused on the 57 ha oyster reef, for which construction was initiated in December 2017 and the reef was completed in June 2018. The artificial reef comprises granite stones or caisson-components consisting of reinforced concrete structures. It provided the habitat for a large number of Mytilus edulis to attach. The water depth in this area is 8.22 m, the average annual water temperature is 13.22°C, the average current velocity is 0.28-0.37 m/s, and there are no macroscopic features.

Ecological model construction 2.2.1 Ecopath model construction
The Ecopath model is commonly used to quantitatively estimate mass transfer and energy flow in an ecosystem (Christensen and Walters, 2004;Christensen et al., 2005). According to thermodynamic theories, the energy input and output of a certain biological functional group should be balanced, with productivity equaling the sum of mortalities (Equation 1). Another equation was used to determine the inner mass balance of the group (Equation 2), where the consumption (Q) of each group was equal to the sum of the production (P), respiration (R i ), and unassimilated food (U i ) of the group. The equations are as follows: where B i is the biomass of functional group i, (P/B) i is the production ratio, (Q/B) i is the consumption ratio, EE i is the ecotrophic efficiency (the proportion of production utilized in the system), DC ij is the proportion of prey i in the diet of predator j, Y i is the fishery mortality, BA i is the biomass accumulation rate, and E i is the net migration rate. In the parameter input process, B i , (P/B) i , (Q/B) i , and DC ij are required. If one of the first three values is missing, the corresponding EE value can be designated. In addition, E i and catches are also required.

Ecosim dynamic simulation
Ecosim is a time-scale based dynamic model based on Ecopath (Christensen and Walters, 2004). It incorporates a time-dynamic aspect by varying the Ecopath model over time increments, simulating systemic response changes over the time series. The equation used is as follows: where dB i /d t is the rate of change in biomass, g i is the net growth efficiency, Q ji is the consumption rate of function group j to function group i, I i is the emigration rate, F i is the fishing mortality rate, Mo i is the non-predatory natural mortality, e i is the move-in rate, and B i is the biomass of functional group i (Christensen and Walters, 2004).

Input data 2.3.1 Division of functional groups
We constructed 21 functional groups based on the ecological habits, economic values, and ecological roles of the marine ranching species in the study area. Among them, the bivalve shellfish Geographical location of the artificial oyster reef area in the study marine ranching in Bohai Bay, China. Li et al. 10.3389/fmars.2023.1174235 Frontiers in Marine Science frontiersin.org comprised C. gigas and M. edulis, accounting for 86.00% and 14.00% of the bivalve shellfish biomass, respectively. The constructed functional groups covered the main energy flow in the 57-ha oyster reef area (Table 1).

Source of parameters
The units for system energy flow in our EwE model were all wet-weight units (t/km 2 ). Biomass measurements were obtained via four biological surveys (summer: June 2019, autumn: September 2019, winter: December 2019, spring: March 2020), in which fish and cephalopod biomass were estimated mainly from trap and gill nets, while trawling and underwater video footage were used for corrections. Fishing was conducted for one 72 h period during each of the four surveys; the nets used had a mesh size of 2 cm with the section, length, width, and height dimensions of a net being 5 m, 8 m, 30 cm, and 20 cm, respectively. Fish biomass was derived using the method described by Li et al. (2021). Macrobenthic biomass was estimated by combining the videos and SCUBA grasping. Bivalve shellfish were collected from the oyster reefs by SCUBA diver using 0.5×0.5 m quadrats, with the locations of 3-6 replicates being evenly distributed in each season. Meiobenthos fauna was sampled at each site using box samplers (0.05 m 2 ). Zooplankton samples were obtained via vertical towing of a plankton net with a mesh size of 169 mm. Phytoplankton biomass was calculated from chlorophyll readings as measured using standard procedures described by Parsons et al., 1984. Detritus was obtained as the organic matter in the water column by filtering 1 L of subsurface and bottom seawater through a mesh size of 77 μm. Physical factors were measured in situ using a YSI meter (EXO2, Yellow Springs, OH, United States).
The P/B (production ratio) of fish is equal to the total instantaneous fish mortality value, which is equal to the sum of the natural mortality (M) and fishing mortality (Palomares and Pauly, 1989). The calculation formula is as follows: where Y is the catch, B is the biomass, K is the curvature of von Bertalanffy growth equation, L ∞ is the asymptotic body length (total length, cm), and T c is the average water temperature (°C). The P/B of invertebrates is calculated as follows (Tumbiolo and Downing, 1994).
where B is the biomass and W is the mean individual mass. The formula used to calculate the consumption value (Q/B) is as follows: where W inf is the asymptotic weight (g), T is the average annual water temperature in the water body (T = 1000/K), A is the ratio of caudal fin length to width of fish (A = h 2 c =s, where h c is the height and s is the surface area of the caudal fin), h is the set index of variable food types (1 for omnivorous, 0 for carnivorous and herbivorous), and d is the dummy variable (d = 0 for herbivorous and carnivorous fish, and d = 1 for detritivorous fish).
The P/B and Q/B values of other functional groups were derived from relevant studies conducted at the same latitude and in adjacent seas (Lv et al., 2018). Food composition matrices were obtained by gastric content analysis and by reference to studies in nearby seas (Zhang, 2018) and the FishBase database (www.fishbase.de). (The edibility matrix is provided in Table S1, and the specific parameter sources are listed in Table S2). The vulnerability value in Ecosim determines the trophic control relationship between the predator and prey (i.e., upstream, downstream, and intermediate effects). In this study, we referred to the calculations by Kluger et al. (2016). The formula used to calculate the vulnerability value is as follows: where TL i is the trophic level of functional group i. The vulnerability takes values range from 0 to 1, where 0.1-0.3 indicates a downward control, 0.3 a mixed control, and 0.3-1 an upward control effect (Christensen and Walters, 2004). The range of vulnerability in Ecosim is 1 to ∞; thus, v i must be converted to v new :

ECC estimation and dynamic simulation
To reflect a balance and stability of the marine ecosystem in the Ecopath model, the biomass of economically and ecologically important organisms in the ecosystem is gradually increased until the ecotrophic efficiency (EE) value of any group is 1; at this point, the ECC of that species group is obtained (Byron et al., 2011).
Three scenarios were investigated through dynamic simulation: maintaining the status quo (current) in terms of population size, reaching the ECC, and exceeding the ECC by 20%. The target species stocking expansion was simulated for a period of 30 years under scenarios of differing final enhancement biomass, which was implemented through a linear increase in stock from 2 to 6 years and then held constant for the remaining 24 years. The time series of changes in the relative biomass was extracted for each simulated scenario. The C budget parameters of both C. gigas and M. edulis were obtained from Ren (2014).

Model equilibrium and uncertainty analysis
The Ecopath model equilibrium was evaluated according to the parameters employed by Heymans et al. (2016). A model was considered balanced when it satisfied the following three requirements: (1) The EE values were less than 1; (2) the nutrient conversion efficiency P/Q was between 0.1 and 0.3, with potentially lower and higher values for top predators and fast-growing species, respectively; and (3) the net nutrient conversion efficiency of all functional groups were less than 1.
To test whether the parameters modified by balancing were in accordance with the basic rules of ecosystem ecology, P/B and Q/B pre-balancing between taxa and trophic levels was performed using the PREBAL module (Link, 2010).The prerequisites were: (1) functional group biomass should span 5-7 orders of magnitude; (2) the slope of functional group biomass on a logarithmic scale should decrease by 5-10% across all taxa relative to the trophic arrangement; and (3) P/B and Q/B across taxa and trophic levels should decrease with increasing trophic levels. The results of the pre-balance analysis are shown in Supplementary Figure 1.
For the B, P/B, Q/B, and DC values, the spectrum index can evaluate the overall quality of the ecosystem model by assigning uncertainty indices according to the data source and combining the uncertainty of each functional group. The sensitivity of the Ecosim results to the Ecopath input parameters (B, P/B, and Q/B) was evaluated by performing 500 Monte Carlo simulations of these input parameters, where each parameter was set to vary by 10%.

Ecological network analysis
An ecological network analysis plug-in was built into the EwE software for direct analysis (Christensen and Walters, 2004), and a set of ecological indices described the structure and function of ecosystems. The following indicators were used to characterize the overall ecosystem: total system throughput (TST), total primary production (TPP), total biomass (TB), and total respiration (TR), which describe the overall activity and size of the ecosystem (Latham, 2006;Ortiz et al., 2015). We assessed the ratios of TPP/ TR and TPP/TB as important indicators of ecosystem stability and maturity (Odum, 1969). Finn' s cycling index (FCI), connectance index (CI), systematic omnivory index (SOI), and finn's mean path length (FMPL) expressed the complexity of food webs (Finn, 1976;Christensen and Walters, 2004;Libralato, 2013). Ascendency (A) expressed the development and maturity of the ecosystem, and overhead (O) indicated the ability of the ecosystem to resist disturbances (Ulanowicz, 1986). Entropy (H) indicated the amount of energy interaction between all groups of organisms in the system, and the average degree of system energy flow (Baird et al., 2007). Average interaction information (AMI) represented the determinism of the degree of connectedness (i.e., energy flow), among the biological groups of the system (Ulanowicz, 2004). Trophic transfer efficiency (TE) is the ratio between, on the one side, the sum of the exports and the flow transferred to the next TL, and, on the other side, the throughput on the TL; this study used a mean TE for TL > 2. The Shannon diversity index (SDI) was assessed as a relative measure of biodiversity. Finally, radar plots were drawn to show the magnitude of the changes in the ecological network finger table and to describe the changes in energy transfer before and after oyster value addition using Lindeman energy flow diagrams.

Ecological model quality analysis
The PREBAL diagnostic results were consistent with the general trend displayed by the model (Heymans et al., 2016), with the high biomass of bivalve shellfish being the result of the construction of an artificial oyster reef. Phytoplankton and zooplankton displayed a high P/B ratio because they were the main food sources for the bivalve shellfish (Supplementary Figure 1).
The confidence index of the Ecopath model was 0.603, indicating that the input data was reasonable and provided a good confidence level. The results of the Monte Carlo simulations showed that the optimal input parameters for B, P/B, and Q/B were consistent with our actual input parameters and did not change, except for the bioaccumulation value for Stichopus japonicus, which should be raised to 0.20.

Evaluation of the ECC of bivalve shellfish
The output parameters of the Ecopath model are presented in Table 2. The highest biomass in the marine ranching ecosystem was that of bivalve shellfish (187.76 t/km 2 ), followed by detritus (40.00 t/ km 2 ), while Hexagrammos otakii fish had the lowest biomass (0.20 t/km 2 ). In this ecosystem, gobiid fish and annelids shared the highest EE value of 0.97; Starfish had the lowest EE value (0.00), followed by skeleton shrimp (0.11). In our modeling of bivalve shellfish ECC, the EE of phytoplankton was equal to1 when the bivalve shellfish biomass reached 282.66 t/km 2 . When the bivalve shellfish biomass was increased to 339.19 t/km 2 , the EE value of phytoplankton was greater than 1, and the Ecopath model was unbalanced (Table 3). As a result, it was concluded that the ECC of bivalve shellfish in the study marine ranching is 282.66 t/km 2 .
3.3 Analysis of the overall characteristics of the study marine ranching ecosystem The TLs of each functional group in the study marine ranching ecosystem are presented in Table 4. Cephalopods had the highest TL (4.15), followed by the fish Platycephalus indicus (3.77) and Sebastes schlegelii (3.58). The energy transfer in the marine ranching of the study marine ranching was mainly concentrated in TL I and TL II, with the highest energy throughput of detritus in TL I and the highest energy throughput of bivalve shellfish in TL II.
Radar plots of the ecological network analysis showed that system parameters responded to different degrees as the biomass of bivalve shellfish was simulated to increase ( Figure 2). The TST, TB, TR, and TPP (which characterize the system size) showed an increasing trend, with the most obvious increasing parameter being TB. Among the parameters characterizing system maturity and stability, TPP/TR, TPP/TB, H and O/C decreased, whereas FMPL, FCI, A/C, and AMI increased. SOI and CI did not change. We also found that TE, which describes energy transfer, and SDI, which characterizes system biodiversity, showed decreasing trends (see Table S3 for specific parameters). The Lindeman energy flow diagram (Figure 3) revealed that the total energy transfer efficiency of the system decreased as the bivalve shellfish biomass increased, but the transfer efficiency of primary producers and detritus to TL II increased, where primary producers were elevated more than detritus. There was a slight decrease in TL II III and TL III IV.

Effect of bivalve shellfish proliferation on other functional groups
The Ecosim dynamic simulations showed that the other functional groups were affected to varying degrees by increasing bivalve shellfish biomass in our three scenarios ( Figure 4). As bivalve shellfish proliferated to their ECC, the biomass of H. otakii, pelagic fishes, gobiid fish, skeleton shrimp, annelids, other macrobenthos, meiobenthos, S. japonicus, zooplankton, phytoplankton, and detritus decreased. The most severe negative effect was observed for gobiid fish, which showed a 99.38% reduction in biomass compared in the 10th year of the simulation that continued to decrease. The biomass of phytoplankton, the main food source of bivalve shellfish, did not fall below 85.00% of the initial biomass by the end of the simulation. The biomass with the highest increase was that of snails, with a maximum fold change of 7.18.
When the bivalve shellfish biomass was simulated to exceed their ECC by more than 20%, the biomass of skeleton shrimp, annelids, other macrobenthos, meiobenthos, S. japonicus, and detritus showed the opposite response to the scenario under bivalve shellfish proliferation to ECC; the biomass of bivalve shellfish decreased when they proliferated to ECC but increased when they proliferated to 120%ECC. By the end of the simulation, the biomass of other macrobenthos, meiobenthos, and S. japonicus exceeded that of scenario 1. Biomass changes in the remaining functional groups were consistent with those observed under bivalve shellfish proliferation to ECC.

C budget levels of individuals and populations of bivalve shellfish
During one growth cycle (360 d), one C. gigas individual ingested 14.97 g of C, of which 22.19% was stored on the seafloor through biodeposition, 26.92% was removed from the water through harvesting (with 89.33% of this amount being used for shell formation), and 50.89% was released into the water through respiration, calcification, reproduction, and excretion. Over the same period, one M. edulis ingested 9.19 g of C, of which 38.10% was stored on the seafloor through biodeposition, 11.49% was harvested from the water (with shell formation accounting for 79.25% of the C usage), and 50.41% was released into the water through respiration, calcification, reproduction, and excretion ( Figure 5).
The average wet weight of C. gigas and M. edulis was 42.71 ± 6.08 g and 19.11 ± 3.20 g, respectively. Based on the total bivalve shellfish biomass in the study area (187.76 t/km 2 ), these values correspond with population sizes of 2.16 million and 0.79 million for C. gigas and M. edulis, respectively, in the study marine ranching. In total, these bivalve shellfish ingested 39.52 t of C in 360 d, of which 9.92 t was stored on the seafloor through biodeposition, 9.52 t was removed from the water through harvesting (including 8.42 t being absorbed for shell formation), and 20.08 t was used for respiration, calcification, reproduction, and excretion, and then released back into the water body. When the biomass of the bivalve shellfish reached ECC (282.66 t/km 2 ) in our model, their number increased to3.24 million, and that of M. edulis increased to 1.18 million. This translated to the ingestion of 59.41 t of C during the 360 d growth cycle, of which 14.91 t was biodeposited on the seafloor, 14.32 t was harvested from the water (including 12.67 t being absorbed by shell formation, and 30.18 t was used for respiration, calcification, reproduction, and excretion and was released back into the water body ( Figure 6).

Overall maturity and stability of the study marine ranching ecosystem
In this study, we analyzed ecosystem scale, system maturity, stability, and trophic transfer to predict the potential impact of an increased bivalve shellfish biomass on the study marine ranching ecosystem.
An increase in bivalve shellfish biomass resulted in an increase in the TST, TB, TR and TPP of the ecosystem, indicating that the overall scale of the system gradually expanded. According to Odum (1969) and Christensen et al. (2005), both TPP/TR and TPP/TB are important parameters that respond to ecosystem maturity: when an ecosystem develops maturely, both TBB/TR and TPP/TB approach 1. In this study, TPP/TR and TPP/TB showed a decreasing trend with an increase in bivalve shellfish biomass and were 1.40 and 9.26, respectively, when the bivalve shellfish reached ECC. This indicated that the system matured with the increase in bivalve shellfish biomass. As an ecosystem matures, the food chain changes from linear to reticulate, with increased values of TST, FCI, FMPL, CI and SOI (Libralato, 2013). In our study, FCI and FMPL both increased significantly as bivalve shellfish biomass increased, comparable to the rates described by Heymans et al. (2011) for the 1-10 km 2 ecosystem. This increase in FCI and FMPL can also indicate that the energy utilization of the detritus functional group improved and that the path of material circulation is longer. This marine ranching ecosystem was compared with the models proposed by Wu et al. (2016) and Heymans et al. (2011). The SOI and CI values before and after bivalve shellfish proliferation in this area were relatively low, indicating that the food web of the system was simple and linear. At the same time, A/C and O/C are parameters that characterize the stability and maturity of the ecosystem. The relatively high A/C ratio indicates that the ecosystem is well organized and unlikely to disintegrate spontaneously. However, a decrease in the O/C value may make the system vulnerable to external pressure (Feng et al., 2022b), and the system is the most stable when the ratio of A to C is 45.96 (Ulanowicz, 2012). In our study area, the prevailing A/C value was 23.61%, indicating that the ecosystem was not sufficiently stable. However, as we modeled at bivalve shellfish biomass increased to ECC, the A/C increased to 28.88%, indicating that it had stronger self-regulation ability, but that it was still not sufficiently stable. The A, C to TST ratios, AMI and H are also important parameters for describing the stability of system maturity. According to Odum (1969) and Finn (1976), a mature ecosystem should have a high FCI and low H value. In this study, as the bivalve shellfish biomass increased to the ECC, FCI increased by 1.32 times and H decreased by 7.30%, implying that the ecosystem displayed a higher maturity and stability at this time point. The corresponding increase in AMI (from 0.91 to 1.04bits) indicates that the ecosystem is more restricted at this time and energy flows along a more fixed path.
Through ecological network analysis, we found that the SDI decreased with an increase in bivalve shellfish biomass. As SDI characterizes the diversity of a system, this result indicated that a higher bivalve shellfish biomass would have a negative impact on ecosystem biodiversity. We also found a slight decrease in the total system TE when ECC was reached. Lindeman energy flow diagrams indicated, that the TE from detritus and primary production to TL II increased as the bivalve shellfish biomass reached the ECC, with a more pronounced increase in the transfer efficiency of primary production. The main reason for this is that bivalve shellfish account for more than 70.00% of TL II, and they mainly feed on phytoplankton, which leads to an increase in phytoplankton utilization. Simultaneously, bivalve shellfish effectively store energy in macrobenthos through predation (Wang et al., 2022). Therefore, an increase in bivalve shellfish can better sustain the biomass of economically important macroinvertebrates, which can generate considerable economic benefits. We could conclude that the study marine ranching ecosystem is not yet mature, but the system can be expected to mature and stabilize as the bivalve shellfish biomass increases towards the ECC.

Effect of bivalve shellfish proliferation on other functional groups
According to the results of the Ecosim simulation, an increase in bivalve shellfish biomass can have an indirect or direct effect on other functional groups. We first examined changes in other functional groups as the bivalve shellfish biomass was simulated to increase from its current level to its ECC. The main food of bivalve shellfish is phytoplankton; therefore, phytoplankton availability is the biggest limiting factor and hypothetical contradiction. We found that the biomass of phytoplankton never decreased to less than 85.00% of the initial biomass at any time during simulations. This may be due to the low vulnerability and low EE value of bivalve shellfish, which ultimately limits the Radar plots show the relative changes in the analyzed indices of ecological forgetfulness for the two scenarios (current, ECC) predicted according to the EwE model. ECC, ecological carrying capacity; EwE, ecopath with ecosim.
increase in the predation pressure on phytoplankton. As a result, bivalve shellfish have the opportunity to proliferate without running out of food resources. However, the decrease in phytoplankton biomass affected other functional groups that compete with bivalve shellfish for phytoplankton as a food source. Gobiid fish was the most drastically affected, displaying a 99.99% reduction in their biomass by the end of the simulation. One explanation is the bivalve shellfish proliferation, because gobiid fish and bivalve shellfish have more than 50% of their prey in common, such as skeleton shrimp, annelids, and zooplankton; another explanation is that the construction of oyster reefs may also reduce the gobiid living space. Zooplankton species are some of the most important food competitors of bivalve shellfish, and the biomass of bivalve shellfish is understandably reduced by an increase in zooplankton biomass. However, the biomass of zooplankton in our study never dropped below 85.00% of their original biomass. This may be because that pelagic fish are the main predators of zooplankton in the study area and are also competitors for bivalve shellfish prey. The increase in bivalve shellfish biomass led to a decrease in the biomass of the pelagic fish population, which indirectly reduced their predation pressure on zooplankton, resulting in a relatively small drop in zooplankton biomass. Furthermore, the increase in bivalve shellfish biomass also led to an increase in the biomass of Charybdis japonica, snails, and starfish.
When the bivalve shellfish biomass was expanded to 120%ECC, the biomass of some functional groups showed the opposite trend compared to their response to a 100%ECC bivalve shellfish proliferation. This may be since that the further proliferation of bivalve shellfish caused a large reduction of pelagic and gobiid fish, which mainly prey on phytoplankton and zooplankton and therefore compete with bivalve shellfish for food resources. This could result in a corresponding increase in the biomass of other organisms that feed on plankton, which in turn allows an in the biomass of some higher trophic levels. Starfish no natural enemies in the study area and bivalve shellfish are their main food source. Therefore, an increasing bivalve shellfish biomass causes the starfish population to expand continuously. Starfish have low economic value in this region, but to avoid them overfeeding on bivalve shellfish, they must be adequately harvested. We believe that the appropriate harvesting of higher trophic-level species is required in the later stages of marine ranching management to ensure ecosystem stability.
Bivalve shellfish are considered the most sustainable aquaculture species, mainly because they prey on naturally occurring phytoplankton and do not require external bait inputs (Dumbauld et al., 2009). The proliferation of bivalve shellfish generally requires the addition of settlement structures for them, thus changing environmental conditions while potentially Lindeman energy flow diagram between trophic levels for two scenarios for (A) current and (B) ECC. ECC, ecological carrying capacity; D, detritus; P, primary producers; TL, trophic level; TE, transfer efficiency; TST, total system throughput.
providing habitat for other organisms and contributing to an increased species diversity (Filgueira et al., 2021). However, the expansion of bivalve shellfish biomass can have certain negative effects on ecological structure and species diversity. To alleviate the drastic decline in the biomass of other species during bivalve shellfish stocking, while also increasing primary production in the water column, the appropriate proliferation of macroalgae as a food source may be a solution. Our results further showed that an increase in bivalve shellfish resulted in a moderate increase in the biomass of the main economic species in the area, which could increase the economic benefits of the area, in agreement with the findings of Loṕez-Jamar et al. (1984). Ecosim simulations of future scenarios indicate that the biomass of bivalve shellfish can be expanded to a maximum of 282.66 t/km 2 , above which the impact on other species is relatively large. Our study also shows that it is necessary to model the ecosystem of marine ranching The relative changes in biomass of each functional group under three scenarios (current, ECC, and 120% ECC) of bivalve biomass augmentation in the study marine ranching ecosystem model were predicted. The biomass of bivalves was kept constant after the dashed line. ECC, ecological carrying capacity; 120%ECC, 120% ecological carrying capacity. before multiplication and release activities are implemented.
4.3 C budget at the population level of bivalve shellfish in the study marine ranching Tang and Hui (2016) showed that shellfish function by filtering water, calcification, and biological deposition, and their physiological characteristics make them an important biological C sink. Chauvaud et al. (2003) proposed that shellfish not only achieve C sequestration through biological deposition and harvesting, but also emit C through respiration, calcification, and reproductive excretion; therefore, shellfish are not pure C sequestrators. Generally, if more C is used than released by an organism, it is a C sink; if not, it is a C source. In the process of shellfish calcification, for every 1 mol of CaCO 3 formed, 1 mol of CO 2 is released, but 2 mol of HCO − 3 can be absorbed, which indicates that the removal of C is greater than its release (Zhang et al., 2005).  showed that the CO 2 produced by bivalve shellfish is not directly released into the atmosphere, phytoplankton and macroalgae absorb some of the CO 2 , the bivalve shellfish themselves reuse a proportion to form shells, and the rest is released into the atmosphere. Overall, bivalve shellfish and macroalgae such as phytoplankton can synergize to enhance the C sink capacity of coastal zone ecosystems (Feng et al., 2022c). According to our model, if the biomass of bivalve shellfish is increased to the ECC in the marine ranching, the C intake of the bivalve shellfish will increase by 19.89 t, C biodeposits on the ocean floor will increase by 4.99 t, harvested C will increase by 4.80 t, and the C redischarged into the water by respiration and excretion will increase by 10.10 t. Furthermore, bivalve shellfish proliferation is generally achieved by placing artificial reefs, and the physical deposition of C on reefs populated with bivalve shellfish is an important process (Newell et al., 2005). Reefs generally have a complex physical structure, which helps to slow down the water current reduce wave energy and suspended matter deposition, thus promoting sediment accumulation and organic C burial (Chowdhury et al., 2019).
This study has several limitations. First, due to the short completion time of the Ecosim model, the lack of previous data did not fit the time series, which may induce some uncertainties in our results. Second, we used stomach content analyses to construct the diet matrix in this region; however, such an analysis can easily affect the reliability of the data by only reflecting an instantaneous prey composition (Gee, 1989) of an organism, not its long-term feeding characteristics. In future research, stable C and nitrogen isotopes can completely retain ingestion information to clarify food sources (Fry, 2006).

Conclusion
In recent years, under the multiple pressures of global climate change and human activities, the Chinese marine environment has changed, and the degradation of offshore biological resources has become increasingly serious. Resultantly, the construction of artificial reefs and stocking activities via marine ranching has been accelerated with the of repairing damaged habitats and conserving biological resources. This study used ecotrophic modeling and dynamic simulations to estimate an ECC of 282.66 t/km 2 for bivalve shellfish in the study marine ranching in Bohai Bay, which would be beneficial to the general development of the ecosystem but may inflict a negative impact on some species. In our simulation, negative responses to bivalve shellfish proliferation were A B

FIGURE 5
Carbon budget in one growth cycle of (A) Crassostrea gigas and (B) Mytilus edulis (Modified from Ren, 2014). DIC, dissolved inorganic carbon. Li et al. 10.3389/fmars.2023.1174235 Frontiers in Marine Science frontiersin.org mainly due to the lack of primary productivity, but the addition of macroalgae can alleviate these effects while assisting with C sequestration (Watanabe et al., 2020). Via the synergistic ecological functions of bivalve shellfish and macroalgae, the C sequestration capacity of ecosystems can be effectively improved.

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 authors.