The gross ecosystem product of oyster reefs

Oyster reefs are critical interfacial habitats that underpin coastal and marine bio-diversity, and enhance multiple ecosystem services of blue carbon habitats, including water filtration and shoreline stabilization. However, despite this ecological value, oyster reefs have suffered dramatic global declines. The absence of a systemic economic valuation of oyster reef ecosystem services introduces deep uncertainty in the prioritization of restoration areas, thereby constraining coordinated large-scale transboundary conservation, restoration, eco-engineering efforts, and policy. Here, we novelly introduce and estimate the gross ecopotential, or Gross Ecosystem Product (GEP), as a systemic index of reef ecosystem services grounded in ecological fitness. Fitness is assessed by habitat suitability (HS) at a ∼9.2 km resolution globally. The ecopotential is calculated as a multicriteria socio-ecological value function combining potential oyster trade value with monetary scenarios for the quantified shell-carbon component, water filtration and nutrient removal. Potential effects on adjacent vegetated blue-carbon habitats are discussed as conceptual pathways but are not included in the numerical carbon valuation. We highlight ten hotspot countries – the USA, Australia, Chile, China, Japan, South Korea, Argentina, France, Germany, and Italy – where HS and restoration opportunities are the highest. Seabed light availability and minimum primary productivity, reflecting phytoplankton biomass, emerged as the dominant global drivers of oyster reef HS and GEP. At the reference threshold (HS ≥ 0.65), approximately 23% of predicted high-HS habitat overlaps the evaluated marine and coastal protected-area network, and approximately 10% of high-HS cells in the hotspot analysis coincide with occurrence records. Because the presence-only archive is incomplete and unevenly sampled, high-HS cells without an occurrence record are interpreted as candidates for further survey rather than confirmed oyster absences. Hydrogeomorphic analyses indicate that ecosystem morphology drives GEP through distinct routing of ecological flows: ∼72% of reefs occur in estuaries, 21% in lagoons or inlets, and 7% in river deltas. Our findings highlight the considerable GEP of oyster reefs when conservation, restoration, and regeneration are interconnected actions with systemic and global impact. The proposed GEP approach integrates publicly available socio-ecological and economic data to generate decision-relevant priority insights for restoration planning. Digital Eco-Information layers, identifying existing and potential reef habitats with high GEP, are provided, including their hydrogeomorphic flows and cascading ecosystem services. The model is readily tunable to higher-quality data and broadly extensible to other biogenic species and habitats, providing a scalable baseline for blue-carbon ecosystem management and future scenario development.

Systemic pressure and nature-positive design

Coastal and marine ecosystems, and in particular ecotones, are experiencing rapid global decline due to local and global pressure, such as natural resource overexploitation, habitat modifications, systemic pollution due to land-to-ocean efflux, and global climate change. This decline threatens biodiversity and ecosystem stability, yet the benefits that these ecosystems provide to society include fisheries, coastal protection, and water quality. The United Nations Sustainable Development Goal 14 (SDG 14, “Life Below Water”) explicitly calls for conserving and sustainably using coastal and marine resources to safeguard human wellbeing in a broad ecosystem health framework (Duarte et al., 2020). However, a comprehensive ecosystem management approach that embeds spatial habitat socio-ecological networks and flows is still lacking, along with a robust economic valuation of ecosystem services [i.e., Gross Ecosystem Product (Ouyang et al., 2020)] needed to prioritize critical areas for enhanced conservation, restoration, and habitat formation within a nature-positive design value framework.

Nature-positive design also embraces the self-emergence of nature, such as the one of habitat-forming species (forming biogenic structures like oyster reefs) as critical ecosystem engineers with high collective dynamics. Within coastal systems, oyster reefs are foundational habitats that deliver multiple ecological and socio-economic services. Dense collective aggregations of oysters create complex three-dimensional structures that enhance biodiversity (with cascading food-web regulation), water filtration, and sediment regulation resulting in shoreline stabilization (Coen et al., 2007; Fodrie et al., 2017). By improving water clarity and modifying sediment and hydrodynamic conditions, oyster reefs may also support adjacent seagrass, mangrove, and salt-marsh habitats; these interactions are context dependent and are treated conceptually here rather than as quantified carbon gains (Fodrie et al., 2017; Macreadie et al., 2021; Wu and Convertino, 2026). Despite this socio-ecological importance, the global extent of oyster reefs has declined by more than 85% in the past century (Beck et al., 2011), with particularly severe losses in estuaries and deltas that are ecosystems subjected to intensive human use, transformation and climate change extremes (Lau et al., 2024).

Oyster reefs as blue carbon ecotones: a flow symbiosis

Despite increasing recognition of oyster reef importance, three major gaps hinder large-scale oyster reef restoration and conservation.

First, global assessments identifying where suitable habitats for oyster reefs currently exist, and where they could sustain or support self-emergent reef formation, are lacking. Most existing studies focus on individual species or local estuaries (Cake, 1983; Chowdhury et al., 2019; Linhoss and Mickle, 2022), overlooking broad-scale biogeo-graphic patterns and eco-environmental drivers.

Second, the ecological and economic values of oyster reefs are rarely quantified beyond local contexts, and previous valuations typically considered only harvest value rather than the broader portfolio of ecosystem services (Beseres Pollack et al., 2013; Grabowski et al., 2012). These gaps limit the ability of policymakers to prioritize investments and to integrate oyster reefs into climate adaptation and blue-carbon strategies.

Third, oyster reefs are rarely considered explicitly within broader blue-carbon systems (Zhang and Convertino, 2024, 2025; Wu and Convertino, 2026; Zhao et al., 2026) (Figure 1). Oyster reefs operate as dynamic coupling nodes, or distributed ecotones, within coastal systems that include mangroves, seagrasses, and salt marshes. Mangroves and marshes export dissolved and particulate organic carbon, while seagrasses fix carbon in situ through photosynthesis. Oyster reefs intercept material flows, trans-form particles into biomass, carbonate shell, and biodeposits, attenuate hydrodynamic energy, promote sediment retention, and improve water clarity. These processes may support adjacent vegetated habitats, but oyster reefs do not function as uniform net carbon sinks (Fodrie et al., 2017). In the present study, adjacent-habitat pathways provide conceptual context and are not included in the numerical carbon valuation.

Oyster reef ecopotential assessment

To quantify the suitability and value of oyster reefs worldwide, we introduce the concept of ecopotential, defined as the integration of ecological habitat suitability (HS) with the socio-economic value of ecosystem services associated with HS. Using global presence-only data and environmental covariates, we apply MaxEnt (Phillips et al., 2006; Elith et al., 2011) to map oyster reef suitability and then quantify associated benefits, including oyster trade potential, water filtration, nutrient removal, and a scenario-based valuation of the quantified shell-carbon component. MaxEnt has been applied to regional oyster habitat-suitability modelling (Linhoss and Mickle, 2022), whereas field-validated habitat-suitability-index frameworks have linked oyster performance to region-specific environmental conditions (Theuerkauf and Lipcius, 2016; Chowdhury et al., 2019). These applications remain species-specific or regional. We treat HS as a first-order screening indicator of oyster-reef potential, recognizing that actual reef development depends on fine-scale biological, geomorphic, and socio-economic conditions. Our goal is therefore to estimate global ecopotential for assessing broad conservation and restoration priorities. We identify ten hotspot countries with the greatest modelled restoration potential and evaluate overlap with Marine Protected Areas (MPAs). While not a mechanistic model of reef dynamics, this approach synthesizes available ecological and economic information to provide decision-relevant insights based on HS-linked ecopotential patterns. By linking ecological suitability with socio-economic valuation, we highlight oyster reefs as strategic assets for restoration, conservation prioritization, and future blue-carbon scenario development.

An occurrence record denotes a georeferenced observation reported by GBIF or iNaturalist and does not delineate a reef footprint. HS denotes the continuous relative-suitability output from MaxEnt; high-HS cells are cells in the highest post-processing suitability class (HS ≥ 0.65) and are candidate areas for assessment rather than confirmed reefs or guaranteed restoration sites. Ecopotential denotes the HS-linked ecosystem-service scenario framework, whereas GEP denotes its monetary valuation component.

The analysis comprised six stages. First, oyster occurrence records were acquired from GBIF and iNaturalist. Second, coordinates were standardized and records lacking coordinates were removed. Third, nine benthic Bio-ORACLE predictors were prepared on a common coastal grid. Fourth, candidate MaxEnt configurations were calibrated and evaluated internally. Fifth, the continuous HS surface was classified and intersected with country and protected-area boundaries for spatial prioritization. Sixth, high-HS areas were linked to standardized ecosystem-service and economic-valuation scenarios.

Global occurrence of oyster reefs and biogeochemical determinants

Because no global inventory of oyster-reef extent is available, we used public presence-only oyster observations as a proxy for broad habitat occurrence. Records were retrieved from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/) and iNaturalist (iNat, https://www.inaturalist.org/) with the R package spocc. The active call was spocc::occ(query = c(“Ostreidae”),from = c(“gbif”, “inat”),limit = 6000), following the public rspatialdata Species Occurrence workflow; no observation-date interval was active. The returned dataset contains 10,442 records, comprising 5,977 GBIF and 4,465 iNaturalist records. The retained snapshot includes the data provider, coordinates, observation date, provider-specific record identifier, and returned taxon-name field. Ostreidae is a family-level query, not a single-species name. These records are reported oyster observations and do not provide verified reef boundaries, reef condition, or natural-versus-cultivated status (Figure 2).

The verifiable processing sequence converted longitude and latitude to numeric values, removed records with missing coordinates, converted the remaining observations to spatial point features, and assigned country attributes by a left spatial join. The occurrence map is a visualization of record distribution, not a statistical correction for sampling bias. The workflow did not apply a complete iNaturalist research-grade restriction, taxonomic-backbone reconciliation, automated geographic-outlier screening, coordinate-based deduplication, or spatial thinning. These procedures can materially affect presence-only models (Beck et al., 2014; Boria et al., 2014; Aiello-Lammens et al., 2015; Zizka et al., 2019); consequently, the pooled family-level output is interpreted as a global screening layer rather than a species-specific restoration prescription. A species- or region-resolved analysis would require harmonized taxonomic validation, bias-aware sampling, and independent field data.

We assembled nine present-day benthic predictors from Bio-ORACLE v2 (Tyberghein et al., 2012; Assis et al., 2018): minimum primary productivity (g m −3 day −1), mean temperature (° C), temperature range (° C), mean light at the seafloor (E m −2 yr −1), maximum phytoplankton carbon concentration (µ mol m −3), minimum dissolved oxygen (µ mol m −3), minimum nitrate (µ mol m −3), minimum salinity (practical salinity scale), and minimum chlorophyll concentration (mg m −3). The layers summarize monthly conditions during 2000–2014 at 5-arcmin resolution (approximately 9.2 km at the equator). Values were retained in their native units, aligned to the common Bio-ORACLE grid, and clipped to the same coastal analysis mask; no additional mathematical transformation was applied before model fitting.

Delineation of coastal ecotones hosting oyster reefs

Oyster Reefs are primarily associated with coastal ecotones (the transitional interface between land and sea). Globally, due to the continuous nature of coastlines and the wide variation in the location of oyster reefs within estuaries or lagoons, it is impossible to precisely define the offshore distance threshold for oyster distribution at a local level. To address this challenge, we applied a standardized coastal buffer zone.

The spatial resolution of the environmental prediction layer provided by Bio-ORACLE is approximately 9.2 kilometers. Therefore, we defined our working area as all ocean cells within 45 kilometers of the global coastline, corresponding to approximately five grid cells.

This choice balances computational feasibility with ecological viability. In particular, it ensures that semi-enclosed systems such as estuaries, lagoons, and large bays are consistently included in the analysis, even if they extend tens of kilometers inland from the open coast. For example, the Beibu Gulf (located between the coasts of Tonkin (northern Vietnam) and South China) has extensive shallow-water environments suitable for oyster growth (Xu et al., 2020). The prediction layer was clipped to a 45 km mask using SDMtoolbox (Brown et al., 2017) to ensure that all environmental covariates were modeled with the same range and resolution.

Prediction of habitat suitability and distribution

To estimate the broad-scale environmental suitability of oyster reefs, we employed the Maximum Entropy (MaxEnt) model (version 3.4.4) (Phillips et al., 2017). MaxEnt is a widely used machine-learning model for predicting relative habitat suitability from presence-only and background data (Phillips et al., 2006; Elith et al., 2011; Merow et al., 2013). Its logistic out-put is interpreted here as relative suitability rather than a calibrated probability of reef occupancy. MaxEnt has previously been applied to regional oyster habitat-suitability modelling (Linhoss and Mickle, 2022); complementary field-validated habitat-suitability-index studies have assessed oyster performance under region-specific environmental conditions (Theuerkauf and Lipcius, 2016; Chowdhury et al., 2019). Here, we use MaxEnt as a macroecological screening tool to identify broad patterns of relative suitability across global coastal ecotones. The use of MaxEnt to characterize habitat-forming aggregative processes was introduced by Convertino et al. (2013). Our goal is not to mechanistically reproduce reef dynamics or to prescribe a local restoration species, but to identify coastal cells whose globally avail-able environmental conditions resemble those associated with pooled Ostreidae occurrence records.

The habitat suitability (HS) (note that here we adopt a macroecological characterization of the probability (Convertino et al., 2013), i.e. the average suitability of a community to host oyster reefs) H S = P (y = 1 | c (t)) is calculated as:

where f (c) is the probability density function of covariates c, and η (c) = α + ρ h (c). α is a normalizing constant that ensures that f 1 (c) integrates to one (f 1 being the probability density function (pdf) of the reef occurrences), and ρ is the constant (Lagrangian multiplier) of the MaxEnt features h (c) (Phillips et al., 2006; Phillips and Miroslav, 2008; Elith et al., 2011; Merow et al., 2013). The Lagrangian multiplier, that multiplies all environmental features, finds the optimal trade-off between model complexity (defined as the number of environmental variables used as predictors) and model accuracy (that is, the probabilistic distance between predictions and data) (Convertino et al., 2013; Merow et al., 2013). The set of parameters is identified by minimizing the prediction error between the observed and modeled floods. Features are transformations of the covariates in the covariate space, and this allows a faster and more precise computation rather than operating in the geographical space (Elith et al., 2011; Convertino et al., 2013).

Model calibration was conducted using the ENMeval package in R (Muscarella et al., 2014). To reduce overfitting, we tested combinations of six feature-class settings (L, LQ, H, LQH, LQHP, and LQHPT) and regularization multipliers of 1.0–5.0 (Radosavljevic and Anderson, 2014). Seventy percent of occurrence records were used for training, with four-fold cross-validation used for internal evaluation. Candidate models were compared with Akaike’s Information Criterion corrected for small sample sizes (AICc), and cross-validated area under the receiver-operating-characteristic curve (AUC) summarized discrimination between presence and background samples. The AUC is therefore an internal model-evaluation metric, not independent field validation of reef occurrence or restoration success (Roberts et al., 2017).

For ecosystem service quantification, we focused on the top quantile of HS values, representing highly suitable areas. Specifically, we defined

The continuous HS surface was classified into five classes using Jenks natural breaks. The lower boundary of the highest-suitability class was 0.6439; the rounded value 0.65 was retained as a transparent post-processing rule for priority mapping. It is not a universal biological occupancy boundary and does not represent the upper 20% of all predictions. Published marine applications use different purpose-dependent criteria, including occurrence-percentile, discrimination-based, fixed, and natural-break classifications (Freeman and Moisen, 2008; Liu et al., 2013; Davies and Guinotte, 2011; Guinotte and Davies, 2014; Yesson et al., 2017; Adamu et al., 2024; Chi et al., 2025).

We assessed threshold sensitivity by applying fixed values of 0.50–0.80 at 0.05 intervals to the same continuous HS raster and valid global coastal domain. At every threshold, the same geodesic pixel-area calculation and September 2023 World Database on Protected Areas (WDPA) overlay were used. We retained protected areas with a marine or coastal designation and a status of Designated, Established, Inscribed, or Adopted. The protected fraction was the high-HS area overlapping these sites divided by the total area exceeding the corresponding threshold; therefore, only the HS threshold varied among scenarios. Results are reported in Supplementary Figure S1 and Supplementary Table S1.

Habitat area calculation and country-specific species inference

We calculated the extent of highly suitable oyster reef habitats (H S ≥ 0.65) using the Tabulate Area function in ArcGIS. This tool calculates the surface area of raster cells that meet specified criteria within each geographic zone. Formally, for each zone z and suitability class c, the area is calculated using Equation 3:

where a i is the surface area of raster cell i, and I (·) is an indicator function equal to 1 if cell i belongs to zone z and class c, and 0 otherwise. In this study, zones were defined as (i) the global coastline mask, (ii) the exclusive economic zones of the ten hotspot countries, and (iii) the World Database on Protected Areas (WDPA) (UNEP-WCMC and IUCN, 2023). The suitability class was defined as all cells with HS ≥ 0.65, consistent with the threshold introduced in Equation 2. This approach allowed us to calculate global totals, country-level areas, and the overlap of highly suitable habitats with existing marine protected areas.

In addition, we identified the most frequently reported taxon name in the occurrence records for each hotspot country to provide descriptive taxonomic context for Figure 3.

Ecosystem service quantification: shell carbon and water filtration

We quantified two standardized physical indicators for comparing hotspot countries: the shell-carbon component associated with an assumed standing oyster population and its daily water-filtration capacity. The common parameter set was retained as a transparent global baseline because harmonized species- and region-specific measurements are not available for every hotspot. Combining isolated local estimates obtained for different species, size classes, seasons, reef conditions, and protocols would reduce cross-country comparability and could create false precision. These inputs are therefore scenario parameters, not universal biological constants. We assumed 15 oysters m −2 for 30% of reef area and 10 oysters m −2 for 70% of reef area, following the standardized reef-density framework described for Chesapeake Bay restoration assessment (Oyster Metrics Workgroup, 2011); the area-weighted value (11.5 oysters m −2) was rounded to B = 12 oysters m −2 in the global scenario. We further assumed a mean shell weight of 100 g, of which 12 g (12%) is elemental carbon, and an individual filtration capacity of 189 L day −1. Actual filtration varies with oyster size, temperature, salinity, food concentration, suspended solids, and measurement protocol (Ehrich and Harris, 2015).

Based on these parameters, we calculated the number of oysters associated with high-HS area (HS ≥ 0.65) in each hotspot country and estimated the corresponding shell-carbon component and filtration capacity. The total shell-carbon component C r was defined as:

and the total potential water filtration W r as:

where B is the standardized oyster-density input (oysters m −2), A (H S t o p) is high-HS area (Equation 2) expressed in m 2, C is shell carbon expressed in metric tonnes C oyster −1, and W is daily filtration (L oyster −1 day −1). Areas tabulated in km 2 were converted before density-based scaling using 1 km 2 = 10 6 m 2; consequently, C r is in metric tonnes C and W r is in L day −1. Here, C r is the carbon contained in the assumed shell component; it is not a complete net carbon budget, a verified carbon credit, an annual sequestration rate, or direct net atmospheric CO 2 removal. A net budget would need to resolve calcification-related CO 2 release, organic-carbon burial, shell persistence, sediment dynamics, and an appropriate counterfactual (Fodrie et al., 2017). Potential effects on adjacent seagrass beds, mangroves, and salt marshes are discussed conceptually but are not included in Equation 4 or in the numerical carbon valuation.

Economic valuation of oyster reef services

We evaluated monetary scenarios for three components: (i) oyster trade, (ii) the quantified shell-carbon component, and (iii) water filtration and nutrient removal. Oyster market values were derived from the country-level export-price inputs reported in SELINAWAMUCII (https://www.selinawamucii.com/insights/prices/). Reported prices ranged from $6.08 kg −1 in Italy to $227.50 kg −1 in Chile. For each hotspot country, the potential trade value E r was calculated as:

where 〈 b 〉 = 0.1134 kg oyster −1 is the mean individual mass used in the trade scenario, E is the country-specific unit export price, and division by 2 is the original conservative adjustment for uncertain restoration feasibility. The 12-g value above is shell carbon per oyster and is not the trade-mass input. Table 1 applies Equation 6 to large national high-HS areas after converting km 2 to m 2. It is an illustrative static gross standing-biomass price scenario, not annual revenue, profit, net present value, or a guaranteed financial benefit. It does not imply that every mapped high-HS cell can or should be restored or harvested.

For the shell-carbon valuation scenario, we multiplied C r (Equation 4) by a country-specific social-cost-of-carbon input derived from published country-level estimates (Ricke et al., 2018). In Equation 7, E c denotes the monetary coefficient expressed on the same mass basis as C r; when an input is reported per tonne CO 2, the equivalent per-tonne-C coefficient is obtained with the stoichiometric factor 44 / 12. The resulting scenario value E c r was defined as:

where halving again represents the original feasibility adjustment. This is a scenario-based valuation of the quantified shell-carbon component; it is not a carbon-market price, a verified offset, or a valuation of indirect carbon storage in adjacent vegetated habitats.

For water filtration, we considered both direct purification and nutrient removal services. The potential water filtration value E w r was defined as:

where E w is the filtration-service valuation coefficient used in the original scenario (US$ L −1 on the daily filtration-capacity basis), so E w r retains the time basis of W r. Similarly, nutrient-removal services were valued as:

where N is the dimensionless standardized nitrogen-removal scenario multiplier, E n is the area-based economic value of nitrogen removal, taken as the mean ($4,050.5 ha −1 yr −1) of the published range (Grabowski et al., 2012), and A (H S t o p) is expressed in hectares for this equation (1 km 2 = 100 ha). Thus, E n r is reported in US$ yr −1.

Together, these calculations retain the original GEP valuation framework (Ouyang et al., 2020). The components have different accounting and temporal boundaries: Table 1 is a static standing-biomass price scenario, whereas the nitrogen-removal coefficient is annual. Stock-based and annual-flow values are therefore reported as separate scenario components and should not be summed without temporal harmonization. Realized outcomes require local assessment of substrate, hydrodynamics, larval supply, mortality, time to maturity, permitting, protected-area objectives, competing coastal uses, restoration costs, governance, funding, and market demand.

Oyster reef habitat ecopotential for conservation and restoration

Using MaxEnt, we modelled global relative environmental suitability associated with pooled Ostreidae occurrence records (Figure 4). The selected model yielded an AUC of 0.90. This value indicates high internal discrimination between the available occurrence and background samples under the cross-validation design; it is not independent field validation, a calibrated probability of reef occupancy, or evidence of restoration success. We therefore use the output as a relative global screening layer (Roberts et al., 2017).

Analysis of variable contributions revealed that light availability at the seafloor was the dominant predictor, contributing over 40% to model performance (Figure 5). Minimum primary productivity and temperature range were the next most important factors, followed by maximum phytoplankton. These findings emphasize that oyster reef distribution is tightly linked to benthic light regimes, productivity gradients, and thermal variability, all of which are fundamental drivers of coastal ecosystem dynamics. Response curves further highlighted thresholds of tolerance to temperature and dissolved oxygen, reflecting the physiological constraints of oysters in sustaining viable reef structures (Figure 6).

The global pattern of predicted suitability shows strong concentration in ten hotspot nations: the United States, Australia, Chile, China, Japan, South Korea, Argentina, France, Germany, and Italy (Figure 3). Across these countries, high-HS area (HS ≥ 0.65) totalled 295,691 km 2, representing about 6% of their delineated coastal zones. Only 10.04% of that country-domain high-HS area coincided with cells containing an archived occurrence record (Figure 7). Because record collection is presence-only and spatially uneven, the remaining cells cannot be classified as unoccupied or interpreted as evidence that recorded reefs persist in suboptimal habitat.

Global oyster reefs have declined by more than 85% over the past century (Beck et al., 2011), but the present point archive cannot distinguish historical loss, non-detection, unsuitable substrate, or true absence in a cell without a record. Such high-HS cells are candidates for subsequent survey and feasibility assessment rather than confirmed vacant restoration sites. The separate, globally consistent September 2023 WDPA analysis estimated that 168,592 of 744,113 km 2 above the reference threshold overlapped the evaluated marine/coastal protected-area network (22.66%), leaving 77.34% outside it. This limited representation supports further local assessment of candidate areas, while any proposal to expand or redesign protection must consider verified reef presence, tenure, ecological condition, feasibility, and governance.

The threshold-sensitivity analysis supports this directional conservation-gap conclusion. As the cutoff increased from 0.50 to 0.80, global high-HS area declined from 1,428,737 to 289,493 km 2, whereas the protected fraction changed from 20.16% to 26.47% (Supplementary Figure S1, Supplementary Table S1). Thus, absolute area and quantities scaled directly from area are threshold-dependent, but 73.53–79.84% of high-HS habitat remained outside the evaluated network throughout the tested range.

Our global predictor ranking differs in emphasis from field-validated regional oyster habitat-suitability studies, which commonly identify salinity, temperature, dissolved oxygen, depth, substrate, and hydrodynamic setting as important controls on oyster performance and restora-tion success (Theuerkauf and Lipcius, 2016; Chowdhury et al., 2019; Linhoss and Mickle, 2022). In contrast, seafloor light and minimum primary productivity contributed strongly to the present global model. This difference likely reflects the broader spatial grain, the pooled family-level occurrence proxy, and the restricted set of globally consistent predictors; it does not imply that local physiological and geomorphic controls are unimportant. The global layer should therefore be combined with local substrate, salinity, oxygen, flow, recruitment, disease, and engineering data before restoration decisions are made.

Top-ten countries for oyster suitability and hotspots geomorphic characterization

For descriptive context, Figure 3 displays the most frequently returned taxon-name string in the occurrence archive for each hotspot country: Ostrea chilensis (Chile), Ostrea puelchana (Argentina), Ostrea edulis (France and Italy), Crassostrea virginica (United States), Magallana sikamea (China), Magallana gigas (South Korea and Germany), and Ostrea angasi (Australia). Because the retained occurrence dataset was not subjected to complete tax-onomic reconciliation and sampling effort was uneven across regions, these labels do not demonstrate that the reported taxa are locally dominant, native, or locally adapted, nor do they indicate species-specific habitat suitability. The selection of oyster species and stocks for restoration must therefore be evaluated locally and is beyond the scope of this family-level model. But these associations highlight that hotspot suitability is not only a function of environmental conditions but also of locally adapted taxa. Incorporating dominant native species into restoration designs is critical for long-term ecological persistence and functional success (Howie and Bishop, 2021).

Beyond species identity, the geomorphological context of hotspots provides further insight into reef distribution patterns (Figure 8). Across the top ten countries, oyster reefs were concentrated in three main coastal geomorphic settings: (i) lagoons and bays, often partially enclosed and dominated by wind-induced waves; (ii) estuaries, characterized by extensive tidal channels and mixed tidal–fluvial dynamics; and (iii) river deltas, where freshwater inputs and sediment delivery dominate hydrodynamics. Globally, approximately 72% of oyster reefs were associated with estuaries, 21% with lagoons or bays (including human-modified systems such as the Venice Lagoon in Italy or artificial bays in China’s Greater Bay Area), and 7% with deltas.

These geomorphic categories are not mutually exclusive: many reef systems occupy transitional or composite environments, such as Apalachicola Bay (USA), southwest South Korea, and coastal Zhejiang (China). Hydrodynamic conditions in deltas, in particular, can at times exceed thresholds for larval settlement, underscoring the importance of matching reef restoration strategies to local flow regimes.

Shell-carbon component and water-filtration capacity

Under the standardized parameter set in Equations 4, 5, the United States, Australia, and Chile have the largest area-scaled indicators (Figure 9). Their estimated shell-carbon components are 4.92, 4.43, and 2.95 million metric tonnes C, respectively. These values represent the carbon contained in the assumed standing shell component under the common scenario; they are not observed stocks, annual sequestration rates, verified carbon credits, or estimates of net atmospheric CO 2 removal.

This accounting distinction is important because oyster reefs can be carbon sources or sinks depending on geomorphic setting and the balance among calcification, shell preservation, and organic-carbon burial (Fodrie et al., 2017). Possible reef-mediated effects on sediment retention, water clarity, and adjacent seagrass, mangrove, or salt-marsh carbon remain ecologically relevant, but they are conceptual pathways in this study and are not included in the reported numerical carbon values.

The same standardized scenario yields daily filtration capacities of 77.58 quadrillion L for the United States, 69.92 quadrillion L for Australia, and 46.58 quadrillion L for Chile. These are area-scaled capacities under a common individual rate, not measured national water-quality improvements; realized clearance varies with biological and environmental conditions (Ehrich and Harris, 2015). Meta-analysis of eastern-oyster restoration on the United States Atlantic and Gulf coasts nevertheless shows that restored reefs can enhance multiple services, including oyster production, habitat provision, nitrogen removal, and sediment properties (Smith et al., 2023). The comparison supports the direction of the service pathways while also emphasizing that realized rates cannot be transferred uniformly to every global high-HS cell.

Potential economic value of reef-based oyster farming: benefits and perils

The economic assessment of reef-based oyster farming serves a dual purpose: to highlight the potential benefits of restored habitats as sources of livelihood and trade, and to identify risks of ecological pressure arising from international demand concentrated in biodiversity hotspots. This framing also opens the possibility of ecological payoff schemes, whereby trade-related revenues are reinvested in local restoration programs, similar to discussions of virtual water trade in global markets (Dalin et al., 2012).

Table 1 gives the largest illustrative gross value for Chile (approximately US$6.35 trillion), followed by Argentina (US$2.74 trillion) and the United States (US$1.31 trillion); Italy has the smallest value among the ten countries. The large totals arise because the density and price assumptions are scaled over extensive national high-HS areas, with 1 km 2 = 10 6 m 2. They are static gross standing-biomass price scenarios, not forecasts of realized trade, annual revenue, profit, net present value, or funds that governments or communities will necessarily receive.

Latin America plays a substantial role in the global seafood industry, where export markets have historically influenced fisheries development (De las Mercedez Jaime Muniz et al., 2022). China is also the largest producer and consumer of farmed oysters (Peng et al., 2021). These market contexts help explain why price and mapped area jointly affect the cross-country scenario, but they do not establish that high-HS cells are biologically available, legally harvestable, or capable of supplying the assumed quantity without changing price or demand.

Accordingly, Table 1 is an upper-level screening scenario rather than a restoration or investment plan. Realized outcomes require site-level assessment of substrate, hydrodynamics, larval supply, mortality, maturity time, permitting, protected-area objectives, competing coastal uses, restoration and operating costs, governance, funding, and market demand. The table should therefore be read alongside, but not automatically added to, the separately bounded ecosystem-service scenarios.

Global ecosystem product of oyster reefs

Figure 10 retains the original monetary scenarios for the quantified shell-carbon component, water purification, and nitrogen removal across the ten hotspot countries. The United States and Australia have the largest area-scaled values. For the United States, the component values shown are approximately US$13.8 billion for nitrogen removal, US$2.5 billion for the shell-carbon scenario, and US$8.8 billion for water purification. Only the nitrogen-removal coefficient is explicitly expressed as an annual flow in the present framework; the components should therefore not be summed or interpreted as a single annual national return without temporal harmonization. They exclude restoration costs and foregone harvest opportunities (Grabowski et al., 2012).

Cross-country variation reflects both mapped area and the standardized economic inputs. In particular, the shell-carbon value uses a social-cost-of-carbon scenario, not a carbon-credit price; it does not demonstrate credit eligibility or market revenue. The values are intended to compare potential scale under common assumptions rather than predict locally realized benefits.

The direction of these service pathways is consistent with restoration syntheses, but realized recovery depends on project goals, site selection, reef design, monitoring, and environmental context (Smith et al., 2023; Smith and Pruett, 2025). Likewise, location-specific oyster-reef valuations vary with ecological production, beneficiaries, costs, and the spatial setting of the service (Grabowski et al., 2012). These comparisons support using GEP as a transparent global screening framework while reinforcing that site-specific biological, economic, and governance information is required before implementation.

Limitations and future directions

The principal limitations define how the maps and valuations should be used. First, GBIF and iNaturalist records may contain identification and coordinate errors and reflect uneven observation effort. The present workflow did not apply complete research-grade filtering, taxonomic reconciliation, automated outlier screening, coordinate deduplication, or spatial thinning. The family-level model is consequently appropriate for broad screening, not for choosing a local restoration species. Future versioned workflows should retain provider meta-data and apply taxonomic and coordinate quality control, bias-aware thinning or background selection, and species- or region-specific modelling where sufficient verified records are available (Beck et al., 2014; Boria et al., 2014; Aiello-Lammens et al., 2015; Zizka et al., 2019). A fully versioned end-to-end release should also archive the model code, software versions, and machine-readable parameter settings.

Second, the AUC of 0.90 is an internal presence–background discrimination measure. Spatially structured cross-validation and independent regional reef inventories or field observations are needed to evaluate transferability and restoration relevance (Roberts et al., 2017). The 5-arcmin Bio-ORACLE climatologies summarize 2000–2014 conditions and do not resolve substrate, larval connectivity, disease, fine-scale hydrodynamics, coastal engineering, or watershed-derived nutrients and pollutants. The model is static and does not forecast climate change, sea-level rise, marine heatwaves, deoxygenation, storm exposure, carbonate chemistry, land-use change, or future pollution. Coupling HS with projected oceanographic, sea-level, and watershed-stressor scenarios would allow future studies to identify persistent refugia, emerging opportunities, and areas where present suitability may decline.

Third, the 0.65 decision rule affects all area-scaled quantities. Across the tested global thresholds, protected representation remained comparatively stable (20.16–26.47%), but absolute high-HS area changed substantially (Supplementary Figure S1; Supplementary Table S1). The sensitivity analysis uses the global valid coastal domain, whereas the economic comparison uses the ten-country hotspot domain; these totals should not be interchanged. The continuous HS surface remains the primary model output.

Finally, the common oyster-density, shell-carbon, filtration, and valuation inputs support comparability but do not represent all species or regions. Future work should develop harmonized species- and region-specific parameter datasets and propagate distributions for occurrence error, model predictions, thresholds, ecological rates, prices, costs, and feasibility rather than relying only on point scenarios. It should also extend the framework to multi-objective conservation–production and mariculture planning; the water-use and carbon-neutrality optimization of Tan et al. (2025) provides one relevant example. These extensions are future objectives rather than capabilities claimed for the present climatological baseline.

Blue carbon ecosystems (BCEs)—including mangroves, seagrasses, and tidal marshes—are increasingly recognized as natural climate solutions due to their capacity to store carbon while providing co-benefits such as water purification and biodiversity support (Macreadie et al., 2021). In this study, we extended this perspective to oyster reefs, which function as ecosystem engineers and can interact with adjacent BCEs. We introduced ecopotential as the integration of continuous relative HS with standardized ecosystem-service scenarios, and GEP as the corresponding monetary valuation framework, to screen priority regions for further conservation and restoration assessment. By linking ecological patterns with explicit economic scenarios, the framework complements other blue-carbon and natural-capital approaches (Convertino and Valverde, 2013); it does not itself establish reef presence, restoration feasibility, net carbon removal, or realized financial return.

The analysis identifies where the standardized shell-carbon, water-filtration, nutrient-removal, and standing-biomass price scenarios are largest. These values are indicative comparisons, not predictions. Shell carbon is not net atmospheric CO 2 removal, conceptual effects on adjacent vegetated habitats are not numerically valued, and the Table 1 totals are not annual revenue or guaranteed benefits. Restoration outcomes will depend on local substrate, hydrodynamics, larval supply, species interactions, basin-scale stressors, costs, governance, and sustained investment.

Future research and policy efforts should focus on three directions:

Basin-scale frameworks. Oyster reefs should be considered as natural infrastructure at the land–sea interface, where stability and function depend on hydrological dynamics and reef structure. Quantifying how reef networks influence wider coastal systems is essential to secure their long-term functionality.

Integrated valuation of natural habitats. Beyond ecological metrics, oyster reef restoration should be assessed in terms of cascading benefits for other habitats, local economies, and climate regulation. Integration into national accounts or blue-carbon frameworks requires locally validated service rates, explicit counterfactuals, compatible temporal units, restoration costs, and a complete net carbon budget.

Climate adaptation. Scenario-specific climate, sea-level, oceanographic, and watershed forcing should be coupled to the present climatological baseline before it is used to infer future refugia or restoration trajectories. Multi-objective extensions could also evaluate conservation and sustainable mariculture jointly.

Finally, we stress that local restoration cannot succeed in isolation. Many suitable oyster reef sites are located in estuaries and deltas heavily affected by nutrient loading, pollutants, and basin-scale stressors. Effective conservation and restoration solutions must therefore combine ecological restoration with watershed or basin management, including pollution control, and education to reduce human pressures. Only through such integrated basin approaches can oyster reefs realize their full potential as ecological and economic assets, reinforcing coastal resilience and supporting global sustainability goals.