Ancient oysters provide far-field sea-level constraints on ice melt timing and sources during the early–mid Holocene climate transition

The early- to mid-Holocene represents a critical period when Earth’s climate system transitioned from rapid post-glacial warming to the relatively stable conditions that characterized much of human civilization¹. During this interval, ongoing deglaciation drove continued sea-level rise and reshaped coastal environments worldwide². Understanding the timing and sources of this ice loss is important for constraining ice sheet responses to warming and improving projections of future sea-level change³. Far-field sea-level records—those located distant from former ice sheets—provide valuable insights into this process because they primarily reflect changes in global ice volume rather than local deformation from ice unloading⁴. However, questions remain about when major deglaciation ended and which ice sheets contributed most significantly^5,6, with recent studies presenting conflicting evidence about whether ice loss ceased by 7–6 ka or continued until 4–2 ka^7,8.

Records from the South China Sea (SCS) have been used in the analysis of ice volume and global mean sea-level change^7,8. However, data constraints from this region remain debated, particularly regarding the existence of a relative sea-level (RSL) highstand^9,10,11, which refers to higher-than-present RSLs between 7 and 4 ka. These debates stem partly from poor preservation of sedimentary indicators typically used to constrain RSL on intermediate-field coastlines, which lack accommodation and are prone to reworking on coastlines under stable to falling RSL¹². This issue is compounded by recent syntheses^9,10,11 that rely on earlier compilations⁹ without incorporating subsequently published data and consequently include older records with greater uncertainties, which have not been subject to rigorous quality control protocols.

Limited evidence for a mid-Holocene highstand has been documented in the Pearl River Delta (PRD)^13,14, the most intensely debated of which comes from an ancient oyster bed originally reported to record a highstand of ~ 1.9 m^15,16,17. Its interpretation has been questioned due to uncertainties in modern analogues, potential influence from waves and storm surges, beach erosion, growth in tidal pools, imprecise surveying, and the lack of corroborating evidence^15,17. The literature also lacks detailed descriptions of the reef’s lithology, biology, and dating, with the original interpretation based on a single radiocarbon age¹⁴.

Here, we address the controversy surrounding mid-Holocene RSL change in the SCS by integrating new and existing data from the PRD (Fig. 1). Our approach combines ecological, morphologic, and geochronologic evidence from an ancient oyster reef with a region-wide, quality-controlled database. By comparing our findings with other far-field sites across Southeast Asia, Australia, the Mediterranean, Africa, and the Americas, we establish a broader context for evaluating ice-melt scenarios. This synthesis provides a robust test of competing hypotheses for Holocene RSL evolution and enables reassessment of the timing and sources of global ice-volume changes.

Environmental and climatic controls on the vertical distribution of Saccostrea cuccullata

S. cuccullata (Born, 1778) is a dominant rock oyster thriving across intertidal rocky shores in brackish to fully marine waters¹⁸. However, quantitative assessments of its vertical distribution across different tidal regimes are limited. Our regional surveys show that S. cuccullata predominantly encrusts sheltered rock faces oriented perpendicular to the shoreline and sheltered from the predominant waves, often forming interlocking, three-dimensional structured habitat (Fig. 2B, C). In contrast, wave-exposed rock surfaces parallel to the shoreline tend to be dominated by barnacles (Chthamalus malayensis, Capitulum mitella) and other sessile organisms better adapted to high-energy conditions¹⁹. The vertical distribution of S. cuccullata at Big Wave Bay (BWB) (25–190 SWLI; standardized water level index units, which normalize tidal elevation across different tidal regimes) matches other sites across Hong Kong (38–184 SWLI units, 2σ)^19,20 and the wider Indo-Pacific^21,22,23,24, where this species is found within the great diurnal range (Mean higher high water [MHHW]-Mean lower low water [MLLW]; 1.44 m in BWB) (Table S1; Fig. 2A). Abundance data confirm that its vertical range can be modelled by a normal distribution in statistical reconstructions of past RSL and its rates of change²⁵ (Fig. S1). This regional consistency in vertical habitat range responds similarly to tidal regimes reinforces the reliability of S. cuccullata as a sea-level indicator. Accordingly, we define a 2σ indicative range of 38–184 SWLI units (0.85–1.86 mPD; − 0.34–0.56 mMSL in BWB) for S. cuccullata (Fig. 3; Table 1), which is used to reconstruct RSL (“Methods” section, Fig. 6). In the absence of local survey data, a conservative indicative range of MHHW to MLLW may be adopted.

A concern is whether oysters colonize higher elevations within tidal pools or due to “splashing effects” from typhoons and monsoon influence^15,17. We observed S. cuccullata occupying higher elevations in tidal ponds at BWB (Fig. 2D, E), a phenomenon known as ‘ponding’ previously described for coral microatoll indicators²⁶. Ponding occurs when water trapped from ebbing tides in natural depressions along rocky shores form microhabitats that allow organisms tolerant to temperature and salinity fluctuations, like S. cuccullata, to survive above their typical upper limit on the shore²⁶. In Hong Kong, S. cuccullata commonly encrusts within tidal pools frequently replenished by strong wave action¹⁹. At BWB, oysters were found in tidal pools between 169 and 267 SWLI units (Fig. 2A), with 62% encrusted above MHHW (1.98 mPD; 200 SWLI units), but none above HAT (3.15 mPD; 368 SWLI units). Including ponded oysters broadens S. cuccullata indicative range to 38–267 SWLI units and raises reference water level (RWL) by 58% to 115 SWLI units (Fig. 3; Table 1). However, even with ponding, no oysters were observed above 2.45 mPD, which remains well below the base of the ancient oyster bed.

We further deemed episodic typhoon surges and monsoonal wave action insufficient to facilitate sustained oyster colonization above their typical range¹⁷. Hong Kong experiences 5–6 typhoons annually generating wave heights up to 10 m²⁷. Although sedimentary proxies from the northern SCS suggest increased tropical cyclone frequency and stronger monsoons during the mid-Holocene²⁸, these events are too infrequent for sustained colonization at higher elevations. Sessile organisms, like S. cuccullata, have low tolerance to desiccation stress and require consistent inundation (e.g., from wave splashing) to survive desiccation stress; when irrigation ceases, colonizers at higher elevations typically disappear within weeks²⁹. Factors such as reduced immersion, limited feeding time, and increased physiological stress further hinder sustained colonization above MHHW level³⁰. Typhoon-induced surges occasionally reach the height of the ancient reef, but their rarity (≥ 10.5-years return period recorded at tide gauge records from Tai Miu Wan and Waglan Island; Table S3) precludes sufficient exposure to allow for S. cuccullata survival above MHHW. Similarly, monsoon-driven winter wave action allows S. cuccullata to occur ~ 30 SWLI units higher during the dry season¹⁹ (Fig. S2; Table S2). However, even under peak monsoon conditions, unponded oysters remained below MHHW. If enhanced mid Holocene monsoon-driven wave exposure were sufficient to enable reef formation at the ancient oyster bed elevation, contemporary S. cuccullata would be expected to colonize similar heights during strong monsoon seasons. Instead, no S. cuccullata specimens (including those in tidal pools) were observed above 266 SWLI units (2.45 mPD) (Fig. 2A), which is still 53 SWLI units (0.96 m) below the base of the ancient oyster bed (Fig. 4). This gap between maximum modern occurrence and ancient reef elevation indicates that typhoon-driven and seasonal wave splashing, even under enhanced monsoon conditions, cannot explain S. cuccullata formation above its modern range.

In situ formation of the ancient oyster reef: morphological, geochemical, and chronological evidence

Since the discovery of the ancient oyster bed in BWB¹⁴, its origin has been debated, particularly whether extreme events have transported and reworked the deposit to elevations above its living counterparts^15,16,17. To address these questions, we analysed the reef’s morphology, structure, composition, and chronology to determine whether it formed in situ under conditions comparable to today.

Petrographic evidence indicates that the ancient oyster bed is a filled frame reef³¹ (Text S1; Fig. S3), constructed in place by S. cuccullata skeletons. The domical framework shows upward and lateral growth bands from vertical and lateral accretion (Fig. 4B, 5G), with semi-open cavities infilled by detrital material, providing protection from wave damage and bioerosion. Limited cementation (< 20%) (Fig. 5C–G) includes marl within oyster cavities and crystalline calcium carbonates between intact skeletons (Fig. 5H, I), suggesting periodic wave splash during formation. The matrix contains subangular (roundness: 0.52 ± 0.24), poorly sorted sediment grains (Fig. S4) of 68% calcite from biolithic fragments (e.g., bivalves, coralline algae, and reef skeletons) damaged by wave activity and 32% quartz typical of mature coastal sediments in Hong Kong (Fig. 5K). These characteristics differ from local beach sands (higher roundness of ~ 0.8)³², suggesting that detrital components were transported by wave action but shielded from further abrasion and sorting once trapped within reef cavities³¹. Collectively, these features confirm in situ reef development under dynamic wave conditions near the beach face during higher-than-present RSL.

The hypothesis that the reef was deposited at higher elevations due to storm activity is inconsistent with the evidence. The BWB reef has been compared to storm-generated tempestites found as high as 20 m above present sea level¹⁵. Reefs reworked during extreme wave events typically exhibit tilted, overturned, or irregular orientations and incorporate reef debris from different ages³³. Multiple lines of evidence instead support in-situ reef formation. First, the vertically-oriented calcite lamellae of S. cuccullata skeletons (Fig. 5D, G) and flat base of the BWB reef indicate in-situ growth on a stable substrate, unlike storm-reworked reef boulders that display irregular shapes due to sliding or breakage³⁴. Although the reef’s upper surface show signs of wave- and bio-erosion³⁴, its preserved basal morphology confirms original depositional position. Second, six radiocarbon ages sampled across the reef support prolonged formation: older ages occur in the inner and upper portions and closely spaced samples (e.g., BWB_OYSTER_04 & 05) yield concordant ages (Fig. 4D, E; Table 2). Oyster reefs develop through layered lateral and vertical growth within their living range. The agreement among ages indicates systematic lateral accretion from 5.4 to 3.3 ka, tracking falling RSL. This ordered growth pattern contrasts with storm deposits, which would show randomly mixed fossil ages throughout the reef. Third, petrographic analysis of proposed tempestite deposits¹⁵ revealed rounded pebbles and higher degrees of abrasion, features absent in the BWB reef, further ruling out storm-related deposition.

Finally, it is unlikely that the ancient reef formed within a tidal pool, given its morphology and depositional context within the bedrock cliff. Unlike contemporary ponded oysters observed in elevated tidal pools occupying depressions on rock platforms (Fig. 2D), the ancient reef is situated on a concave slope covered by an overhanging slab (Fig. 4A). This slope geometry would prevent water retention, unless blocked by a large obstacle. Such an obstacle would have left geomorphic evidence, for example a large boulder, which is absent at the outcrop location. Furthermore, if the reef had formed within a tidal pool, it would consist of patchy oyster clusters (Fig. 2D, E) rather than a continuous reef. Instead, the concave shape of the slope, sheltered by the overlying rhyolite, likely provided a low-energy, stable environment that enabled development of a much more extensive reef structure than the present-day oyster encrustations¹⁶.

Timing and magnitude of mid-Holocene highstand in the PRD

The timing and magnitude of the mid-Holocene RSL highstand in the PRD remains debated. Although previous studies have proposed a mid-Holocene highstand in the PRD of at least 2 m between 5.9 and 4 ka^13,14 with higher levels potentially persisting up to 1.2 m until 1.6 ka^35,36. However, these indicators were either omitted from earlier RSL compilations without explanation^9,10,11 or dismissed due to potential reworking by storm processes¹⁵, which our analyses have ruled out.

To address the limitations of existing compilations, which lack rigorous quality control and omit published data, we updated and reinterpreted the available RSL data for Hong Kong and the PRD region^9,10,11 (Fig. 1B) from 11 ka to present following international standard protocol (“Synthesis and modeling of RSL data” section). L1 data satisfy all quality criteria, including reliable geographic location, calibrated age (ka), elevation, and well-constrained indicative meaning (Table S4 & S5), whereas L2 data may be affected by unresolved uncertainties, including potential errors not fully captured in the uncertainty estimates. This database contains 48 sea-level index points (SLIPs), which delimit the past position of RSL with corresponding age and elevational uncertainties and 244 terrestrial and marine limiting data points that provide upper or lower RSL constraints, respectively. Most data are intercalated sedimentary samples (Fig. 6, Supplemental information 2), complemented by six new SLIPs derived from the BWB ancient oyster reef. We applied a temporal empirical hierarchical model (ET-GP)²⁵ to this standardized, quality-controlled RSL data compilation. Our analysis confirms that a mid-Holocene highstand was achieved in the region by at least 5.6–5.3 ka, with RSL reaching + 2.3 ± 1.2 m (2σ) above present-day levels. A notable temporal gap exists in the record between 7 and 6 ka, with no SLIPs available to constrain RSL during this interval. Consequently, we do not attempt to define the precise timing or peak magnitude of the highstand; rather, the value of + 2.3 ± 1.2 m at 5.6–5.3 ka represents a minimum magnitude and timing estimate, with the ET-GP reconstruction reflecting this data gap through characteristically wider uncertainties in that interval. After the highstand, RSL declined gradually to + 1.6 ± 1.8 m between 3.5 and 3.2 ka due to hydro-isostatic adjustment³⁷ and equatorial ocean siphoning³⁸, before stabilizing at present-day levels between 1.2 and 0.5 ka.

Additional evidence for higher-than-present RSL in the broader PRD corroborates the BWB ancient oyster reef, although with lower reliability. A single Corbicula sp. shell encrusted on a tidal notch at Shijie¹³ (Fig. 6B) provides the most notable supporting evidence, indicating RSL reached + 3.6 ± 1.2 m above present between 5.9 and 4.5 ka, which overlaps with the age and RSL indicated by the ancient oyster at BWB. However, several factors constrain the reliability of the Shijie evidence. The genus Corbicula consists of both fresh and brackish water species³⁹. Without species-level identification or sufficient understanding of the genus’s modern vertical distribution, precise RSL interpretation remains uncertain. The specimen’s presence on a hard surface also deviates from its typical burrowing behaviour in soft substrates³⁹, adding further uncertainty. Elevated beach rock in the PRD and Hong Kong^35,36 (Fig. 1B) also supports the BWB oyster reef data, suggesting RSL was 0.9 m higher at 2.9 ka (Fig. 6B). However, these beach rock samples are considered lower-quality indicators (L2 ranking; Table S4 & S5) because they lack diagnostic supporting evidence (e.g., seaward-dipping planar laminations; calcareous cement observed in petrographic analysis)³⁵ and their ¹⁴C ages were obtained from shells that may potentially be reworked³⁵. Only two other low-elevation L2 sedimentary indicators (Fig. S5) contradict the highstand evidence, but these intercalated mangrove and tidal flat indicators^9,36 are missing supporting microfossil, geochemical, and sedimentary evidence and are subject to potential influence from compaction. The spatial coherence among the Shijie outcrop, PRD and Hong Kong beach rock, and BWB reef strengthens the regional evidence for a mid-Holocene highstand in the PRD.

The BWB ancient oyster reef underscores the importance of fixed biological indicators in defining highstand magnitude in far-field locations like the PRD. Fixed biological indicators constrain RSL exclusively along the falling limb during the mid to late Holocene, with no evidence of their formation during the early-Holocene rise (Fig. 6A). This absence reflects the rapid early-Holocene rise (up to 30 m ka⁻¹)^9,10, which outpaced the vertical accretion capacity of oyster reefs (5–39 m ka⁻¹, mean ≈ 10 m ka⁻¹)⁴⁰ and enhanced intertidal erosion, preventing the stability required for reef development. As RSL decelerated between 8 and 6 ka, oyster reef accretion became possible. Sedimentary indicators bracket the highstand through terrestrial (e.g., freshwater-swamp deposits) and marine (marine deposits) limiting constraints (Fig. 6), but do not produce more precise SLIPs because intertidal facies are often reworked by fluvial erosion during subsequent progradation⁴¹. Our findings highlight how fixed biological indicators compensate for the limitations of sedimentary indicators and fill temporal gaps during the mid to late Holocene under a regime of falling RSL in far-field settings⁶.

Implications for the source and timing of deglacial barystatic input

As a far-field location, the SCS coast is central to estimating the barystatic signal from RSL records through glacial isostatic adjustment (GIA) modelling⁷. Using our RSL data compilation and GIA model runs, we provide new constraints on early- to mid-Holocene RSL and barystatic changes, the timing and magnitude of which remain poorly understood due to chronological uncertainties and limited geological evidence for paleo-ice sheet coverage in near-field regions^7,42. We first assessed how well two widely used independently-developed models—ANU-ICE and ICE-6G_C—capture early- to mid-Holocene barystatic signals. This analysis reveals substantial misfits between RSL records and GIA predictions during the early Holocene, with the ANU-ICE pairing achieving a better fit to RSL data compared to ICE-6G_C (Fig. 7A).

To provide context for the factors contributing to the better fit of the ANU-ICE pairing, we conducted sensitivity analyses examining both ice-equivalent functions and Earth rheological structure (Fig. 7B, C). When combining ICE-6G_C with four different solid Earth models with different mantle viscosity structures, data-model misfits remain consistently high (Fig. 7E). Thus, varying Earth models alone does not improve the fit with RSL estimates in this region (Fig. 7D). Varying ice models, however, reveals that model pairings using ANU-ICE and ICE-6G_C with a delayed global barystatic input exhibit lower misfits than the rest of the ICE-6G_C model runs (Fig. 7D, E). These better-performing ice models share a delayed timeline of ice melt deceleration and termination, producing lower RSL predictions during the early Holocene and improved agreement in mid-Holocene highstand timing and magnitude. During this period, RSL rise was decelerating due to the final deglaciation of the Scandinavian (~ 9.5 ka)⁴³ and Laurentide Ice Sheets (~ 7 ka)⁴⁴ and reduced barystatic input from the Greenland and Antarctic Ice Sheets^7,42.

While a delay in the barystatic signal improves the fit, regional processes like fault activity and basin extension^10,45 could lower early Holocene SLIPs in the PRD by 1–3 m, with further contributions from sediment compaction¹⁰, leading to the misfit with non-delayed ICE-6G_C scenarios. Further investigation is needed to disentangle local and global signals across the SCS. However, a 1 ka delayed global contribution better fits far-field sites across Asia, Oceania, Africa, and South America^{46,47,48,49,50,51,52,53,54} (Fig. 8). Exceptions include St. Croix, where deformation from the collapse of the Laurentide Ice Sheet forebulge and proximity to Barbados, a site to which the ice-equivalent function of ICE-6G_C is tuned⁵⁵, favour this model pairing. Tunisia, lacking early-Holocene (pre-7.5 ka) data, shows similar fits across all models for mid- to late-Holocene RSL records despite additional water loading effects modulated by the coastal geometry within the Gulf of Gabès⁵⁶. Overall, the consistent preference for a 1 ka delayed contribution across far-field sites strongly supports its role as a key factor in resolving early Holocene RSL discrepancies.

Comparing the two best-fit models shows that ANU-ICE performs slightly better than ICE-6G_C with 1 ka delayed global deglaciation (Fig. 7A, D), suggesting that differences in total global meltwater (58 vs. 49 m sea-level equivalent [SLE] at 11.7 ka) may influence model performance^7,55,57. Further sensitivity tests demonstrate that delaying deglaciation and increasing barystatic input from Antarctica alone does not substantially improve model-data misfit (Fig. 7C, D), indicating that additional Northern Hemisphere ice contributions are required. Recent estimates of Holocene ice volume budgets suggest Antarctic contributions of 5.3–10.2 m SLE⁵⁸, North American ice sheets 33–53 m SLE⁵⁹, and Eurasian ice sheets 2.1–4.3 m SLE^5,43,60, totalling 40–66 m SLE during the Holocene⁵. These estimates align with independent reconstructions of total Holocene global mean sea-level change of 45.1 m (35.4–54.6 m 2σ range)⁵ and ~ 43 m (33–53 m, 90% credible interval) with ~ 11 m (4–17 m SLE) from Antarctica⁸. Ice sheet reconstructions consistently point toward dominant Northern Hemisphere contributions (~ 80%) relative to Antarctica (~ 15%), contrasting with ANU-ICE’s roughly equal Antarctic (~ 26 m SLE) and North American (~ 26 m SLE) contributions⁶. Integrating evidence from ice sheet reconstructions with our finding that Antarctic-only modifications fail to improve model fit, these results support a relatively modest Antarctic contribution and favor delayed deglaciation from larger Northern Hemisphere ice volumes—particularly the North American ice sheet complex—during the early to mid-Holocene.

This study advances knowledge of mid-Holocene RSL change in the SCS and provides new constraints on global deglaciation patterns during the Holocene climate transition. By establishing a quantitative indicative meaning for S. cuccullata (111 ± 73 SWLI units) and demonstrating its ecological validity through species distribution analysis, we expand the toolkit for RSL reconstruction in far-field coastal environments where traditional sedimentary indicators may be absent due to falling RSL.

Our morphological, petrographic, and mineralogical analyses identified the ancient oyster bed as a filled-frame reef that formed in situ under wave action, ruling out alternative formation hypotheses related to intensified wave action, more frequent typhoons, or a stronger monsoon. Six new SLIPs fill a temporal gap in the PRD RSL record and help to constrain the timing and magnitude of the mid-Holocene highstand, showing RSL reached its peak by at least 5.6 ka. These results are consistent with other regional highstand estimates, confirming a sustained period of elevated RSL during the mid-Holocene.

Comparison of the PRD data with GIA model predictions reveals better agreement with delayed deglaciation scenarios, supporting growing evidence for protracted Northern Hemisphere ice sheet decay and providing insights into ice sheet-climate feedbacks and the timing of barystatic contributions to global sea-level rise. Future integration of PRD records into broader regional syntheses will be essential for disentangling local processes from regional and global signals and enhancing constraints on ice sheet contributions to past RSL change.

Site selection

BWB is situated on the southeastern coast of Hong Kong Island in Hong Kong (Fig. 1B, C), within a tectonically stable setting with no evidence of large-scale vertical displacement from coseismic events during the Holocene^61,62. The region features a wet season from June to September when it is under influence of the warm and moist oceanic air from the Pacific, and a dry season from December to March when cool and dry air from the Siberian anticyclone in the north persists. BWB is characterised by a wave-dominated environment, with a barrier system flanked by rocky headlands and cliffs. The east-facing beach at the centre of the bay is bordered by a beach ridge that separates it from a small creek flowing into the embayment (Fig. 1B). Human infrastructure is located on and behind the ridge, where the local geomorphology has been altered. The nearest tide gauge at Tai Miu Wan indicates local tides are semidiurnal, with a great diurnal range from mean higher high water to mean lower low water (MHHW-MLLW) of 1.44 m.

The oyster bed remains (S. cuccullata) are preserved within a low cliff along the southern side of the beach (Fig. 1C). This deposit contains a substantial accumulation of shells attached to a rhyolite outcrop, shielded from erosion by an overhanging rock slab. The present environment at BWB is a moderately exposed rocky shoreline, situated above the intertidal zone. Unlike the dense ancient bed, modern S. cuccullata specimens appear as individual encrustations on rock surfaces rather than forming thick accumulations. Other species found along the rocky shoreline include the gastropod Echinolittorina malaccana in the supratidal fringe and the barnacle Chthamalus malayensis spanning the intertidal zone¹⁴.

Contemporary distribution of Saccostrea cuccullata

We conducted a survey of contemporary S. cuccullata living near the oyster remains. Our survey included two different S. cuccullata habitats: 1) rock surfaces in the intertidal zone regularly exposed to daily tidal inundation and wave action, and 2) tidal pools formed from water trapped in rock benches from breaking waves at the highest tides (Fig. 1D; Table S1; Supplementary information 3). The upper and lower bounds of the encrusting oyster bands at each survey location were measured using a Leica GS18 global navigation satellite system (GNSS) to record their elevation relative to the Hong Kong Principal Datum (PD) (1.3 m below mean sea level [MSL]), with an accuracy of ± 5 cm.

In addition, we used ecological survey data to examine the vertical living range of S. cuccullata across Hong Kong (Fig. S2; Table S2). The dataset documented the percentage cover of S. cuccullata (expressed as relative abundance) from transects at sixteen coastal locations across two years. At each location, six transects were established across the intertidal zone. The upper and lower bound of the transects were defined by biological effective tidal height³⁰. The upper end of the transect was positioned at the tidal height where the splash zone gastropod Echinolittorina malaccana was found, while the lower end was set 25 cm below the upper limit of crustose coralline algae, or at 1 m above Principal Datum in the absence of crustose coralline algae. The elevation of transects was measured using a South Survey Galaxy G1 Plus GNSS. Photographs were taken of ten 25 × 25 cm quadrats placed randomly along each transect, and the percentage cover of S. cuccullata was estimated using the Coral Point Count software. The relative abundance at each elevation was averaged for comparison across sites.

To allow comparison of the vertical distribution of oysters across different tidal regimes in Hong Kong, where tidal range varies from 1.03 to 2.11 m (Fig. S2), we expressed the elevation of modern living oysters using standardised water level index (SWLI) units⁶³:

where a value of 100 corresponds to mean tide level (MTL) and a value of 200 corresponds to MHHW. Tidal datums at each site were estimated from the nearest tide gauge (Fig. 3).

We quantified vertical distribution patterns of S. cuccullata from the regional distribution data using skewness and Quantile–Quantile plots within the pingouin package in Python⁶⁴.

Surveying and sampling the oyster reef remains

We scanned the ancient oyster bed with an Emesent Hovermap ST light detection and ranging (LiDAR) system (Fig. 4C–E), with an accuracy of ± 2 cm. The point clouds were pre-processed using the Emesent Aura software and georeferenced to four ground control points that were tied into Principal Datum using a Leica GS18 GNSS system. Digital elevation models (DEM) of the reef remains were constructed using the georeferenced point clouds and the lidR⁶⁵ package in R.

We collected six samples from the ancient oyster bed for radiocarbon dating using a hammer and chisel in February 2024 (Fig. 4A). The elevation of the samples were surveyed by a Topcon OS-200 total station from a temporary benchmark georeferenced by the Leica GNSS system. An additional 15 cm × 8 cm × 5 cm block was collected for petrographic and mineralogical analysis. After collection, the samples were soaked in freshwater and cleaned by brushing to remove microbial mats in preparation for further analyses.

Sedimentary petrology and mineralogy of the oyster reef remains

To characterize the ancient oyster bed and interpret its depositional context, we described its macro-constituents and performed petrographic and mineralogical analyses. Half of the 15 cm × 8 cm × 5 cm block was fractured for examination of its macro-constituents (Fig. 4A). Oyster species were identified following the taxonomy outlined in local survey¹⁹. We prepared five thin sections—three from the side parallel to the shore (hereafter referred to as the ‘lateral’ axis), and two from the side perpendicular to the shore (the ‘horizontal’ axis) (Fig. 5A–B). Each 30-µm thin section was scanned using a uSCOPE MXII microscope under cross-polarised light to describe shell structure and orientation and the biolithic constituents forming the reef matrix. We used ImageJ to assess the roundness and circularity of sediment grains in the reef matrix⁶⁶.

The term ‘reef’ describes calcareous deposits essentially built in place by sessile organisms, creating topographic relief higher than surrounding sediment or rock substrate⁶⁷. These aquatic biosedimentary structures are formed through sediment binding and baffling, and the biological activities of the aquatic organisms that construct them⁶⁸. We followed the reef classification system of Riding³¹ to classify the ancient reef and evaluate evidence of its formation. This system categorizes reefs into three types—cluster, segment, and frame—based on the relative composition of its matrix, skeleton, and cement (MSC) components. The skeleton component consists of the surficial biomineralized builder of the reef, which can vary in morphology, orientation, attachment, and support³¹. The skeletons may be supported by a matrix derived from particulate sediment and/or by cement, which precipitates from water in pores partially enclosed within or between skeletons and matrix grains. The proportion of the three components can be represented on a MSC ternary plot to classify the reef type (Text S1; Fig. S3). This approach allows for the identification of depositional and constructional processes that shaped the reef.

We also conducted X-ray diffraction (XRD) on the oyster shells and reef matrix within the ancient reef to examine the mineral composition of oyster reef remains. Both the shell and matrix samples were ground and homogenised using a mortar and pestle. Their powder X-ray diffraction (PXRD) patterns were recorded using a Rigaku MiniFlex 600-C X-ray diffractometer with Cu Kα radiation. The samples were scanned with a step size of 0.01°, and their PXRD patterns were compared with reference databases using JADE software for estimation of mineral composition.

Radiocarbon dating ancient oysters

We prepared samples for ¹⁴C dating by scraping off the outer, weathered portion of the reef overgrown by microbial mats (Table 2). The shells of S. cuccullata within the sampled section were cleaned by Milli-Q water before they were sent to the Australian National University Radiocarbon Laboratory and Yale Analytical and Stable Isotope Centre for Accelerator Mass Spectrometry (AMS) radiocarbon dating. The six new dates, along with the previous date reported¹⁴, were calibrated using the IOSACal Python Library⁶⁹, which incorporates the latest Marine20 calibration dataset⁷⁰, assuming a local marine reservoir offset (∆R) of − 12 ± 71 years⁷¹.

Reconstruction of relative sea level

We constructed six new SLIPs, which constrain the unique position of past RSL over time, from the six new dated oyster samples and the dated sample reported in previous studies¹⁴. RSL was estimated following the equation:

where Elevation_i is the elevation of each ¹⁴C dated oyster sample i and RWL_i is its reference water level (RWL), the midpoint of the vertical range of S. cuccullata estimated from modern surveys, both expressed relative to PD (Table 1). The vertical uncertainty (indicative range uncertainty) is derived from half the 2σ range that S. cuccullata can occupy in the contemporary environment⁷².

Synthesis and modeling of RSL data

SLIPs from BWB were compared to other standardized RSL data from Hong Kong and the PRD region. We updated the existing PRD dataset^9,10,11 and incorporated dated samples from recent publications since 1996. For a complete reference list of compiled publications, see Supplementary information 2. Database construction followed international community standards⁷³ of the HOLSEA protocol (https://www.holsea.org; Supplementary information 2). For each sample, we critically assessed its: 1) geographic coordinates, 2) age calibrated to thousands of years before present (ka), 3) elevation, and 4) indicative meaning based on its stratigraphic context and diagnostic features (Table S4). We incorporated additional uncertainty from changes in paleotidal range⁷⁴. Rigorous quality screening was conducted to systematically evaluate the reliability of each data point. The screening criteria (Table S5) were applied based on: 1) the reliability of the sample’s location, age, and elevation data; and 2) the strength of the supporting evidence for determining each sample’s indicative meaning. L1 samples satisfied all criteria, whereas L2 samples could possibly be influenced by unresolved sources of error that are not fully reflected in the uncertainty estimates, and therefore should be interpreted with caution. We excluded L3 samples—for which the location, age, elevation, or indicative meaning could not be reasonably estimated—from further analysis.

To evaluate whether inferences from the SCS are supported by far-field locations globally, we synthesized Holocene RSL records from regions including Australia (North Queensland⁴⁶, New South Wales⁴⁷), Southeast Asia (Singapore⁴⁸), Africa (Langebaan Lagoon, South Africa⁴⁹), the Mediterranean (Tunisia⁵⁰), the Caribbean (St. Croix⁵¹) and South America (Suriname and Guyuana⁵², Rio Grande do Norte State, Brazil^53,54).

To quantitatively estimate temporal trends in Holocene RSL evolution with uncertainties, we applied the ET-GP²⁵ model to L1 data in the Hong Kong and PRD and to the global records following ref⁴. This model accounts for the vertical and chronological uncertainties of the RSL data to provide a probabilistic assessment of past RSL changes. Briefly, RSL over time is characterized by a zero-mean Gaussian process prior⁷⁵ represented by a Matérn covariance function with smoothness parameter ν = 3/2, defined by a prior amplitude and timescale of RSL variability. The proxy records provide noisy observations of RSL with geochronological uncertainty⁷⁶, assuming both age and RSL errors are normally distributed. The vertical error variance includes a supplemental white noise term to account for variations in the data that are not captured by the Matérn covariance function. Hyperparameters, which characterize timescales of variability in the Matérn covariance and the way in which information is shared between time points, are tuned to maximize the marginal likelihood of the model conditional upon the RSL data (Table S6). The 2σ uncertainties in RSL and rates are provided in “Results and discussion” section. To assess the sensitivity of the minimum highstand estimate to model specification, sensitivity tests were performed on the ratio of amplitude to timescale hyperparameters and the choice of Matérn covariance function (Fig. S6; Table S7).

Sensitivity test on GIA model predictions

We compared our RSL data with GIA predictions from two ice models with different ice-equivalent functions that developed independently—ICE-6G_C⁵⁵ and the ANU-ICE model⁷—to assess how well they capture early- to mid-Holocene barystatic signals in regional sea-level records. The ice-equivalent function of ICE-6G_C was constrained by near-field data (e.g., ice margin chronologies, GPS measurements, and RSL data) and was tuned to fit the Barbados tectonically-corrected coral RSL record, while ANU-ICE was tuned to fit RSL data from far-field locations including the coast of the SCS such as the PRD.

Both models were paired with one 3D Earth model, HetML140⁷⁷. The HetML140 model includes a laterally varying elastic lithosphere, a 40 km-thick uniform sub-lithosphere, two layers of 3D Earth structure in both the upper and lower mantle. Previous regional studies have shown that 1D viscosity models cannot simultaneously fit RSL data from both China and Malay-Thailand regions, which require different upper mantle viscosity values due to local tectonic influences⁷⁸. The incorporation of a 3D Earth structure allowed the model to account for lateral heterogeneity in mantle viscosity, and the model HetML140 was selected because it yields a superior fit to near-field RSL data compared to VM5a when paired with ICE-6G_C⁷⁷.

To better understand the factors contributing to data-model agreement, we conducted a sensitivity analysis to assess the influence of deglaciation timing and 3D Earth rheology on predictions of highstand timing and magnitude from GIA models. The ICE-6G_C ice model was paired with the one-dimensional (1D) VM5a Earth model⁷⁹ and adopted as a reference model in the analysis. The reference model was adjusted to 1) delay global deglaciation by 1 ka, 2) delay Antarctic deglaciation only by 1 ka, 3) employ 3D upper mantle viscosity only, and 4) employ 3D lower mantle viscosity only, allowing an evaluation of influence from these individual factors on GIA model predictions.

We calculated the misfit between the RSL estimated from the ET-GP, RSLd, and RSL predicted by GIA models, RSLm using the following misfit statistic⁸⁰:

where Φ(M) represents the data-model misfit for model M, J is the total number of time steps for which RSLd was calculated from the available SLIPs in a given study area, RSLd_j and RSLm(M)_j are the observed and predicted RSL for the jth time step, and σ_j is the 1σ uncertainty associated with the observed RSL from the ET-GP. The higher the Φ(M) value, the greater the misfit between the predicted RSL from the GIA models and the estimated RSL. A key distinction from χ² analyses typically used in data-model comparison is that instead of comparing GIA-modelled RSL predictions directly with SLIPs, instead we use probabilistic RSL estimates from the ET-GP. This reduces bias from anomalous data by incorporating all observations into a unified framework that accounts for uncertainties. We calculated data-weighted average misfit (Φ_t) for each ice and Earth model combination by weighting site-specific Φ(M) values by the number of data points available at each site, ensuring that sites with more data had a proportionally greater influence on the overall model performance assessment.