# Foraminiferal environmental DNA reveals late Holocene sea-level changes

Zhaojia Liu, Nicole S Khan, Howard K Y Yu, Arthur Chung, Magali

Schweizer, Celia Schunter

<figure><img src="image_1.png" /></figure>

► To cite this version:

Zhaojia Liu, Nicole S Khan, Howard K Y Yu, Arthur Chung, Magali Schweizer, et al. Foraminiferal environmental DNA reveals late Holocene sea-level changes. Communications Earth & Environment, 2025, 6 (1), pp.1029. 10.1038/s43247-025-03001-w. hal-05444752

HAL Id: hal-05444752

https://hal.science/hal-05444752v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

# Foraminiferal environmental DNA reveals late Holocene sea-level changes

<figure><img src="image_4.png" /></figure>

Check for updates

Zhaojia Liu©<fnref n="1" />, Nicole S. Khan©<fnref n="1" />✉, Howard K. Y. Yu©<fnref n="1" />, Arthur Chung©², Magali Schweizer©³ & Celia Schunter©²,⁴

Reconstructing past relative sea level provides critical insight into mechanisms driving sea-level change and informs future projections. Foraminifera are widely used sea-level proxies, but their application is often limited by poor preservation. Here, we demonstrate that foraminiferal environmental DNA and sedimentary ancient DNA provide a complementary approach to traditional morphological methods for relative sea-level reconstruction. By analyzing surface sediments and a core from subtropical intertidal environments in the Pearl River Delta, we found a clear vertical zonation in the environmental DNA assemblage consistent with morphological results. An environmental DNA-based transfer function enabled reconstruction with decadal temporal and decimeter vertical resolution for two periods: 290–1703 CE and 1956–present. Notably, sedimentary DNA preservation extended the reconstruction beyond morphological methods, which was limited by taphonomic processes. The environmental DNA reconstruction closely matched tide-gauge and geological records, underscoring its potential as a robust tool for reconstructing past relative sea level and its driving mechanisms.

Holocene relative sea level (RSL) records extend knowledge of mechanisms driving RSL change beyond the limited duration instrumental record<fnref n="1" />⁻⁴. These records enable comparison with other reconstructed climatic variables and model-based predictions to identify drivers of RSL change and their interactions within the climate system⁵⁻⁸ and thus inform future sea-level projections⁴,⁹.

Holocene RSL reconstructions are derived from geological proxies that exhibit systematic and quantifiable relationships to tidal elevation and thus RSL¹⁰. Among these, salt-marsh foraminifera are widely used due to their (1) ubiquity in low-energy intertidal environments across the globe; and (2) low-diversity assemblages found in high numbers¹¹,¹². Their sensitivity to subaerial exposure (inundation frequency) results in distinct vertical zonation within the tidal frame¹³⁻¹⁵. Knowledge of their present-day distribution and ecological preferences can be used to interpret assemblages of fossil foraminifera preserved in sediment archives to reconstruct past tidal elevation and produce high-resolution (centimeter to decimeter-scale) records of RSL change¹⁶⁻¹⁸. However, a major limitation of this approach is the poor or selective preservation of foraminiferal tests, particularly in (sub) tropical settings where elevated temperatures and precipitation promote test degradation or dissolution¹⁹⁻²¹. This problem in part arises from the traditional morphological methods used to identify and enumerate foraminiferal

assemblages, which rely on the identification of individuals under a ste-
reomicroscope based on test morphology.

Environmental DNA (eDNA), DNA obtained from environmental samples (soil, air, and water) without isolating or culturing target organisms, offers a promising alternative²². eDNA metabarcoding using a fragment of the small subunit of ribosomal RNA gene (SSU rDNA) is increasingly being used in foraminifera studies, including environmental biomonitoring²³⁻²⁶ and biodiversity assessments²⁴,²⁵,²⁷. Compared to traditional morphological methods, eDNA reveals greater diversity in modern benthic foraminifera, especially for organic-walled or finely-agglutinated monothalamids²⁸. DNA released after cell lysis can bind to environmental compounds, such as clay minerals, humic acids and other charged particles, providing protection from degradation by microbial DNases (deoxyribonucleases) in the environment²⁹⁻³². Additionally, the relatively anoxic and radiation-free environment in buried sediments facilitates the preservation of this sedimentary ancient DNA (sedaDNA) for hundreds to hundreds of thousands years³³⁻³⁵, even in tropical settings³⁶,³⁷. SedaDNA has therefore been used in foraminiferal paleoecological and paleoceanographic reconstructions³⁸⁻⁴³. However, these reconstructions have primarily focused on open-ocean or fjord environments, where foraminiferal assemblage patterns are influenced by environmental variables related to water depth, sea-surface temperature,

<footnote n="1">¹Department of Earth and Planetary Sciences and the Swire Institute of Marine Science, The University of Hong Kong, Hong Kong SAR, China. ²Swire Institute of Marine Science, School of Biological Sciences, The University of Hong Kong, Hong Kong SAR, China. ³Université d'Angers, Nantes Université, Le Mans Université, CNRS, Laboratoire de Planétologie et Géosciences, LPG UMR 6112, 49000 Angers, France. ⁴State Key Laboratory of Marine Environmental Health and Department of Chemistry, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong SAR, China. ✉e-mail: nskhan@hku.hk</footnote>
or distance to land. In contrast, low-energy intertidal environments—which produce precise RSL records—are characterized by benthic foraminiferal assemblages shaped by inundation frequency and salinity gradients¹³⁻¹⁵, and potential for eDNA preservation may be poorer under more oxic conditions and greater UV exposure¹⁴,⁴⁵. Despite the potential of foraminiferal eDNA to overcome limitations associated with test preservation and reveal greater diversity, foraminiferal eDNA remains unexplored as a sea-level proxy.

Here we evaluate the utility of foraminiferal eDNA and sedaDNA as a sea-level proxy in intertidal environments from subtropical Hong Kong. We characterize contemporary foraminifera vertical distribution using eDNA metabarcoding and morphological methods from three transects and examine the environmental factors controlling their distribution. To better reconcile taxonomy between eDNA and morphological methods, we develop an in-house DNA sequence database of dominant taxa identified through morphological analysis. We use these contemporary assemblages to

interpret a radiocarbon-dated mangrove core. A Bayesian transfer function (BTF) is developed from the modern eDNA training set, incorporating prior information from the elevation range of environmental zones identified by stable carbon isotope geochemistry, to reconstruct paleo-mangrove elevation and RSL change. Despite limitations with the preservation quality of sedaDNA in sediment layers lacking foraminiferal tests, this approach provides robust reconstructions of paleo-mangrove elevation (PME) and RSL change that are consistent with regional tide gauge and RSL datasets. These findings establish foraminiferal eDNA and sedaDNA as a valuable tool for paleoenvironmental and sea-level studies.

## Results

Samples were collected from a site where foraminifera, diatoms, geochemical proxies and chronological analyses were previously examined for paleoenvironmental and RSL reconstructions⁴⁶,⁴⁷ (Fig. 1). The site was

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Location of the study area.** The study area is located in the Pearl River Delta (PRD) (a) and in Deep Bay (b). Location of surface sampling transects (MP_A, MP_B, and SK), sediment core MPSC01, the nearest tide gauge at Tsim Bei Tsui, and the land surface elevation as represented by the LiDAR digital elevation model</figcaption></figure>

provided by the Hong Kong S.A.R government in Mai Po and Sha Kiu (c). Historical boundaries of mangrove in Mai Po Nature Reserve inferred from aerial images⁴⁷ are indicated by the black dashed line.
chosen due to its proximity to a tide gauge and availability of historical aerial
imagery showing mangrove and mudflat extent, enabling validation of
reconstructions. Previous investigations identified foraminifera through
morphological analyses conducted across three transects and sediment core
MPSC01<sup>47</sup>. Six environmental zones—subtidal, mudflat, mangrove fringe,
mangrove, terrestrial transition, and terrestrial—were differentiated based
on δ<sup>13</sup>C, total organic carbon (TOC), and carbon to nitrogen (C/N) ratios
geochemistry (Supplementary Fig. 1; Supplementary Data 1). Together with
radiocarbon dating of the core, these proxies revealed a transition from tidal
flat organic mud to mangrove peaty mud and muddy peat at 1957 CE
(Supplementary Data 2).

Foraminiferal eDNA/sedaDNA assemblages

Analysis of foraminiferal eDNA from modern samples and sedaDNA from core samples revealed substantial genetic diversity, encompassing 824 operational taxonomic units (OTUs) derived from 20,897,581 reads from 60 modern and 28 core samples after filtering out non-foraminiferal, low-quality, and chimeric reads. Taxonomic resolution varied, with 19 OTUs identified at the species level (8% of total sequences), 3 at the genus level (0.1%), 327 at the family level (48%), and 4 at the order level (0.2%) based on GenBank comparisons. An additional 9 OTUs were assigned to species level using an in-house database developed by sequencing abundant local taxa absent in public databases (Supplementary Data 3). The remaining 462 OTUs were identified as foraminifera (Phylum Foraminifera), but could not be assigned to higher taxonomic levels (Supplementary Data 4).

Of the assigned OTUs, 169 belonged to calcareous taxa, 97 to agglu-
tinated taxa, and 96 to monothalamids. The most abundant hard-shelled
families were Ammoniidae and Miliamminidae, collectively representing
33% of the total sequences, while other hard-shelled families (e.g. Tro-
chamminidae, Elphidiidae) contributed fewer (7%) sequences (Supple-
mentary Figs. 2, 3; Supplementary Data 4). Monothalamous OTUs (e.g.,
Allogromiidae and Saccamminidae) contributed 19% of total sequences.
Undetermined foraminiferal OTUs accounted for 41%, with two dominant
OTUs (OTU3 and OTU6) contributing 12% (Supplementary Figs. 2, 3;
Supplementary Data 4).

Distinct assemblages were observed in modern samples across the elevation gradient, expressed in standardized water level index units (SWLI) that normalize tidal elevations across sites with different tidal characteristics. At high elevations above 100 SWLI (equivalent to mean tide level [MTL]), Miliamminidae dominated (relative abundances of sequences ranging from 1–93%) along with undetermined OTUs (e.g. OTU3: ≤43%; OTU11: ≤17%), Ammoniidae (≤42%), and Trochamminidae (≤14%) (Fig. 2). Ammoniidae were prevalent (1–55%) in upper mudflat and mangrove environments (57–200 SWLI), while Nummulitidae (1–23%) and OTU6 (≤73%) dominated low-elevation mudflat and subtidal environments (<57 SWLI). Monothalamous families exhibited broad elevation tolerance, occurring across subtidal to mangrove environments (Allogromiidae: ≤15%; Saccamminidae: ≤62%) (Fig. 2c).

Morphological foraminiferal assemblages

A total of 8660 individuals representing 24 morphospecies were identified from modern and core samples, including 14 agglutinated and 10 calcareous taxa across 12 families (Supplementary Figs. 2, 3). Seven families from the morphological dataset linked to 191 OTUs that contribute 35% of the total reads (33% in modern samples and 44% in core samples). Five morphospecies from five different families—*Acostata* sp., *Bruneica* sp., and *Glo-mospirella* spp. (lacking reference sequences in GenBank), as well as *Haplophragmoides wilberti* and *Rosalina* sp.—were not detected in the eDNA/sedaDNA assemblage (Supplementary Fig. 2, 3).

Similar to eDNA and sedaDNA dataset, Miliamminidae and
Ammoniidae were dominant, with Miliamminidae present in 92% of
modern samples (Fig. 2a), contributing 41% and 53% of total abundance in
modern and core samples, respectively (Supplementary Figs. 2, 3). However,
Miliamminidae sequences were rare below MTL. Ammoniidae exhibited
similar distribution patterns in mangrove environments across both

datasets (Fig. 2b). Some taxa abundant in the morphological assemblages (e.g., *Ammonia sobrina*, absent from GenBank, and *A. tepida* and *A. falso-beccarii*) were not found in eDNA data, causing discrepancies in low-elevation environments (Fig. 2b).

Low-elevation subtidal and mudflat environments (<57 SWLI) were dominated by calcareous taxa, including Ammoniidae (e.g., *A. falso-beccarii*: ≤76% in each sample; *A. sobrina*: ≤53%), Rosalinidae (≤43%), and agglutinated Lituolidae (≤81%), collectively accounting for 68% of counts in these environments (Fig. 2a, b). In upper mudflat and mangrove fringe environments (60–136 SWLI), Lituolidae (≤56%), Ammoniidae (*Ammonia* spp.; 5–32%; *A. tepida*: ≤34%) dominated, together contributing 71% of total counts (Fig. 2a, b). At higher elevations (>136 SWLI), agglutinated Trochamminidae (*A. mexicana*: ≤41%; *T. inflata*: ≤29%) and Haplophragmoididae (≤40%) became increasingly abundant (Fig. 2a).

Environmental variables shaping morphological and eDNA
assemblages

Canonical correspondence analysis (CCA) and partial CCA (pCCA) were
performed to assess the environmental controls on eDNA and morpholo-
gical assemblages, revealing salinity, δ¹³C, and elevation as key environ-
mental factors shaping eDNA assemblages, which explained 20%, 16%, and
15% of the variance, respectively (Fig. 3a, b; Supplementary Data 5, 6). In
comparison, elevation (23%), followed by sand content (17%) and salinity
(13%) explained the variability in morphological assemblages (Supple-
mentary Data 6). Non-metric multidimensional scaling (NMDS) using
Bray-Curtis dissimilarities of community composition revealed elevation
preferences of taxa from both morphological and eDNA assemblages
(Fig. 3c, d). The community composition was significantly influenced by
elevation, with strong effects observed for both morphological (PERMA-
NOVA; F = 33.43, R² = 0.366, *p* = 0.001) and eDNA assemblages (F = 27.49,
R² = 0.325, *p* = 0.001).

Downcore signatures of foraminiferal sedaDNA and morpholo-
gical assemblages

Three lithostratigraphic units (Unit I-III) were identified in core MPSC01 (Fig. 4). Foraminiferal sedaDNA was successfully extracted and amplified from sediments above 0.85 m Hong Kong Principal Datum (PD), indicating sufficient preservation to permit paleoenvironmental interpretation. Nevertheless, the observed OTUs in the sedaDNA assemblage was significantly lower than those in the modern eDNA assemblage (ANOVA; *p* < 0.05). Units II and III displayed similar α-diversity metrics, as measured by observed OTUs, Shannon, and Chao1 indices (Supplementary Data 9). In contrast, the sedaDNA assemblage in Unit I exhibited significantly lower observed OTUs and Chao1 indices (ANOVA; *p* < 0.05), suggesting the potential influence of taphonomic processes with depth.

Foraminiferal sedaDNA reveals distinct changes in core assemblages with depth. The sedaDNA assemblage in the basal lithostratigraphic unit (Unit I; 0.85–1.44 m PD)—deposited between 486 (290–652, 2σ uncertainty range) CE and 1957 (1957–1958) CE—resembled the modern mudflat eDNA assemblage, characterized by low Miliamminidae abundance (≤17%), moderate Ammoniidae (≤30%), and substantial Saccamminidae (6–69%) contributions (Fig. 4). These patterns align with linear discriminant analysis (LDA) applied to core geochemistry, which also indicated a mudflat environment (Fig. 4).

In Unit II (1.45 to 1.89 m PD), which was deposited between 1957 and 1991 CE, sedaDNA assemblages were dominated by Miliamminidae (14–82%), with moderate contributions from Ammoniidae (2–16%) and Saccamminidae (2–18%) (Fig. 4). This shift, characterized by increased Miliamminidae sequences, reflects a transition from a low-elevation mudflat to a mangrove environment and is consistent with LDA predictions (Fig. 4). An exception at 1.85–1.86 m PD showed Saccamminidae dominance (91%), which lacked modern analogues and likely indicates a localized disturbance or depositional change at this boundary between lithostratigraphic units.
<figure><img src="image_6.png" /><figcaption>**Fig. 2** | Vertical distribution of major taxa in the morphological and eDNA assemblages. Panels show agglutinated taxa (a), calcareous taxa (b), and taxa in the eDNA assemblage without counterparts in the morphological assemblage (c). Both morphological and eDNA assemblages are shown in relative abundance (%).</figcaption></figure>

Elevation is presented as Standardized Water Level Index units (SWLI), of which a value of 100 represents mean tide level (MTL) and a value of 200 represents mean higher high water (MHHW).

Similarly, the sedaDNA assemblages in Unit III (1.89 to 2.15 m PD), deposited after 1991 CE, were consistently dominated by Miliamminidae (46–81%), a pattern characteristic of modern mangrove environments and supported by LDA results (Fig. 4). However, in samples near the top of the

core (above 2.15 m PD), Miliamminidae abundance dropped (<8%), while Saccamminidae (≤38%) and undetermined OTUs (e.g., OTU25 ≤ 58%) increased. This assemblage lacks a modern analogue in the training set and may reflect influence by bioturbation.
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Multivariate ordination analyses conducted on morphological and DNA assemblages reveal the influence of examined environmental variables on taxa composition. Taxa-environment Canonical Correspondence Analysis (CCA) biplots of morphological (a) and DNA (b) taxa that exceed 5% of the assemblage in at least one sample. The lengths of the environmental arrows on the taxa-environment biplot approximate their relative importance of influence on the foraminiferal assemblage. Arrow direction reflects its approximate correlation with the ordination axes and other environmental factors, and the perpendicular projection of taxa position on the environmental arrow indicate their weighted average optima along each environmental variable. The first two axes of eDNA and morphological CCA</figcaption></figure>

explained 16% and 37% of total variance, respectively. Non-metric Multi-Dimensional Scaling (NMDS) ordination of morphological (c) and DNA (d) assemblages of all modern surface sampling stations and core samples based on Bray-Curtis similarities. Samples collected from MP_A, MP_B and SK transect were marked with different colors. Core samples are shown as depth below surface (cm). Grey lines indicate the relative standardized water level index units (SWLI) elevations associated with sample positions in NMDS ordination space. Some major undermined OTUs and assigned taxa of eDNA assemblage are shown. Taxa names are written in abbreviation (Supplementary Data 7, 8).

Compared to the sedaDNA assemblage, the morphological assemblage preserved in the core exhibited a consistent composition across all lithos- tratigraphic units (1.30–2.31 m PD), being dominated by agglutinated taxa, including Miliamminidae (8–76%) and Trochamminidae (A. mexicana, 10–39%; *T. inflata* ≤ 20%) (Fig. 4). This composition reflects a mangrove environment, aligning with interpretations from the sedaDNA assemblage. However, although no significant difference in α-diversity metrics was observed between modern and core morphological assemblages (ANOVA; *p* > 0.05), taphonomic issues in the core were more evident than with the sedaDNA assemblage. Foraminiferal tests were entirely absent below 1.30 m PD in Unit I (deposited before 1956 CE), and tests were sparse (<100 individuals per 1 cm³ sediment) between 1.55 and 1.61 m PD (Fig. 4), preventing paleoenvironmental interpretations of sedimentary archive based solely on morphological assemblages.

## Performance of eDNA and morphological training sets

We developed Bayesian transfer functions (BTFs) from eDNA and morphological training sets to estimate response curves (mean with 95% credible interval) for each taxon in relation to tidal elevation (SWLI). Many taxa in the eDNA assemblage exhibited unimodal distributions with tidal elevation, including several abundant OTUs (e.g. OTU26 and OTU29) and monothalamous species (e.g. *Vellaria pellucida*). However, some taxa in both assemblages displayed bimodal or multimodal distributions, including abundant families Saccamminidae and Ammoniidae (Fig. 5b, c).

To evaluate BTF performance, a 10-fold cross-validation was conducted on both eDNA and morphological datasets (Fig. 5; Table 1; Supplementary Data 10). Both BTFs showed strong correlations between observed and predicted elevation (R² = 0.95 for eDNA, 0.93 for morphology). After removing outliers with cross-validation residuals exceeding two
<figure><img src="image_8.png" /><figcaption>**Fig. 4** | **Core MPSC01.** Presented are the core lithology, chronology established by Bacon age-depth model (age outliers were removed but shown in Supplementary Fig. 4), geochemistry variables and associated linear discriminant analysis (LDA) prediction of environmental zones, paleo-mangrove elevation (PME) prediction inferred from LDA, foraminiferal morphological and sedaDNA assemblages, reconstructed PME estimated by foraminiferal eDNA- and morphology-Bayesian transfer function (BTF), dissimilarity index based on Bray-Curtis (BC) distance, and goodness-of-fit statistics. Morpho- and eDNA taxa are grouped at the family level for</figcaption></figure>

comparability. “Others” includes OTUs assigned to order level. The envelop of BTF predictions represent $\pm 2\sigma$ uncertainties. For the dissimilarity index, the 10th percentile of the sedaDNA datasets and the 20th percentile of the morphological datasets was normalized to 1, with values > 1 suggesting poor analogy to a modern counterpart. For goodness-of-fit statistics in core samples, the 95% threshold of squared residual lengths was normalized to 1. Values below 1 indicate a poor fit of the transfer function to elevation.

standard deviations from the mean, the morphology-BTF slightly outperformed the eDNA-BTF in accuracy (98% vs. 94%), but the eDNA-BTF had lower uncertainty ($1\sigma$ = 18 vs. 34 SWLI; 0.17 vs. 0.32 m) and RMSEP (25 vs. 32 SWLI; 0.24 vs. 0.3 m). Elevated uncertainty and outliers were mainly found in the low and high ends of the SWLI range for the eDNA-BTF, and the morphology-BTF tended to overpredict elevation in lower-mudflat environments (Fig. 5a).

## Reconstructed paleoenvironmental and paleo-mangrove elevation change

The eDNA-BTF predicted an increase in paleo-mangrove elevation (PME) from 131 to 210 SWLI from core base (0.85 m PD) to top, with an average $1\sigma$ uncertainty of 12 SWLI (equivalent of 0.11 m). In comparison, the morphology-BTF estimated a rise in PME from 179 to 204 SWLI from 1.3 m PD to core top, with an average $1\sigma$ uncertainty of 18 SWLI (0.17 m) (Fig. 4).

To evaluate the reliability of our reconstruction, we compared it with annual mean sea level (MSL) changes in Hong Kong recorded by tide gauges since 1963 CE and other sea-level data from the Pearl River Delta (PRD) (Fig. 6). Excluding core samples exceeding dissimilarity index thresholds, the eDNA-BTF estimated a RSL rise of -0.01 to 0.08 meters between 1978 and 2016 CE (Fig. 6a), capturing all tide gauge records within $2\sigma$ uncertainties. This trend was consistent with the morphology-BTF, but with a slightly higher mean squared error (MSE: eDNA = $0.027 \text{ m}^2$; morphology-BTF = $0.016 \text{ m}^2$) (Fig. 6c). The eDNA-BTF extended the RSL record back to 486 (290–652) CE, supported by index and limiting points from the PRD (Fig. 6d).

## Discussion

The reconstruction of RSL using transfer functions assumes that tidal elevation, a surrogate variable for the frequency and duration of tidal inundation, is an important ecological determinant of taxa distribution, or linearly related to another ecologically significant factor$^{48,49}$. Indeed, elevation consistently influenced both eDNA and morphological assemblages, with samples from all three transects clustering by elevation in NMDS space (Fig. 3c, d). However, additional environmental factors, particularly salinity, also shaped eDNA assemblages. Salinity, which decreases with increasing elevation, was identified as a primary control on foraminiferal eDNA composition, in agreement with previous studies of estuarine foraminiferal

eDNA$^{25}$. Certain undetermined OTUs and monothalamids, for example OTU29, were strongly associated with high-salinity environments (Fig. 3b), which helps explain why elevation accounted for less variation in eDNA assemblages compared to morphological data. In contrast, tidal elevation remained the dominant variable explaining morphological assemblages (Supplementary Data 6), consistent with previous research$^{50–52}$.

The divergence in dominant environmental controls likely originates from differences in the populations represented by the eDNA and morphological assemblages. Morphological assemblages used to develop transfer functions are typically derived from tests of dead benthic adult foraminifera to minimize temporal bias$^{53}$. In contrast, eDNA inevitably integrates DNA from multiple sources, including living adults, juveniles, propagules, allochthonous DNA, and DNA preserved in dead foraminiferal tests$^{25,54,55}$. For example, family Nummulitidae, despite being mineral-shelled, was identified only in the eDNA assemblage, suggesting the presence of propagules rather than adult individuals. The mobility of eDNA, similar to diatom assemblages in intertidal zones$^{56–58}$, allows it to capture signals from both autochthonous and allochthonous taxa (e.g., planktonic taxa Globigerinitidae; Supplementary Data 4). In Mai Po and Sha Kiu mudflat samples, the low proportion of stained (live) foraminiferal tests (Supplementary Data 11) suggests limited in situ living populations, implying a greater contribution of propagules and allochthonous DNA to the eDNA assemblage. The strong influence of $\delta^{13}\text{C}$-indicative of a marine or terrestrial source of sediment-on the eDNA assemblage further underscores its multi-source nature. High $\delta^{13}\text{C}$ values suggest greater contributions from marine organic matter sources$^{46}$ and hence stronger contributions of allochthonous eDNA in these lower-elevation environments more open to external influences$^{59}$. Although dissolution of calcareous tests under low pH conditions$^{19,20,47,60}$ may contribute to differences between morphological and eDNA assemblages, it is unlikely to account for the divergence in dominant environmental controls because pH explained only a small proportion of the variance in both eDNA (3%) and morphological assemblages (4.7%) (Supplementary Data 6) and calcareous taxa were abundant across both mangrove and mudflat zones, indicating minimal taphonomic bias.

Another factor contributing to the observed discrepancies is the inclusion of monothalamous taxa, often overlooked in morphological analyses due to their fragile tests. Monothalamids, which are abundant in eDNA assemblages from coastal and deep-sea settings$^{25,26,61}$, may be better
<figure><img src="image_9.png" /><figcaption>**Fig. 5** | Ten-fold cross validation and taxa response curves modeled by the Bayesian transfer function (BTF). BTFs predicted elevation by 10-fold cross-validation versus observed elevation measured at sampling stations of foraminiferal eDNA and morphological assemblages modern training set are shown in (a). Elevation is expressed as standardized water level index units (SWLI). Residuals are calculated by subtracting predicted value from observed elevation. For both</figcaption></figure>

assemblages' BTFs, samples that were screened out due to residuals exceeding two standard deviations from the mean residual in cross-validation are also shown, but highlighted in red. Taxa response curves of major taxa from eDNA (b) and morphological (c) assemblages. Circles represent the relative abundance of corresponding taxa at observed sample elevations, the red solid line and the grey shaded envelope represent the model mean and the 95% credible interval of prediction.

represented in eDNA datasets because DNA from their soft tests is more easily extracted⁶² and the short DNA barcodes of some taxa (e.g., saccaminids) improve PCR amplification and sequencing efficiency⁶³. As shown in a previous eDNA-based study²⁵, monothalamous taxa are sensitive to short-term physiochemical environmental changes. Their inclusion in foraminiferal eDNA assemblages likely amplifies the influence of environmental variables like salinity. Nevertheless, tidal elevation is a key ecological driver of foraminiferal zonation and exhibits strong correlations with other significant environmental variables that influence eDNA assemblage composition. Meeting this key assumption permitted the development of BTFs.

The lower uncertainty and bias of the eDNA-BTF compared to the morphology-BTF (Table 1) underscore its potential for precise and accurate elevation predictions in intertidal environments. While the morphology-BTF exhibited systematic overprediction in mangrove environments (170–200 SWLI) (Fig. 5a), likely due to the linear increase in *Miliammina* spp. abundance in higher elevations (Fig. 5c), the eDNA-BTF showed no discernible pattern in the residuals. However, relatively greater uncertainties in lower mudflat environments (<MTL) (Fig. 5a) suggests that highly-mobile propagules and allochthonous eDNA²⁵ may weaken the mudflat

<table><caption>**Table 1** | Summary of eDNA- and morphology-Bayesian transfer functions (BTFs) performance in 10-fold cross-validation of the modern training set</caption><thead><tr><th>Metric</th><th></th><th>eDNA-BTF</th><th>Morphology-BTF</th></tr></thead><tbody><tr><td>Correlation (R<sup>2</sup>)</td><td></td><td>0.95</td><td>0.93</td></tr><tr><td rowspan="2">1σ Uncertainty</td><td>SWLI</td><td>18</td><td>34</td></tr><tr><td>Meter</td><td>0.17</td><td>0.32</td></tr><tr><td rowspan="2">RMSEP</td><td>SWLI</td><td>25</td><td>32</td></tr><tr><td>Meter</td><td>0.24</td><td>0.3</td></tr><tr><td rowspan="2">Average Residual</td><td>SWLI</td><td>0.7</td><td>-9</td></tr><tr><td>Meter</td><td>0.01</td><td>-0.08</td></tr><tr><td>Accuracy (samples within 2σ uncertainty range of average residual)</td><td></td><td>94%</td><td>98%</td></tr></tbody></table>
<figure><img src="image_10.png" /></figure>

**Fig. 6** | Comparison of reconstructed relative sea level (RSL) with tidal gauge records and regional RSL records from the Pearl River Delta (PRD). Comparison of RSL reconstructions from foraminiferal eDNA- (a) and morphology- (b) Bayesian transfer functions (BTFs) with tide gauge data from Hong Kong. Boxes represent ±2σ vertical uncertainties and ±2σ age uncertainties for the reconstructions. RSL reconstructed from core samples with poor modern analogue were shown with transparent colors. RSL estimated by the core sample that has poor fit with tidal

elevation determined by goodness-of-fit were shown by red box. Comparison of eDNA- and morphology-BTF reconstructed RSL with tidal gauge records is shown in (c). Comparison of RSL estimated by eDNA-BTF with Hong Kong tide gauge records and regional index points—delimiting the position of RSL—and limiting points—providing upper (terrestrial) or lower (marine) bounds on RSL—from the PRD is shown in (d).

elevation signal. This pattern aligns with findings from diatom-based studies, where low-elevation assemblages are often influenced by planktonic taxa and allochthonous sources⁵⁷,⁶⁴⁻⁶⁶. Consistent with morphological studies⁴⁷,⁶⁷,⁶⁸, the clearest elevation signals—and thus the most reliable predictions—are observed in higher-elevation mangrove assemblages, where tidal elevation exerts a more dominant ecological control.

In addition to the correlation of modern foraminiferal eDNA with tidal elevation and its robust performance in BTF, the preservation potential of sedaDNA underpins the application of foraminiferal eDNA for RSL reconstruction. The sedaDNA assemblage showed generally good agreement with depositional environments inferred from LDA predictions across lithostratigraphic units (Fig. 4). Basal Unit I contained mudflat assemblages with low Miliamminidae abundance and mudflat-specific taxa (e.g., undetermined OTU29), while Units II and III showed mangrove assemblages indicated by high abundance of Miliamminidae (Fig. 4). However, the core sedaDNA assemblage had more samples (16/28) falling outside the similarity threshold than the morphological assemblage (3/20; Fig. 4). For the morphological assemblage, poor analogues mainly resulted from high abundances of *Ammoastuta inepta* and *Ammonia convexa*—rare in modern environments, likely due to anthropogenic impacts⁴⁷. Mangrove management policies at Mai Po (e.g., trimming of seedlings)⁴⁷ may have altered modern mudflat assemblages, creating poor analogues for both morphological and sedaDNA assemblages, especially pre-1965 CE core samples (Fig. 4). In addition to human influences, poor sedaDNA analogues may reflect: (1) high taxonomic diversity, (2) eDNA's sensitivity to subtle assemblage changes, (3) potential influence from infaunal populations, or (4) taphonomy and molecular preservation of sedaDNA.

First, the high taxonomic diversity of eDNA/sedaDNA assemblages increases heterogeneity and the proportion of rare taxa compared to morphological assemblages⁵⁹,⁷⁰, raising the likelihood of encountering taxa underrepresented in modern training sets⁴⁷, similar to other sea-level

proxies like diatoms⁶⁹,⁷¹. For example, OTU87, primarily found in mudflats, reached 14% abundance in a core sample (1.15–1.16 m PD, Unit I), but never exceeded 5% in modern samples (Supplementary Data 4).

Second, an abrupt increase of Saccamminidae at 1.85–1.86 m PD, coinciding with the lithostratigraphic boundary between Unit II and III (Fig. 4), also resulted in a poor modern analogue. This shift may reflect changes in hydrodynamics, productivity, or other environmental factors during facies transitions⁷²,⁷³, similar to the pattern observed at the mangrove-terrestrial transition in transect MP_B (Fig. 2). Monothalamous foraminifera like Saccamminidae are highly adaptive to fluctuating conditions and can respond rapidly to environmental changes, including seasonal cycles⁷⁴. Temporal variation in living populations may further contribute to discrepancies with modern analogues, highlighting the need for research on the seasonal dynamics of foraminiferal eDNA⁷⁵⁻⁷⁷.

Third, eDNA assemblages in the bioturbated mixed zone (>2.15 m PD) may have been influenced by living infaunal populations such as OTU11 and Elphidiidae, which lack preserved morphological counterparts either because they do not possess hard shells or are represented by propagule DNA⁴². Similar patterns have been reported in deep-sea sediments³⁸,⁷⁸, where the more abundant and pristine DNA from living foraminifera dominated the bioturbated zone⁴¹. Because DNA from living cells decays rapidly after cell death—primarily through endogenous nuclease activity (i.e., breakdown of nucleic acids originating from within the cell) rather than microbial activities⁴⁴,⁷⁹,⁸⁰—eDNA below the bioturbated zone primarily represents extracellular sources and the total population flux³⁸.

Lastly, both sedaDNA and morphological assemblages were affected by post-depositional processes. Below 1.31 m PD, foraminiferal tests were absent—likely due to dissolution under low pH conditions¹⁹,⁶⁰—while sedaDNA remained preserved (Fig. 4), suggesting greater long-term stability for DNA than foraminiferal tests in such environments. The sedaDNA adsorption onto clay minerals shields genetic material from chemical and
biological degradation³¹,³²,⁸¹, whereas foraminiferal tests are more susceptible to dissolution and disintegration. However, a statistically significant decrease in α-diversity was observed in core samples with increasing depth (Supplementary Data 9), indicating that the sedaDNA assemblage was also affected by taphonomic processes. Lower TOC (< 1.3%) in Unit I compared to the overlying units (Fig. 4) and lack of tests suggest reduced sedimentary shelter and protection for sedaDNA³¹,³⁸,⁴⁴. Additionally, the lower TOC—similar to the modern surface samples from mudflat and subtidal environments (Supplementary Fig. 1; Supplementary Data 1)—may reflect more oxygenated depositional conditions that promote intensified microbial activity⁸²,⁸³, thereby accelerating eDNA degradation during early deposition⁴⁵,⁸⁴. Monothalamous taxa with short barcodes (e.g., Saccamnidae) are more resilient to DNA fragmentation⁶³, which may increase their apparent abundance in degraded samples from Unit I and potentially lead to overestimation of PME, resulting in offsets between BTF and LDA estimates in core samples below 1.21 m PD (Fig. 4) and underprediction of RSL (Fig. 6d). While post-mortem damage and test loss preclude morphological analysis, eDNA metabarcoding can recover information from fragmented sequences as long as amplifiable DNA persists³³,³⁸,⁸⁵, making the eDNA method more applicable under extreme conditions (e.g. low pH) where morphological approaches are limited¹⁹,⁶⁰,⁸⁶⁻⁸⁸.

Foraminiferal sedaDNA and tests were absent, however, from a short interval between 1.25–1.31 m PD, likely resulting from a hiatus in sedimentation. This interpretation is supported by shifts in radiocarbon-determined ages and sedimentation rates at 1.31 m PD (Fig. 4). This hiatus, occurring around 800 years ago, is attributed to a lack of accommodation space⁸⁹ caused by late Holocene RSL fall¹⁶,⁹⁰, during which sediment accumulation ceased and sedimentary eDNA was exposed to bioturbation and weathering (prolonged warm, oxidative, and moist conditions in these subtropical wetlands) for at least 250 years, leading to the complete degradation or severe fragmentation of DNA⁴¹,⁴⁴,⁹¹⁻⁹³. The subsequent post-industrial RSL rise created accommodation for sediment deposition⁹⁴,⁹⁵, resulting in the rapid sedimentation rate (average >1.5 cm/yr) above 1.3 m PD (Supplementary Fig. 4).

Importantly, the presence of core samples with poor modern analogues does not necessarily compromise reconstruction reliability, particularly for heterogeneous datasets with high diversity⁴⁸. Applying the eDNA-BTF to reconstruct RSL from the core showed strong agreement between the tide gauge and other geologic RSL records from the PRD (Fig. 6), even for samples lacking close modern analogs (Fig. 4), supporting the approach's validity⁹⁶,⁹⁷. Our eDNA-based approach generated high-resolution RSL reconstructions with decadal temporal and decimeter vertical precision for two distinct chronological intervals (2σ range: 290–1703 CE and 1956 CE–present) and improved resolution for the period 290–1703 CE compared to existing geologic RSL records (Fig. 6d). Superior eDNA preservation enabled RSL reconstruction as far back as 290 CE, whereas morphological analyses could not be applied to sediments older than 1956 CE. The eDNA-BTF reconstruction demonstrated slightly lower accuracy in RSL estimation (MSE of 0.027 m² excluding samples with poor analogues) compared to the morphology-BTF (0.016 m²) and a morphology-BTF established from a salt-marsh core in New Jersey (0.014 m²)⁹⁸ but higher precision (1σ uncertainty: 12 vs. 18 SWLI). These results demonstrate that foraminiferal eDNA is a robust proxy for reconstructing late Holocene RSL, especially where morphological methods are not feasible.

Some limitations observed, such as the impact of sedimentary hiatuses, reflect broader challenges in developing high-resolution far-field RSL records⁴⁷,⁵²,⁸⁷ rather than inherent shortcomings of the eDNA proxy itself. The eDNA approach may have greater potential at sites characterized by rising late Holocene RSL, including mangrove archives from the Caribbean⁹⁹⁻¹⁰¹ and mid-Pacific¹⁰², or salt-marsh archives from across the north Atlantic⁶⁹,⁹⁶,¹⁰³⁻¹⁰⁵, where more continuous sedimentary records are preserved.

The comparison between eDNA and morphological assemblages reveals critical differences influencing their utility as proxies for reconstructing past RSL. The eDNA assemblages incorporate DNA from

propagules and juvenile stages, providing a broader genetic signal than morphological methods and offering valuable ecological insights. However, this broader signal introduces complexity when interpreting correlations with tidal elevation, as the detected DNA may not originate from in situ communities. To disentangle the influence of propagules, future studies could employ methods such as 10-g extraction kits²⁴ or focus on the >63 µm sediment fraction by sieving sediment prior to extraction to minimize propagule contributions to the measured intracellular eDNA pool¹⁰⁶,¹⁰⁷.

Although we found the eDNA method provides more precise tidal elevation signals than traditional morphological approaches, its application requires stringent field and laboratory procedures to preserve sample integrity and minimize cross-contamination²²,¹⁰⁸,¹⁰⁹. These challenges can be particularly acute when conducting remote fieldwork. Practical alternatives include sterile disposable tools, portable freezers, or well-insulated coolers with ice-salt mixes to maintain cold temperatures for extended periods. While eDNA sample preparation is more expensive than traditional morphological analyses, once optimized, it requires no more expertise or time than morphological species identification. Additionally, collaboration with commercial laboratories can mitigate costs and procedural complexity while ensuring reliable sequencing. These adaptations make eDNA methods more feasible for a wider range of research contexts, ensuring their broader applicability for reconstructing past RSL.

A notable property of eDNA methods is the inclusion of undetermined OTUs within assemblages. In this study, 53% of OTUs and 37% of sequences remained undetermined, a proportion consistent with previous foraminiferal eDNA research²³,²⁶,¹¹⁰. These undetermined OTUs likely represent species whose genetic information is not yet documented in reference databases, constraining the ability to interpret their ecological preferences and hindering cross-study comparisons. Without comprehensive taxonomic identification, it is challenging to relate OTUs across molecular and morphological studies, limiting the broader application and interpretability of eDNA-based reconstructions and their integration with morphological approaches. Although these undetermined OTUs have demonstrated the ability to indicate tidal elevation, expanding sequencing efforts for foraminifera—particularly from low-energy intertidal environments—remains essential for enabling the universal application of eDNA as a robust paleoenvironmental tool.

Finally, eDNA assemblages exhibit significantly higher diversity than morphological assemblages, with 824 OTUs identified, representing both species- and subspecies-level diversity, compared to 24 morphospecies (Supplementary Data 4). An additional key advantage of eDNA is its ability to detect taxa that are often underrepresented in traditional studies⁴⁰,⁴¹, such as monothalamids—the inclusion of which has been shown to enhance BTF performance (Supplementary Data 12). This ability to incorporate and detect genetic diversity at such high resolution could enhance the predictive power of eDNA in paleoenvironmental applications, particularly if this diversity reflects ecological sensitivity to environmental variables.

# Methodology

## Site selection

The study area is situated in southeast Deep Bay on the eastern Pearl River Estuary (Fig. 1b). This subtropical region experiences an average maximum temperature of ~30 °C in July and minimum of 17 °C in January. Average annual precipitation totals 2,396 mm (2019–2023), and most rainfall is recorded from June to September¹¹¹. Brackish salinities (9.7–29.7) are observed in Deep Bay¹¹² due to the discharge of the Pearl River and nearby Shenzhen River, and its average water depth is ~2.9 m¹¹³. Tides display a semi diurnal pattern in this area, with an amplitude of 1.96 m (mean higher high water [MHHW]-mean lower low water [MLLW]) recorded by the nearest tide gauge station at Tsim Bei Tsui (Fig. 1c).

Two study sites, Mai Po and Sha Kiu, were chosen to conduct sampling of modern sediment samples. The Mai Po Nature Reserve, which encompasses a large area of intertidal mudflat and mangrove, has been under management of the World Wide Fund for Nature (WWF) since 1975. Mangrove flora in Mai Po are dominated by *Kandelia candel*, *Aegiceras*
*corniculatum*, and *Avicennia marina*¹¹⁴. We chose Mai Po to conduct modern sample collection because the long-term management of this site ensures minimal direct human influence. The sediment core (MPSC01) was collected from the upper-mangrove of Mai Po (Fig. 1c) due to the long-standing presence of mangrove forest since at least the mid-20th century¹⁶. Sha Kiu, an area encompassing mudflat and mangrove environments, with flora dominated by *Aegicera corniculatum* and *Sonneratia apetala*, was also selected for modern sample collection due to its distinctive transition from mudflat to mangrove habitats, which we considered experienced relatively less human activity compared to other areas along the Deep Bay coast.

## Sample collection

To investigate the vertical zonation of foraminiferal environmental DNA (eDNA) assemblages and to construct the modern training set for relative sea-level (RSL) reconstruction, a total of 60 surface samples were collected along three transects: two in Mai Po (MP_A and MP_B) and one in Sha Kiu (SK), covering the subtidal, mudflat, mangrove and terrestrial forest areas (Fig. 1c). To evenly sample the elevation and environmental gradient, sampling stations were spaced at 5-cm elevation intervals, and in the flat mangrove forest of MP_A, every 30 m in distance. At each sampling station, surface (1-cm deep) sediment samples were collected individually for analysis of eDNA, morphology, and other environmental variables from a 10-cm² area. Individual disposable spatulas were used to place samples for eDNA analysis in 50 mL Falcon centrifuge tubes. Sediment samples from the subtidal zone were collected using a Van Veen grab sampler¹⁵, from which ~2 cm³ of the top 1 cm of surface sediments were sampled using spatulas and stored in 15 mL Falcon tubes. To minimize cross-contamination, sample collection was conducted with disposable sterile gloves, and all non-disposable equipment used in sample collection (e.g., spatulas, Van Veen grab sampler) were sterilized with 10% bleach prior to use and between each sampling station. All eDNA samples were kept on ice in the field prior to storage in a -20 °C freezer back in the laboratory on the sampling day¹⁶. Samples for traditional morphological analysis were preserved in an ethanol-calcite buffer solution and stained using rose Bengal in the field to differentiate between living and dead foraminifera¹⁷. Both samples for morphology and analysis of environmental variables were also kept in a cool box on ice in the field and stored in a 4 °C refrigerator in the laboratory.

Two 290-cm long replicate sediment cores (MPSC01) were collected from the Mai Po mangrove using an Eijkelkamp peat sampler in overlapping 50-cm intervals to minimize influence of compaction and contamination (Fig. 1c). Cores were sealed with plastic wrap in PVC pipes and kept on ice before returning to the laboratory. For the collection of the first replicate core used for sedimentary ancient DNA (sedaDNA) analysis, all materials and equipment (e.g., Eijkelkamp peat sampler, PVC pipes) were sterilized with 10% bleach before each 50-cm interval was sampled to prevent cross-contamination. Subsampling—conducted in a sterilized laboratory—was performed every 5 cm downcore using a disposable spatula. ~1 cm³ of sediment was sampled from the inner part of the core to further minimize contamination¹¹⁸ (because the sample's exterior surfaces came into contact with the peat sampler and PVC pipe during coring and storage) and stored in 15 ml Falcon tubes in a -20 °C freezer. The second core replicate (collected adjacent to the first) was subsampled for morphological, geochemical and chronological analyses. For morphological analysis, samples from a 1-cm thick interval from the second replicate core were placed into 50 ml centrifuge tubes at the same sampling interval as sedaDNA, and preserved in an ethanol-calcite buffer solution and stored at 4 °C on the sampling day.

To construct an in-house database of foraminiferal DNA, surface (top 1 cm) sediment samples were collected from the upper-mangrove, lower-mangrove and upper-mudflat environments along MP_A transect and from the mangrove fringe environment along transect SK (Fig. 1c; Supplementary Fig. 1). Samples were collected following the same sampling protocol as surface samples for eDNA metabarcoding. Collected samples were temporarily stored in a cool box with ice and stored in a 4 °C refrigerator in the laboratory on the sampling day.

The elevation of surface samples and the core top were measured and related to a temporary benchmark using an auto level. The benchmark elevation was measured relative to the Hong Kong Principal Datum (PD) using a Leica GS18 GNSS system. Given the reduced tidal amplitude recorded by the Solinst Levellogger® 5 at the MP_B transect⁴⁷, elevation of each surface sampling station and core top was converted to standardized water level index units (SWLI)¹¹⁹, where a value of 100 represents mean tide level (MTL) and a value of 200 represents MHHW, based on observations from the Tsim Bei Tsui tide gauge (Fig. 1c). The location of each sampling station was recorded using a handheld GPS.

## Single-cell foraminifera DNA isolation and Sanger sequencing

To construct an in-house database of foraminiferal DNA to complement the current public database, we conducted DNA barcoding with single individuals of foraminifera. For isolating single individuals, sediment samples were wet-sieved (63–500 µm) using artificial sea water (ASW, Red Sea Salt, salinity = 34) within a week of collection²⁴. The 63–500 µm size fraction, the same size fraction that we used for morphological analysis of foraminifera, was examined in ASW under a stereomicroscope for signs of vitality. Foraminifera exhibiting pigmented cytoplasm within the test and an empty last chamber were considered to be viable and placed with a bleach-sterilized brush into a picking dish with ASW for vitality assessment for maximum of 8 hours¹¹⁶,¹²⁰. Active foraminifera (exhibiting reticulopodial activity and movement) were cleaned and photographed using a Leica stereomicroscope for morphological identification. After imaging, samples were placed individually into 50 µl DOC (Deoxycholate) buffer, smashed using a sterile disposable crusher and then incubated at 60 °C for an hour¹²¹. Polymerase chain reaction (PCR) amplification was conducted using foraminifera-specific primers s14F1-sB and s14F1-N6¹²¹,¹²² following protocol described in Darling et al. (2016)¹²². For individuals that did not show any visible bands in the agarose gel electrophoresis of the nested PCR, primers s14F1-s17¹²¹ were used in a 2nd nested PCR amplification instead of s14F1-N6. Successfully amplified samples were sequenced with the Sanger method (BGI, Guangzhou). Sequences from 13 foraminiferal individuals were successfully obtained by single-individual DNA barcoding and used to construct the in-house database. Four new *Ammonia* phylotypes and one new *Cribroelphidium* phylotype were identified (Supplementary Data 3). The full sequences were uploaded to the NCBI Genbank database.

## eDNA/sedaDNA extraction, amplification, and sequencing

All eDNA/sedaDNA laboratory procedures were performed in a molecular laboratory that had no prior history of foraminiferal DNA research. All DNA extraction, PCR, and library construction procedures were conducted within a safety cabinet, with all materials pre-exposed to UV irradiation for ≥15 minutes¹²³. Except for disposable sterile or pre-autoclaved equipment, all reusable materials (e.g., pipettes) were sterilized with 70% ethanol and/or 10% bleach before use and between sample processing steps to minimize cross-contamination. DNA was extracted from surface sediment and core samples using the 0.5-g PowerSoil® DNA Isolation Kit (QIAGEN) in a sterile laboratory environment following the manufacturer protocol, and stored in a -80 °C freezer. DNA extraction was conducted in duplicate (2 technical replicates) for all surface and core samples to maximize genetic diversity²⁵. Specifically, core sample extractions were performed from the deepest to shallowest layers to minimize potential cross-contamination from upper, typically higher-concentration sediment layers. We performed the PCR by amplifying the specific hypervariable region 37 f of foraminiferal 18 s SSU rDNA using primers s14F1 and s15¹¹⁰. To multiplex PCR products for sequencing, both primers were synthesized with the overhang adapter sequences for compatibility with Illumina index and sequencing adapters following the Illumina MiSeq System sequencing protocol: 5' Adaptor—s14F1 3' and 5' Adaptor—s15 3'. For each extracted sample, PCR was conducted in duplicate (2 technical replicates) to minimize potential PCR biases¹²⁴. Positive (DNA using synthetic *Miliammina fusca* clone 4 small subunit ribosomal RNA gene as template, gBlocks®, 787 bp) and negative (blank) controls were applied during PCR amplification and library
construction to indicate positive amplification and to detect possible cross-contamination. The synthetic sequence used as a positive control was artificially modified during synthesis through nucleotide substitutions in our target barcode regions to enable contamination detection. PCR mix was prepared in 25 µl reaction, containing 2.5 µl of samples, 12.5 µl of Taq PCR Master Mix (QIAGEN), 4 µl of each primer (1 µM), 1.25 µl of Mg²⁺ (1 µM) and 0.75 µl of MilliQ water. PCR was performed as follows: an initial denaturation at 94 °C for 3 min, followed by 45 cycles of 94 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s, with a final extension at 72 °C for 5 min. The amplified products underwent a first clean-up using the AMPPure XP beads (Beckman Coulter, Singapore). A 20 µl aliquot of technical PCR replicates were pooled before first clean-up. Index PCR was then performed on the cleaned products to attach dual indices on Illumina sequencing adapters using the Nextera XT Index. The reaction mix comprises 25 µl × KAPA HiFi HotStart ReadyMix, 5 µl of each index primer (10 µM), 10 µl of PCR grade water, and 5 µL cleaned PCR product. Index PCR amplification condition was set up as follows: initial denaturation at 95 °C for 3 min, 8 cycles of denaturation for 30 s at 95 °C, annealing for 30 s at 55 °C and extension for 30 s at 72 °C and a final extension for 5 min at 72 °C. The 2nd PCR clean-up was conducted consecutively, and in total 177 cleaned foraminiferal eDNA/sedaDNA samples (including negative control, even when no visible bands were observed during agarose gel electrophoresis of the PCR products) were then pooled to 20 µM using 10 mM Tris (pH 8.5) and sent to Novogene Co, Ltd for sequencing on the NovaSeq System PE 150 platform. To minimize potential cross-contamination of the sedaDNA from modern eDNA assemblage, all laboratory procedures before the second PCR clean-up were performed separately for modern and core samples.

## Bioinformatics

Raw sequence data were demultiplexed, and primer sequences were removed by Cutadapt v2.8¹²⁵. The resulting pair-end reads were then merged by PEAR v0.9.11 with a quality threshold of 26¹²⁶. The concatenated sequence file underwent subsequent denoising using VSEARCH-UNOISE3 with default alpha diversity parameters and a minimum size threshold of 4, followed by chimera removal with VSEARCH-UCHIME¹²⁷. Sequences with length shorter than 90 bp and longer than 230 bp were discarded using VSEARCH. Sequences exhibiting >97% similarity were clustered into operational taxonomic units (OTUs)¹²⁸ by VSEARCH¹²⁷, and a resulting OTU table (numbers of OTUs per sample) containing 36,144,049 sequence number was constructed using VSEARCH-USEARCH. Quality control was applied to discard OTUs that: (1) had fewer than 10 reads¹¹⁰, or (2) appeared in negative controls¹²⁹ with read counts equal to or exceeding those in any sample¹²³, to minimize index hopping artifacts¹³⁰. Technical extraction replicates from the same surface or core sample were pooled prior to subsequent processing. Taxonomic assignments were made using BLASTn v2.12.0¹³¹ with e-value of 0.001 referencing the Genbank nucleotide (nt) database (downloaded on 13 September 2024) from the National Center for Biotechnology Information (NCBI)¹³², following previous foraminiferal eDNA studies²⁵,³⁸. The in-house database was also used as a complementary reference²⁵. The python script taxonomy_assignment_BLAST_V1.py¹³³ was used for processing the BLAST results to assign final taxonomic classifications to the species, family, and phylum level based on identity thresholds of 98%, 90%¹³⁴, and 80%¹³⁵ respectively, utilizing the top BLAST hits and best consensus taxonomy that met the predefined identity percentage thresholds. Compared to assigning OTUs to species, genus, family, and phylum levels at 99%, 95%, 90%, and 80% identity thresholds, respectively, this approach demonstrated superior performance in RSL reconstruction (Supplementary Data 12). Where BLAST results showed conflicting taxonomic assignments, OTUs were assigned at higher taxonomic ranks (genus or order level) to establish consensus identification¹²³. Where taxonomic assignments differed between GenBank and in-house databases, we adopted the in-house database classification¹²³. OTUs taxonomic assignments were conducted referring to the NCBI Taxonomy database, followed by manual verification and correction against taxonomic records in the World Register of Marine Species (WoRMS)¹³⁶. OTUs with identities between 80 to 90%

thresholds (assigned to Foraminifera phylum) and OTUs that were assigned to a foraminifera sequence with unknown taxonomy level (e.g. sequences derived from environmental sequencing) were classified as “undetermined OTUs” of foraminifera¹³⁵ and were treated as individual taxa in the subsequent analysis. OTUs with identity thresholds <80% were not considered to be foraminifera and were removed²⁴.

Although clustering sequences into amplicon sequence variants (ASVs) has become increasingly popular due to its ability to retain more genetic diversity¹³⁷, we did not use this approach. This decision was based on the large proportion of undetermined sequences in our dataset and our use of individual undetermined sequences for the analysis. Introducing higher taxonomic diversity would involve taxa with low fractional abundance, which could reduce statistical reliability in RSL reconstruction given the fixed total count¹³⁸⁻¹⁴⁰.

## Environmental variables and δ¹³C and C/N geochemistry

To assess influence of environmental parameters on foraminiferal eDNA, a suite of environmental variables was measured from modern samples: salinity, pH, and grain size. Sediment samples were centrifuged at 3000 rpm for 30 min to measure porewater salinity and pH using a Thermo Scientific™ A3255 pH/conductivity multimeter¹⁴¹. For samples that yielded insufficient porewater for measurement, 5 mL of Milli-Q water was introduced to the samples before centrifugation and measurement⁶⁴. Sediment samples for grain size measurement were treated with 30% hydrogen peroxide to remove organic carbon and then dispersed using Calgon solution. Grain size analysis was performed using a Beckman-Coulter LS 13 320 laser particle size analyzer. The software program GRADISTAT 9.1 was used to calculate grain size statistics¹⁴².

To evaluate factors influencing the composition of modern foraminiferal eDNA assemblages and to provide constraints for RSL reconstructions, stable carbon isotope geochemistry parameters—including total organic carbon (TOC), total nitrogen (TN), and δ¹³C—were measured in both modern and core samples. Prior to δ¹³C, TOC, and TN analysis, sediment samples were pretreated with 5% hydrochloric acid (HCl) for 24 h to remove carbonate material, followed by rinsing with ~1500 ml Milli-Q water⁴⁷. δ¹³C, TOC, and TN were measured at the Stable Isotope Laboratory at the University of Hong Kong using a EuroVector EA3028 Nu Horizon isotope-ratio mass spectrometer (IRMS) and a EA Isolink Elemental Analyser. δ¹³C values were normalized to the Vienna Pee Dee Belemnite (VPDB) scale with reference to triplicate analysis of certified reference materials supplied by the United States Geological Survey, USGS40 (average ¹⁵N = -4.52 and ¹³C = -26.39) and USGS41a (average ¹⁵N = 47.55 and ¹³C = 36.55).

## Morphological foraminiferal analysis

To provide context for foraminiferal assemblages inferred through eDNA and sedaDNA analyses, traditional morphological analysis was performed on surface and core sediment samples. The 63–500 µm sediment fraction was isolated through wet-sieving prior to wet-picking of foraminiferal tests under a stereomicroscope¹⁴³,¹⁴⁴. Where possible, a minimum of 100 individuals were enumerated from each sample to ensure robust transfer function development¹⁴⁰,¹⁴⁵. Taxonomy followed earlier work from the study site⁴⁷ (Supplementary Data 13). Poor specimen preservation prevented species-level identification in several genera. Therefore, species of Ammobaculites (A. exiguous, A. agglutinans), Miliammina (M. fusca, M. obliqua, M. petila), and two Ammonia species (A. beccarii, A. confertitesta) were combined under Ammobaculites spp., Miliammina spp., and Ammonia spp., respectively⁵¹,⁸⁶,⁹⁹,¹⁴⁶. Further attempts to combine A. tepida within Ammonia spp. to avoid misidentification resulted in anomalies in PME prediction (Supplementary Fig. 5), and were therefore abandoned. The morphological dataset comprises only the dead assemblage (unstained by rose Bengal) to avoid the influence of temporally variable living foraminiferal populations⁵³. To minimize bias introduced by rare taxa, only the taxa with relative abundance >5% in at least one surface or core subsample were included in the analysis⁴⁷.
¹⁴C dating

The chronology of core MPSC01 was constructed from radiocarbon-dated plant fragments and fine-fraction organic mud from 15 core samples (Supplementary Data 14). Plant fragments, including mangrove bark, branches and twigs deemed to be deposited horizontally on the former sediment surface with no clear evidence of reworking, were chosen as dating materials²¹,¹⁰³. Plant fragments were examined and cleaned with Milli-Q water under a stereomicroscope to remove all visible rootlets and adhering material. Fine-fraction organic mud samples for radiocarbon dating were prepared by isolating the <63-µm fraction on a 0.47 µm quartz fiber filter under filtration, followed by oven drying for 24 h at 40 °C. Dried organic mud samples were weighed, wrapped in aluminum foil, placed in clean glass jars, and stored in a desiccator until further processing. To minimize modern carbon contamination and cross-contamination, all equipment was pre-treated before use. Quartz fiber filters and aluminum foil were baked at 550 °C for 4 h, then cooled in a desiccator. All sieves and processing tools were cleaned between samples using a four-step protocol: (1) washing with soapy water; (2) rinsing with tap water; (3) soaking in 10% HCl for ~30 seconds; and (4) triple rinsing with Milli-Q water. Additionally, the 63-µm sieve was sonicated for 15 min before processing each sample. Samples were submitted either to Beta Analytic in Miami, FL, USA, using accelerator mass spectrometry (AMS) dating, or to the Yale Analytical and Stable Isotope Center (YASIC) in New Haven, CT, USA, using the MICADAS AMS system.

An age model for core MPSC01 was constructed using the Bacon age-depth model¹⁴⁷, incorporating all radiocarbon-dated samples to estimate the age with $\sigma$ uncertainty of each eDNA/morphological core sample (Supplementary Fig. 4; Supplementary Data 14). The Intcal20 calibration curve¹⁴⁸ was used for calibrating the conventional ages to calendar years, and the post-1950 CE samples were calibrated using the bomb spike dataset¹⁴⁹.

Alpha diversity analysis

To estimate the $\alpha$-diversity of eDNA and sedaDNA assemblages, rarefaction¹⁵⁰ was performed to standardize each sample to the lowest sequencing depth of the dataset (10,591 reads), a depth at which most samples reached a plateau in rarefaction curves (Supplementary Fig. 6). Two samples (one modern and one core sample) with insufficient read counts were excluded from further analysis. For the remaining samples, three $\alpha$-diversity metrics were calculated: observed OTUs¹³⁷, Shannon index¹⁵¹ and Chao1 index¹⁵². Differences in $\alpha$-diversity were assessed using one-way ANOVA with 1,000 permutations (*p*-value threshold = 0.05)¹²³,¹⁵³, focusing on: (1) variation among modern samples across six identified environments; (2) differences between modern and core samples; and (3) differences between core samples from Units II/ III and Unit I.

Similarly, $\alpha$-diversity metrics—including taxa richness, Shannon index, and Chao1 index—were calculated and then assessed with one-way ANOVA for the morphological assemblage, with a focus on: (1) differences among samples from six modern environments, and (2) comparisons between modern and core samples. Levene's test confirmed homogeneity of variances across all groups (*p* > 0.05) prior to ANOVA¹⁵⁴.

Multivariate Ordination analysis

Non-Metric Multidimensional Scaling (NMDS) analysis was conducted on modern foraminiferal eDNA and morphological assemblages to highlight dissimilarities in terms of community composition between sampling stations with different elevations¹⁵⁵. NMDS was performed using a Bray-Curtis dissimilarity matrix based on taxa relative abundance data. The resulting ordination plot displays the first two NMDS dimensions, where the spatial arrangement of points reflects the compositional dissimilarity among stations. A Permutational Multivariate Analysis of Variance (PERMANOVA) with 999 permutations, based on Bray-Curtis distance, was performed to assess whether the observed grouping patterns in the NMDS ordination were significantly explained by elevation differences¹⁵⁶. Analysis was conducted using the “vegan” package (v2.6-4) in R (v4.3.1)¹⁵⁷.

Canonical correspondence analysis (CCA) was performed with relative abundance data to examine the correlation between the contemporary foraminiferal eDNA and morphological assemblages with environmental variables across all modern samples¹⁵⁸, also conducted using the “vegan” package (v2.6-4) in R¹⁵⁷. For both eDNA and morphological assemblages, we included only taxa that met two criteria in the CCA: (1) > 5% relative abundance in at least one modern or core sample, and (2) presence in more than one sample within the modern training set⁴⁷,⁵⁰. In CCA and partial CCA (pCCA), the degree of variation explained by environmental variables are measured by canonical eigenvalues⁵⁰,⁶⁵.

Transfer function construction

A Bayesian transfer function (BTF) was constructed using modern training sets of foraminiferal eDNA and morphological assemblages. Taxa with >5% abundance of the assemblage in any sample were included in the model to minimize uncertainty associated with low-abundance taxa⁹⁸,¹³⁸,¹³⁹. The BTF approach uses penalized spline functions to model non-parametric species response curves, enabling the representation of complex, multimodal, and non-Gaussian responses to elevation. Informative priors derived from elevation-dependent environments based on stable carbon isotope geochemistry⁹⁸ were incorporated to constrain predictions.

For the morphological assemblage, counts of foraminiferal taxa served as model input, excluding the core sample at 1.55-1.56 m PD due to insufficient counts. This approach reduces uncertainty associated with larger count sizes while maintaining flexibility to capture multi-modal and non-Gaussian species responses⁹⁸. For the eDNA assemblage, raw sequence counts were not used due to the discrepancies in sequencing depth between samples¹⁵⁹. Instead, rarefied sequences—subsampled to a uniform sequencing depth¹⁵⁰—were used as model input, as in the $\alpha$-diversity calculation. While rarefying DNA data requires caution¹⁶⁰, compared to normalization methods¹⁶¹, this approach preserves the discrete nature of the dataset and avoids violating the BTF assumption of multinomial taxa distribution⁹⁸,¹⁶².

Ten-fold cross-validation was conducted to evaluate BTF performance for eDNA and morphological assemblages by comparing the root mean squared error of prediction (RMSEP), correlation between predicted and observed values ($R^2$), average $\sigma$ uncertainty, average residual and average absolute residual between predicted and observed values of surface samples⁴⁸. A single-pass approach was adopted to screen out samples with residuals exceeding two standard deviations from the mean during BTF construction¹⁶³.

Given the lack of established best practices for developing BTFs from eDNA data, sensitivity tests were conducted to assess how training set composition influenced BTF performance. These tests included (1) the use of rarefied vs. non-rarefied taxa, (2) using different abundance thresholds for removing rare taxa (e.g., 1, 2, 5%), (3) using different taxonomy assignment criteria (e.g., using 99% and 95% identity for species- and genus-level assignment), (4) inclusion vs. exclusion of undetermined OTUs, (5) inclusion vs. exclusion of soft-walled monothalamids, (6) using training sets with stronger elevation signals (e.g., mangrove-only or mangrove and mudflat subsets), and (7) the use of site-specific training sets. These tests confirmed that excluding low-abundance taxa through the application of an abundance threshold¹³⁸,¹³⁹—for which 5% was determined to be optimal in our study—mitigates the potential impact of rare sequence loss during rarefaction, ensuring robust BTF performance. Using OTUs clustered at the lowest taxonomic level allowed input taxa to represent varying taxonomic levels, while treating undetermined OTUs as individual taxa and including monothalamids further improved model performance. Additionally, the training set incorporating samples across all environments yielded the most precise and accurate RSL reconstructions (Supplementary Fig. 7; Supplementary Data 12). The eDNA-BTF also demonstrated better performance in cross-validation and higher accuracy compared to the widely used weighted-averaging partial least squares regression (WA-PLS) transfer function (Supplementary Data 15), providing an additional rationale for selecting this BTF approach.
## Linear Discriminant Analysis (LDA) Priors

Linear discriminant analysis (LDA) was applied to core samples from MPSC01 to estimate prior elevation ranges used in the BTF⁹⁸. Linear discriminant functions were derived from geochemical data (δ¹³C, TOC, and C/N) collected from modern surface samples across six elevation-dependent environmental zones⁴⁷, including subtidal, mudflat, mangrove fringe, mangrove, terrestrial transition, and terrestrial environments (Supplementary Fig. 1; Supplementary Data 1). They were then applied to core samples to estimate the likelihood of each sample belonging to the pre-defined environmental zones. Core samples were assigned to a single environmental zone if their probability of affiliation exceeded a threshold of >0.9. Samples with probabilities <0.9 were assigned to multiple environmental zones. The observed elevation range corresponding to each environmental zone was used as prior information in the BTF (Supplementary Data 16).

## Application of Bayesian Transfer Function (BTF) to Core Samples

The eDNA- and morphology-BTFs were applied to core MPSC01 to reconstruct paleo-mangrove elevation (PME). The performance and reliability of the BTF models was evaluated using: (1) 1σ uncertainty derived from core prediction results, and (2) dissimilarity and goodness-of-fit thresholds between modern and fossil assemblages and their elevation⁹⁶ (Fig. 4).

Dissimilarity between modern and core samples was assessed using Bray-Curtis metrics⁹⁶, consistent with the approach applied in the NMDS analysis. Thresholds of 20% dissimilarity for morphological assemblages and 10% for eDNA assemblages were applied, reflecting the higher diversity and heterogeneity of eDNA datasets²⁰. Core samples exceeding these thresholds were interpreted cautiously due to the lack of modern counterparts⁹⁶. Goodness-of-fit statistics were calculated by fitting core samples into a constrained ordination (CCA) derived from the modern training set to evaluate their fit to elevation⁹⁶. A 95% threshold of squared residual lengths was used to identify poor fits.

## Validation of RSL reconstructions using regional tide gauge and geologic data

RSL was calculated by subtracting the PME reconstructed by BTFs from the elevation of each core sample⁶⁹, both expressed relative to the Principal Datum of Hong Kong¹⁶⁴. Vertical uncertainties were derived from sample-specific 2σ uncertainties from the BTFs. A chronology for RSL changes was obtained from the age-depth model (Section “¹⁴C dating”).

Reconstructed RSL was compared with tide gauge data from Hong Kong to validate accuracy. Due to the rapid sedimentation rate observed in the upper portion of MPSC01, annual average tide gauge measurements were calculated relative to a 5-year average (2017–2022 CE) from tide gauges at Tai Po Kau, Tai Miu Wan, Tsim Bei Tsui, Shek Pik, and Quarry Bay⁴⁷ (Supplementary Fig. 8). This comparison spans reconstructed RSL and observed measurements from 1963 to 2022 (Figs. 6a, b). The mean squared error (MSE) between reconstructed RSL midpoints and tide gauge observations was calculated to assess reconstruction accuracy (e.g., Cahill et al., 2016)⁹⁸. For older samples in age than the tide gauge records, the eDNA-BTF reconstruction was evaluated by comparison with regional sea-level index points (defining RSL position and age) and terrestrial/marine limiting points (indicating upper and lower RSL limits)¹¹⁹ from the PRD⁴⁶,¹⁶⁵, to confirm its accuracy.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

All Supplementary Data are provided in the Excel file entitled “Supplementary Data” deposited on Figshare (https://doi.org/10.6084/m9 figshare.30405109). All Supplementary Figs. and corresponding legends are presented in “Supplementary Information” file. Raw sequence data of eDNA samples this study have been deposited in National Center for Biotechnology

Information (NCBI) Sequence Read Archive (SRA) under the BioProject ID PRJNA1273657 (https://dataview.ncbi.nlm.nih.gov/object/PRJNA1273657?reviewer=jdt0fihub4idpug78ea1k1dkk9). All barcoded sequences have been deposited in NCBI Genbank with accession numbers PV764368-PV764380.

## Code availability

Codes of bioinformatic pipeline and eDNA Bayesian transfer function are available from the GitHub repository: https://github.com/FrancisLauzj-hku/Forams_meta barcoding_sea_level.

Received: 16 August 2025; Accepted: 7 November 2025;

Published online: 23 November 2025

## References

1. Dutton, A. et al. Sea-level rise due to polar ice-sheet mass loss during past warm periods. *Science* **349**, 6244 (2015).

2. Gehrels, W. R., Belknap, D. F. & Kelley, J. T. Integrated high-precision analyses of Holocene relative sea-level changes: lessons from the coast of Maine. *Geol. Soc. Am. Bull.* **108**, 1073–1088 (1996).

3. Engelhart, S. E., Horton, B. P., Douglas, B. C., Peltier, W. R. & Törnqvist, T. E. Spatial variability of late Holocene and 20th century sea-level rise along the Atlantic coast of the United States. *Geology* **37**, 1115–1118 (2009).

4. Lambeck, K. et al. Sea level change along the Italian coast during the Holocene and projections for the future. *Quatern Int* **232**, 250–257 (2011).

5. Walker, J. S. et al. Common Era sea-level budgets along the US Atlantic coast. *Nat. Commun.* **12**, 1841 (2021).

6. Kopp, R. E. et al. Temperature-driven global sea-level variability in the Common Era. *Proc. Natl. Acad. Sci.* **113**, E1434–E1441 (2016).

7. Gehrels, W. R. & Woodworth, P. L. When did modern rates of sea-level rise start? *Glob. Planet Change* **100**, 263–277 (2013).

8. Garrett, E. et al. Drivers of 20th century sea-level change in southern New Zealand determined from proxy and instrumental records. *J. Quat. Sci.* **37**, 1025–1043 (2022).

9. Rahmstorf, S. A semi-empirical approach to projecting future. *Science* **315**, 368 (2004).

10. Van de Plassche, O. *Sea-level Research: a Manual for the Collection and Evaluation of Data: a Manual for the Collection and Evaluation of Data*. (Springer, 1986).

11. Patterson, R. T. et al. Relative utility of foraminifera, diatoms and macrophytes as high resolution indicators of paleo-sea level in coastal British Columbia, Canada. *Quat. Sci. Rev.* **24**, 2002–2014 (2005).

12. Papaspyrou, S., Diz, P., García-Robledo, E., Corzo, A. & Jimenez-Arias, J.-L. Benthic foraminiferal community changes and their relationship to environmental dynamics in intertidal muddy sediments (Bay of Cádiz, SW Spain). *Mar. Ecol. Prog. Ser.* **490**, 121–135 (2013).

13. Woodroffe, S. A. Recognising subtidal foraminiferal assemblages: implications for quantitative sea-level reconstructions using a foraminifera-based transfer function. *J. Quat. Sci.* **24**, 215–223 (2009).

14. Berkeley, A., Perry, C. T., Smithers, S. G., Horton, B. P. & Cundy, A. B. Foraminiferal biofacies across mangrove-mudflat environments at Cocoa Creek, north Queensland, Australia. *Mar. Geol.* **263**, 64–86 (2009).

15. Scott, D. & Medioli, F. Vertical zonations of marsh foraminifera as accurate indicators of former sea-levels. *Nature* **272**, 528–531 (1978).

16. Rush, G. et al. Development of an intertidal foraminifera training set for the North Sea and an assessment of its application for Holocene sea-level reconstructions. *Mar. Micropaleontol.* **169**, 102055 (2021).

17. Barnett, R. L., Garneau, M. & Bernatchez, P. Salt-marsh sea-level indicators and transfer function development for the Magdalen Islands in the Gulf of St. Lawrence, Canada. *Mar. Micropaleontol.* **122**, 13–26 (2016).
18. Rossi, V. et al. The application of foraminifera to reconstruct the rate of 20th century sea level rise, Morbihan Golfe, Brittany, France. *Quat. Res.* **75**, 24–35 (2011).

19. Berkeley, A., Perry, C. T., Smithers, S. G., Horton, B. & Taylor, K. G. A review of the ecological and taphonomic controls on foraminiferal assemblage development in intertidal environments. *Earth-Sci. Rev.* **83**, 205–230 (2007).

20. Berkeley, A., Perry, C. T. & Smithers, S. G. Taphonomic signatures and patterns of test degradation on tropical, intertidal benthic foraminifera. *Mar. Micropaleontol.* **73**, 148–163 (2009).

21. Khan, N. S. et al. Relative sea-level change in South Florida during the past-5000 years. *Glob. Planet Change* **216**, 103902 (2022).

22. Taberlet, P., Bonin, A., Zinger, L. & Coissac, E. *Environmental DNA: For *biodiversity research and monitoring*. (Oxford University Press, 2018).

23. He, X., Sutherland, T. F., Pawlowski, J. & Abbott, C. L. Responses of foraminifera communities to aquaculture-derived organic enrichment as revealed by environmental DNA metabarcoding. *Mol. Ecol.* **28**, 1138–1153 (2019).

24. Brinkmann, I. et al. Through the eDNA looking glass: Responses of fjord benthic foraminiferal communities to contrasting environmental conditions. *J. Eukaryot. Microbiol.* **70**, e12975 (2023).

25. Singer, D. et al. Unlocking foraminiferal genetic diversity on estuarine mudflats with eDNA metabarcoding. *Sci. Total Environ.* **902**, 165983 (2023).

26. Pawlowski, J. et al. Benthic monitoring of salmon farms in Norway using foraminiferal metabarcoding. *Aquacult Env Interac* **8**, 371–386 (2016).

27. Lejzerowicz, F., Esling, P. & Pawlowski, J. Patchiness of deep-sea benthic Foraminifera across the Southern Ocean: Insights from high-throughput DNA sequencing. *Deep-Sea Res Pt II* **108**, 17–26 (2014).

28. Holzmann, M., Gooday, A. J., Majewski, W. & Pawlowski, J. Molecular and morphological diversity of monothalamous foraminifera from South Georgia and the Falkland Islands: description of four new species. *Eur. J. Protistol.* **85**, 125909 (2022).

29. Crecchio, C. & Stotzky, G. Binding of DNA on humic acids: effect on transformation of Bacillus subtilis and resistance to DNase. *Soil Biol. Biochem.* **30**, 1061–1067 (1998).

30. Blum, S. A., Lorenz, M. G. & Wackernagel, W. Mechanism of retarded DNA degradation and prokaryotic origin of DNases in nonsterile soils. *Syst. Appl. Microbiol.* **20**, 513–521 (1997).

31. Freeman, C. et al. Survival of environmental DNA in sediments: Mineralogic control on DNA taphonomy. *Environ. DNA* **5**, 1691–1705 (2023).

32. Levy-Booth, D. J. et al. Cycling of extracellular DNA in the soil environment. *Soil Biol. Biochem.* **39**, 2977–2991 (2007).

33. Bremond, L. et al. Five thousand years of tropical lake sediment DNA records from Benin. *Quat. Sci. Rev.* **170**, 203–211 (2017).

34. Liu, S. et al. Tibetan terrestrial and aquatic ecosystems collapsed with cryosphere loss inferred from sedimentary ancient metagenomics. *Sci. Adv.* **10**, eadn8490 (2024).

35. Boere, A. et al. Late-Holocene succession of dinoflagellates in an Antarctic fjord using a multi-proxy approach: Paleoenvironmental genomics, lipid biomarkers and palynomorphs. *Geobiology* **7**, 265–281 (2009).

36. Ekram, M. A. E. et al. A 1 Ma sedimentary ancient DNA (sed aDNA) record of catchment vegetation changes and the developmental history of tropical Lake Towuti (Sulawesi, Indonesia). *Geobiology* **22**, e12599 (2024).

37. Ekram, M. A.-E. et al. A Quaternary sedimentary ancient DNA (sedaDNA) record of fungal-terrestrial ecosystem dynamics in a tropical biodiversity hotspot (Lake Towuti, Sulawesi, Indonesia). *Microorganisms* **13**, 1005 (2025).

38. Angeles, I. B. et al. Planktonic foraminifera eDNA signature deposited on the seafloor remains preserved after burial in marine sediments. *Sci. Rep.-Uk* **10**, 20351 (2020).

39. Pawłowska, J. et al. Palaeoceanographic changes in Hornsund Fjord (Spitsbergen, Svalbard) over the last millennium: new insights from ancient DNA. *Climate* **12**, 1459–1472 (2016).

40. Pawłowska, J. et al. Ancient DNA sheds new light on the Svalbard foraminiferal fossil record of the last millennium. *Geobiology* **12**, 277–288 (2014).

41. Lejzerowicz, F. et al. Ancient DNA complements microfossil record in deep-sea subsurface sediments. *Biol. Lett.* **9**, 20130283 (2013).

42. Mojtahid, M. et al. Combining morphological and molecular data to study past foraminiferal communities from a temperate coastal sediment core. *Boreas* (2025).

43. Angeles, I. B. et al. Encapsulated in sediments: eDNA deciphers the ecosystem history of one of the most polluted European marine sites. *Environ. Int.* **172**, 107738 (2023).

44. Schweizer, M. How long after death is DNA preserved in situ in intertidal foraminifera? *J. Micropaleontolgy* **34**, 217–219 (2015).

45. Mejbel, H. S., Dodsworth, W. & Pick, F. R. Effects of temperature and oxygen on cyanobacterial DNA preservation in sediments: A comparison study of major taxa. *Environ. DNA* **4**, 717–731 (2022).

46. Sun, Y., Xiong, H., Lee, M.-T., Brodie, C. & Zong, Y. Geochemical dynamics and depositional history from mangrove sediments within the Pearl River estuary. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **584**, 11 (2021).

47. Yu, H. K. Y. et al. The utility of mangrove foraminifera, diatoms, and stable carbon isotope and C/N geochemistry in relative sea-level reconstruction in the Pearl River Delta, China. *Mar. Geol.* **488**, 107610 (2025).

48. Kemp, A. C. & Telford, R. J. in *Handbook of sea-level research, Transfer functions*, 470-499 (Wiley, 2015).

49. Juggins, S. Quantitative reconstructions in palaeolimnology: new paradigm or sick science?. *Quat. Sci. Rev.* **64**, 20–32 (2013).

50. Horton, B. P. & Edwards, R. J. *Quantifying Holocene sea level change using intertidal foraminifera: lessons from the British Isles*. (Cushman Foundation for foraminiferal research, 2006).

51. Kemp, A. C., Horton, B. P. & Culver, S. J. Distribution of modern salt-marsh foraminifera in the Albemarle–Pamlico estuarine system of North Carolina, USA: Implications for sea-level research. *Mar. Micropaleontol.* **72**, 222–238 (2009).

52. Kemp, A. C. et al. Quantitative vertical zonation of salt-marsh foraminifera for reconstructing former sea level; an example from New Jersey, USA. *Quat. Sci. Rev.* **54**, 26–39 (2012).

53. Horton, B. The distribution of contemporary intertidal foraminifera at Cowpen Marsh, Tees Estuary, UK: implications for studies of Holocene sea-level changes. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **149**, 127–149 (1999).

54. Ellegaard, M. et al. Dead or alive: sediment DNA archives as tools for tracking aquatic evolution and adaptation. *Commun. Biol.* **3**, 169 (2020).

55. Goldstein, S. T. & Alve, E. Experimental assembly of foraminiferal communities from coastal propagule banks. *Mar. Ecol. Prog. Ser.* **437**, 1–11 (2011).

56. Hassan, G. S., Espinosa, M. A. & Isla, F. I. Diatom-based inference model for paleosalinity reconstructions in estuaries along the northeastern coast of Argentina. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **275**, 77–91 (2009).

57. Zong, Y. & Horton, B. P. Diatom-based tidal-level transfer functions as an aid in reconstructing Quaternary history of sea-level movements in the UK. *J. Quat. Sci.* **14**, 153–167 (1999).

58. Horton, B. P., Corbett, R., Culver, S. J., Edwards, R. J. & Hillier, C. Modern saltmarsh diatom distributions of the Outer Banks, North Carolina, and the development of a transfer function for high resolution reconstructions of sea level. *Estuar., Coast. Shelf Sci.* **69**, 381–394 (2006).

59. Reef, R. et al. Using eDNA to determine the source of organic carbon in seagrass meadows. *Limnol. Oceanogr.* **62**, 1254–1265 (2017).
60. Le Cadre, V. R., Debenay, J.-P. & Lesourd, M. Low pH effects on Ammonia beccarii test deformation: implications for using test deformations as a pollution indicator. *J. Foraminifer. Res.* **33**, 1–9 (2003).

61. Angeles, I. B., Nguyen, N.-L., Greco, M., Tan, K. S. & Pawlowski, J. Assigning the unassigned: A signature-based classification of rDNA metabarcodes reveals new deep-sea diversity. *Plos One* **19**, e0298440 (2024).

62. Santos, S. S. et al. Soil DNA extraction procedure influences protist 18S rRNA gene community profiling outcome. *Protist* **168**, 283–293 (2017).

63. Peng, Q., Vijaya Satya, R., Lewis, M., Randad, P. & Wang, Y. Reducing amplification artifacts in high multiplex amplicon sequencing by using molecular barcodes. *BMC genomics* **16**, 1–12 (2015).

64. Hong, I. et al. Diatoms of the intertidal environments of Willapa Bay, Washington, USA as a sea-level indicator. *Mar. Micropaleontol.* **167**, 102033 (2021).

65. Sawai, Y., Horton, B. P. & Nagumo, T. The development of a diatom-based transfer function along the Pacific coast of eastern Hokkaido, northern Japan—an aid in paleoseismic studies of the Kuril subduction zone. *Quat. Sci. Rev.* **23**, 2467–2483 (2004).

66. Zong, Y. Mid-and late-Holocene sea-level changes in Roudsea Marsh, northwest England: a diatom biostratigraphical investigation. *Holocene* **7**, 311–323 (1997).

67. Walker, J. S. et al. Incorporating temporal and spatial variability of salt-marsh foraminifera into sea-level reconstructions. *Mar. Geol.* **429**, 106293 (2020).

68. Gehrels, W. R. Using foraminiferal transfer functions to produce high-resolution sea-level records from salt-marsh deposits, Maine, USA. *Holocene* **10**, 367–376 (2000).

69. Kemp, A. C. et al. The relative utility of foraminifera and diatoms for reconstructing late Holocene sea-level change in North Carolina, USA. *Quat. Res.* **71**, 9–21 (2009).

70. Watcham, E. P., Shennan, I. & Barlow, N. L. Scale considerations in using diatoms as indicators of sea-level change: Lessons from Alaska. *J. Quat. Sci.* **28**, 165–179 (2013).

71. Gehrels, W. R., Roe, H. M. & Charman, D. J. Foraminifera, testate amoebae and diatoms as sea-level indicators in UK saltmarshes: a quantitative multiproxy approach. *J. Quat. Sci.: Published Quat. Res. Assoc.* **16**, 201–220 (2001).

72. Gooday, A. J. Organic-walled allogromiids: aspects of their occurrence, diversity and ecology in marine habitats. *J. Foraminifer. Res.* **32**, 384–399 (2002).

73. Sabbatini, A. et al. Modern benthic foraminifers at Northern shallow sites of Adriatic Sea and soft-walled, monothalamous taxa: a brief overview. *Micropaleontology* **56**, 359–376 (2010).

74. Henderson, Z. Soft-walled monothalamid foraminifera from the intertidal zones of the Lorn area, north-west Scotland. *J. Mar. Biol. Assoc. U. Kingd.* **103**, e18 (2023).

75. Berkeley, A., Perry, C., Smithers, S. & Horton, B. The spatial and vertical distribution of living (stained) benthic foraminifera from a tropical, intertidal environment, north Queensland, Australia. *Mar. Micropaleontol.* **69**, 240–261 (2008).

76. Buzas, M. A., C. Hayek, L.-A., Jett, J. A. & Reed, S. A. *Pulsating patches: History and Analyses of spatial, seasonal, and yearly distribution of living benthic foraminifera*. Vol. 97 (Smithsonian Institution Scholarly Press, 2015).

77. Debenay, J.-P., Bicchi, E., Goubert, E. & Du Châtelet, E. A. Spatiotemporal distribution of benthic foraminifera in relation to estuarine dynamics (Vie estuary, Vendée, W France). *Estuar., Coast. Shelf Sci.* **67**, 181–197 (2006).

78. De Schepper, S. et al. The potential of sedimentary ancient DNA for reconstructing past sea ice evolution. *ISME J.* **13**, 2566–2577 (2019).

79. Hebsgaard, M. B., Phillips, M. J. & Willerslev, E. Geologically ancient DNA: fact or artefact?. *Trends Microbiol.* **13**, 212–220 (2005).

80. Hofreiter, M., Jaenicke, V., Serre, D., Haeseler, A. V. & Pääbo, S. DNA sequences from multiple amplifications reveal artifacts induced by cytosine deamination in ancient DNA. *Nucleic acids Res.* **29**, 4793–4799 (2001).

81. Demanèche, S., Jocteur-Monrozier, L., Quiquampoix, H. & Simonet, P. Evaluation of biological and physical protection against nuclease degradation of clay-bound plasmid DNA. *Appl. Environ. Microbiol.* **67**, 293–299 (2001).

82. Tyson, R. Sedimentation rate, dilution, preservation and total organic carbon: some results of a modelling study. *Org. Geochem.* **32**, 333–339 (2001).

83. Le Moigne, F. A., Cisternas-Novoa, C., Piontek, J., Maßmig, M. & Engel, A. On the effect of low oxygen concentrations on bacterial degradation of sinking particles. *Sci. Rep.-Uk* **7**, 16722 (2017).

84. Corinaldesi, C., Beolchini, F. & DELL'ANNO, A. Damage and degradation rates of extracellular DNA in marine sediments: implications for the preservation of gene sequences. *Mol. Ecol.* **17**, 3939–3951 (2008).

85. Jo, T. et al. Rapid degradation of longer DNA fragments enables the improved estimation of distribution and biomass using environmental DNA. *Mol. Ecol. Resour.* **17**, e25–e33 (2017).

86. Culver, S. J. & Horton, B. P. Infaunal marsh foraminifera from the outer banks, North Carolina, USA. *J. Foraminifer. Res.* **35**, 148–170 (2005).

87. Goldstein, S. T. & Watkins, G. T. Taphonomy of salt marsh foraminifera: an example from coastal Georgia. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **149**, 103–114 (1999).

88. Murray, J. W. & Alve, E. Natural dissolution of modern shallow water benthic foraminifera: taphonomic effects on the palaeoecological record. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **146**, 195–209 (1999).

89. Rogers, K. et al. Wetland carbon storage controlled by millennial-scale variation in relative sea-level rise. *Nature* **567**, 91–95 (2019).

90. Xiong, H., Zong, Y., Qian, P., Huang, G. & Fu, S. Holocene sea-level history of the northern coast of South China Sea. *Quat. Sci. Rev.* **194**, 12–26 (2018).

91. Willerslev, E. et al. Fifty thousand years of Arctic vegetation and megafaunal diet. *Nature* **506**, 47–51 (2014).

92. Baisheva, I. et al. Permafrost-thaw lake development in Central Yakutia: sedimentary ancient DNA and element analyses from a Holocene sediment record. *J. Paleolimnol.* **70**, 1–18 (2023).

93. Harðardóttir, S. et al. Millennial-scale variations in Arctic sea ice are recorded in sedimentary ancient DNA of the microalga Polarella glacialis. *Commun. Earth Environ.* **5**, 25 (2024).

94. Kemp, A. C. et al. Climate related sea-level variations over the past two millennia. *Proc. Natl. Acad. Sci.* **108**, 11017–11022 (2011).

95. Gehrels, W. R. et al. Nineteenth and twentieth century sea-level changes in Tasmania and New Zealand. *Earth Planet. Sci. Lett.* **315**, 94–102 (2012).

96. Kemp, A. C., Telford, R. J., Horton, B. P., Anisfeld, S. C. & Sommerfield, C. K. Reconstructing Holocene sea level using salt-marsh foraminifera and transfer functions: lessons from New Jersey, USA. *J. Quat. Sci.* **28**, 617–629 (2013).

97. Juggins, S. & Birks, H. J. B. Quantitative environmental reconstructions from biological data. *Tracking environmental change using lake sediments: Data handling and numerical techniques*, 431–494 (2012).

98. Cahill, N., Kemp, A. C., Horton, B. P. & Parnell, A. C. Bayesian hierarchical model for reconstructing relative sea level: from raw data to rates of change. *Climate* **12**, 525–542 (2016).

99. Khan, N. S., Vane, C. H., Engelhart, S. E., Kendrick, C. & Horton, B. P. The application of δ13C, TOC and C/N geochemistry of mangrove
sediments to reconstruct Holocene paleoenvironments and relative
sea levels, Puerto Rico. Mar. Geol. 415, 105963 (2019).

100. Khan, N. S. et al. Drivers of Holocene sea-level change in the Caribbean. Quat. Sci. Rev. 155, 13–36 (2017).

101. Macintyre, I. G., Toscano, M. A. & Bond, G. Holocene history of the mangrove islands of Twin Cays, Belize, Central America. Atoll Research Bulletin (2004).

102. Tan, F. et al. Holocene relative sea-level histories of far-field islands in the mid-Pacific. Quat. Sci. Rev. 310, 107995 (2023).

103. Kemp, A. C. et al. Late Holocene sea-and land-level change on the US southeastern Atlantic coast. Mar. Geol. 357, 90–100 (2014).

104. Barlow, N. L. et al. Salt marshes as late Holocene tide gauges. Glob. Planet Change 106, 90–110 (2013).

105. Kemp, A. C. et al. Relative sea-level change in Newfoundland, Canada during the past~3000 years. Quat. Sci. Rev. 201, 89–110 (2018).

106. Nguyen, N. L., Pawłowska, J., Angeles, I. B., Zajaczkowski, M. & Pawłowski, J. Metabarcoding reveals high diversity of benthic foraminifera linked to water masses circulation at coastal Svalbard. Geobiology 21, 133–150 (2023).

107. Alve, E. & Goldstein, S. T. The propagule method as an experimental tool in foraminiferal ecology. Approaches to study living foraminifera: collection, maintenance and experimentation, 1-12 (2014).

108. Armbrecht, L. H. et al. Ancient DNA from marine sediments: Precautions and considerations for seafloor coring, sample handling and data generation. Earth-Sci Rev 196 (2019).

109. Brown, T. A. & Barnes, I. M. The current and future applications of ancient DNA in Quaternary science. J. Quat. Sci. 30, 144–153 (2015).

110. Pawlowski, J., Esling, P., Lejzerowicz, F., Cedhagen, T. & Wilding, T. A. Environmental monitoring through protist next-generation sequencing metabarcoding: assessing the impact of fish farming on benthic foraminifera communities. Mol. Ecol. Resour. 14, 1129–1140 (2014).

111. Hong Kong Observatory. Summary of meteorological and tidal observations in Hong Kong, https://www.hko.gov.hk/en/publica/ pubsmo.htm (2024).

112. Environmental Protection Department of Hong Kong. Marine water quality in Hong Kong in 2023, https://www.epd.hk/epd/english/ environmentinhk/water/data/achi_wq_obj.html (2023).

113. Li, M. & Lee, S. Carbon dynamics of Deep Bay, eastern Pearl River Estuary, China. I: A mass balance budget and implications for shorebird conservation. Mar. Ecol. Prog. Ser. 172, 73–87 (1998).

114. Lee, S. Carbon dynamics of Deep Bay, eastern Pearl River estuary, China. II: Trophic relationship based on carbon-and nitrogen-stable isotopes. Mar. Ecol. Prog. Ser. 205, 1–10 (2000).

115. Stoeck, T. et al. Environmental DNA metabarcoding of benthic bacterial communities indicates the benthic footprint of salmon aquaculture. Mar. Pollut. Bull. 127, 139–149 (2018).

116. Schweizer, M. et al. Trophic strategies of intertidal foraminifera explored with single-cell microbiome metabarcoding and morphological methods: What is on the menu?. Ecol. Evolution 12, e9437 (2022).

117. Walton, W. R. Techniques for recognition of living foraminifera. Contribution Cushman Found. Foraminifer. Res. 3, 56–60 (1952).

118. Foster, N. R. et al. Environmental DNA identifies coastal plant community shift 1,000 years ago in Torrens Island, South Australia. Commun. Earth Environ. 5, 115 (2024).

119. Shennan, I., Long, A. J. & Horton, B. P. Handbook of Sea-level Research. (Wiley, 2015).

120. Schweizer, M., Pawlowski, J., Duijnstee, I. A., Kouwenhoven, T. & Van der Zwaan, G. Molecular phylogeny of the foraminiferal genus Uvigerina based on ribosomal DNA sequences. Mar. Micropaleontol. 57, 51–67 (2005).

121. Pawlowski, J. Introduction to the molecular systematics of foraminifera. Micropaleontology 46, 1–12 (2000).

122. Darling, K. F. et al. The genetic diversity, phylogeography and morphology of Elphidiidae (Foraminifera) in the Northeast Atlantic. Mar. Micropaleontol. 129, 1–23 (2016).

123. Chung, A., Kam, Y., Shea, S. K. & Schunter, C. Detecting Fish Diversity in Urban-Impacted Ecosystems: A Comparative Approach of eDNA Metabarcoding and UVC. Environ. DNA 6, e70048 (2024).

124. Nichols, R. V. et al. Minimizing polymerase biases in metabarcoding. Mol. Ecol. Resour. 18, 927–939 (2018).

125. Martin, M. Cutadapt removes adapter sequences from high- throughput sequencing reads. EMBnet. J. 17, 10–12 (2011).

126. Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30, 614–620 (2014).

127. Rognes, T., Flouri, T., Nichols, B., Quince, C. & Mahé, F. VSEARCH: a versatile open source tool for metagenomics. PeerJ 4, e2584 (2016).

128. Salonen, I. S., Chronopoulou, P.-M., Bird, C., Reichart, G.-J. & Koho, K. A. Enrichment of intracellular sulphur cycle-associated bacteria in intertidal benthic foraminifera revealed by 16S and aprA gene analysis. Sci. Rep.-Uk 9, 11692 (2019).

129. Edmonds, K. & Williams, L. The role of the negative control in microbiome analyses. FASEB J. 31, 940 (2017).

130. Hornung, B. V., Zwittink, R. D. & Kuijper, E. J. Issues and current standards of controls in microbiome research. *Fems Microbiol Ecol.* **95**, fiz045 (2019).

131. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).

132. Sayers, E. W. et al. GenBank 2024 update. Nucleic Acids Res. 52, D134–D137 (2024).

133. Joseph7e. Assign-Taxonomy-with-BLAST, https://github.com/ Joseph7e/Assign-Taxonomy-with-BLAST (2020).

134. Lejzerowicz, F., Voltsky, I. & Pawlowski, J. Identifying active foraminifera in the Sea of Japan using metatranscriptomic approach. Deep Sea Res. Part II: Topical Stud. Oceanogr. 86, 214–220 (2013).

135. Li, Q., Lei, Y., Li, H. & Li, T. Distinct responses of abundant and rare foraminifera to environmental variables in the Antarctic region revealed by DNA metabarcoding. Front. Mar. Sci. 10, 1089482 (2023).

136. WoRMS. World Register of Marine Species, https://www. marinespecies.org (2025).

137. Callahan, B. J., McMurdie, P. J. & Holmes, S. P. Exact sequence variants should replace operational taxonomic units in marker-gene data analysis. ISME J. 11, 2639–2643 (2017).

138. Patterson, R. T. & Fishbein, E. Re-examination of the statistical methods used to determine the number of point counts needed for micropaleontological quantitative research. J. Paleontol. 63, 245–248 (1989).

139. Fatela, F. & Taborda, R. Confidence limits of species proportions in microfossil assemblages. Mar. Micropaleontol. 45, 169–174 (2002).

140. Walker, J. S. & Cahill, N. Influence of Foraminifera Count Size and Rare Species on Transfer Function Results Used in Sea-Level Reconstructions. J. Foraminifer. Res. 54, 107–116 (2024).

141. Sawai, Y. et al. Relationships between diatoms and tidal environments in Oregon and Washington, USA. Diatom Res. 31, 17–38 (2016).

142. Blott, S. J. & Pye, K. GRADISTAT: a grain size distribution and statistics package for the analysis of unconsolidated sediments. Earth Surf. Proc. Land 26, 1237–1248 (2001).

143. Callard, S. L., Gehrels, W. R., Morrison, B. V. & Grenfell, H. R. Suitability of salt-marsh foraminifera as proxy indicators of sea level in Tasmania. Mar. Micropaleontol. 79, 121–131 (2011).

144. Hawkes, A. D., Horton, B., Nelson, A. & Hill, D. The application of intertidal foraminifera to reconstruct coastal subsidence during the giant Cascadia earthquake of AD 1700 in Oregon, USA. Quatern Int 221, 116–140 (2010).
145. Kemp, A. C., Wright, A. J. & Cahill, N. Enough is enough, or more is more? Testing the influence of foraminiferal count size on reconstructions of paleo-marsh elevation. *J. Foraminifer. Res.* **50**, 266–278 (2020).

146. Saunders, J. B. Recent foraminifera of mangrove swamps and river estuaries and their fossil counterparts in Trinidad. *Micropaleontology* **4**, 79–92 (1958).

147. Blaauw, M. & Christen, J. A. Flexible paleoclimate age-depth models using an autoregressive gamma process. *Bayesian Anal.* **6**, 457–474 (2011).

148. Reimer, P. J. et al. The Intcal20 Northern Hemisphere radiocarbon age calibration curve (0-55 Cal Kbp). *Radiocarbon* **62**, 725–757 (2020).

149. Reimer, P. J., Brown, T. A. & Reimer, R. W. Discussion: reporting and calibration of post-bomb 14C data. *Radiocarbon* **46**, 1299–1304 (2004).

150. Hughes, J. B. & Hellmann, J. J. The application of rarefaction techniques to molecular inventories of microbial diversity. *Methods Enzymol.* **397**, 292–308 (2005).

151. Lemos, L. N., Fulthorpe, R. R., Triplett, E. W. & Roesch, L. F. Rethinking microbial diversity analysis in the high throughput sequencing era. *J. microbiological methods* **86**, 42–51 (2011).

152. Chao, A. Nonparametric estimation of the number of classes in a population. *Scand. J. Stat.* **11**, 265–270 (1984).

153. Troth, C. R., Sweet, M. J., Nightingale, J. & Burian, A. Seasonality, DNA degradation and spatial heterogeneity as drivers of eDNA detection dynamics. *Sci. Total Environ.* **768**, 144466 (2021).

154. Larson, M. G. Analysis of variance. *Circulation* **117**, 115–121 (2008).

155. Sakata, M. K. et al. Determining an effective sampling method for eDNA metabarcoding: a case study for fish biodiversity monitoring in a small, natural river. *Limnology* **22**, 221–235 (2021).

156. Anderson, M. J. Permutational multivariate analysis of variance (PERMANOVA). *Wiley statsref: statistics reference online*, 1-15 (2014).

157. Oksanen, J. et al. The vegan package. *Community Ecol. package* **10**, 719 (2007).

158. Ter Braak, C. J. Canonical correspondence analysis: a new eigenvector technique for multivariate direct gradient analysis. *Ecology* **67**, 1167–1179 (1986).

159. Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: and this is not optional. *Front Microbiol* **8**, 2224 (2017).

160. McMurdie, P. J. & Holmes, S. Waste not, want not: why rarefying microbiome data is inadmissible. *PLoS computational Biol.* **10**, e1003531 (2014).

161. Joseph, N., Paulson, C., Corrada Bravo, H. & Pop, M. Robust methods for differential abundance analysis in marker gene surveys. *Nat. Methods* **10**, 1200–1202 (2013).

162. Haslett, J. et al. Bayesian palaeoclimate reconstruction. *J. R. Stat. Soc. Ser. A: Stat. Soc.* **169**, 395–438 (2006).

163. Williams, S., Garrett, E., Moss, P., Bartlett, R. & Gehrels, W. R. Development of a training set of contemporary salt-marsh Foraminifera for Late Holocene sea-level reconstructions in southeastern Australia. *Open Quat.* **7**, 1–29 (2021).

164. Lands Department of Hong Kong. *Hong Kong Geographic Data*, https://www.landsd.gov.hk/en/resources/mapping-information/hk-geographic-data (2025).

165. Fu, S., Xiong, H., Zong, Y. & Huang, G. Reasons for the low sedimentation and slow progradation in the Pearl River delta, southern China, during the middle Holocene. *Mar. Geol.* **423**, 106133 (2020).

## Acknowledgements

This research was supported by the Research Grants Council of Hong Kong (GRF Project no.: 27300221 and 17303925). We would like to thank the

World Wide Fund for Nature for granting us permission to conduct research within the nature reserve. We appreciate WWF's efforts in maintaining natural habitats and conserving wildlife in Hong Kong. Our sincere thanks also go to Miss Qin Yonghui, Mr. Gao Chengcheng, and Miss Eunice Leung (names are listed in no particular order) for their tremendous assistance with fieldwork. We gratefully acknowledge Mr. Gu Yifei for his support in conducting data analysis. We also express our gratitude to Prof. Niamh Cahill for providing insights on Bayesian transfer function applications.

## Author contributions

Zhaojia Liu conducted sample collection, laboratory work, data analysis and wrote the original manuscript. Nicole S. Khan conceptualized this study, provided supervision, assisted with sample collection and laboratory work, and contributed to manuscript writing, reviewing and editing. Howard K.Y. Yu conducted sample collection and laboratory work, and reviewed the manuscript. Arthur Chung assisted with laboratory work and data analysis, and reviewed the manuscript. Magali Schweizer provided support for laboratory work and data analysis, and contributed to manuscript review and editing. Celia Schunter conceptualized the study, provided supervision and laboratory support, and contributed to manuscript reviewing and editing.

## Competing interests

The authors declare no competing interests. Nicole S. Khan is an Editorial Board Member for Communications Earth & Environment, but was not involved in the editorial review of, nor the decision to publish this article.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-025-03001-w.

**Correspondence** and requests for materials should be addressed to Nicole S. Khan.

**Peer review information** *Communications Earth & Environment* thanks Marco Coolen and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Alireza Bahadori. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025