Temporal and spatial variability in mudflat and mangrove foraminiferal eDNA communities in subtropical environments and their implication for sea-level reconstruction

Reconstructing relative sea level (RSL) is essential for understanding coastal evolution and mitigating impacts of climate change. Foraminiferal assemblages are established indicators for past sea levels, requiring modern training sets that link assemblage composition to known tidal elevations. However, seasonal and spatial variability in assemblages can affect reconstruction reliability. Foraminiferal environmental DNA (eDNA) has recently been validated as a robust proxy for reconstructing past sea levels. However, eDNA assemblages may be subject to similar spatiotemporal variability along with additional factors affecting DNA transport and preservation. The spatiotemporal stability of intertidal eDNA assemblages – and their influence on RSL reconstruction – remains uncertain. We conducted a 2-year eDNA monitoring study at three stations of varying tidal elevation, sampling across dry and wet seasons to assess variability in mangrove and mudflat environments of Hong Kong. Mid-mangrove eDNA communities exhibited temporal and spatial stability, whereas mudflat and upper-mangrove assemblages, particularly among monothalamous foraminifera, showed pronounced seasonal shifts. Despite this variability in the upper mangrove, eDNA-based elevation estimates in mangrove environments consistently aligned with observed elevations (within 95 % credible intervals), demonstrating reliable calibration for paleoenvironmental applications. However, mudflat samples exhibited an overprediction bias, especially during the wet season, reflecting their heightened sensitivity to exogenous eDNA inputs. These findings highlight the need to account for seasonal and environmental variability in eDNA-based RSL reconstruction. Stable mangroves are optimal for establishing elevation–assemblage relationships, while transitional or mudflat zones require caution due to higher variability. Our study provides guidance for applying foraminiferal eDNA as a proxy in past sea-level reconstruction in dynamic coastal settings.

1 Introduction

Coastal ecosystems face unprecedented pressures from climate change and sea-level rise, requiring robust approaches to monitor environmental change and guide management decisions. Foraminifera are valuable bioindicators in coastal environments because they are sensitive to environmental stressors, including salinity, pollution, temperature, nutrient availability, eutrophication, and climate-driven changes (Scott et al., 2007; Romano et al., 2021; Koukousioura et al., 2011; Youssef et al., 2021; Frontalini et al., 2009; Bouchet et al., 2021). They also contribute to essential ecosystem services through nutrient cycling and sediment processes (Langlet et al., 2020; Piña-Ochoa et al., 2010; Risgaard-Petersen et al., 2006), making their monitoring valuable for integrated coastal management by assessing both current ecosystem health and historical environmental change. Among their various applications, an important use of foraminifera is to reconstruct past relative sea level (RSL) changes, which is essential for understanding coastal evolution and predicting future changes (Khan et al., 2017; Tan et al., 2023; Xiong et al., 2018). Foraminifera are sensitive to the frequency and duration of tidal inundation, forming vertically zoned assemblages along elevation gradients in intertidal habitats (Berkeley et al., 2009; Scott and Medioli, 1978). Traditional RSL reconstructions rely on morphological identification of foraminiferal tests in sediment cores, calibrated with modern training sets that relate community composition to tidal elevation (Barnett et al., 2016; Rush et al., 2021).

A key challenge is the variability in surface foraminiferal assemblages, which fluctuate both temporally (seasonal and interannual variation) and spatially at sub-meter scales (Berkeley et al., 2008; Buragohain and Ghosh, 2021; Horton and Edwards, 2003; Buzas, 1970; Buzas et al., 2015; Kemp et al., 2011; Richirt et al., 2020). This variability can bias RSL reconstructions, especially because modern training sets are often based on single samples collected at one time point per site (Edwards and Horton, 2000; Kemp et al., 2013; Scott and Medioli, 1978; Yu et al., 2025; Horton and Edwards, 2006; Horton and Culver, 2008). To minimize this effect, RSL studies typically use dead foraminifera assemblages for calibration (Horton, 1999; Yu et al., 2025; Hawkes et al., 2010) because they are more stable over time and space (Hayward et al., 1996; Kemp et al., 2011; Murray and Alve, 2000; Scott et al., 2007; Walker et al., 2020). In contrast, live foraminiferal assemblages are more responsive to short-term environmental variability. Nevertheless, recent work shows that accounting for the temporal and spatial variability in dead assemblages can further improve RSL reconstruction accuracy (Walker et al., 2020).

Environmental DNA (eDNA) analysis has recently emerged as a promising alternative to conventional morphological approaches for studying foraminifera (Pawlowski et al., 2014). The application of eDNA techniques to foraminiferal community analysis enhances the ability to detect subtle ecological shifts in transitional ecosystems before visible habitat alterations occur (Singer et al., 2023). Furthermore, foraminiferal eDNA shows distinct vertical zonation correlated with tidal elevation, enabling the use of transfer functions for late Holocene RSL reconstruction (Liu et al., 2025). eDNA can detect taxa that do not readily fossilize and are often missed in traditional studies, providing a more comprehensive and high-resolution view of foraminiferal diversity (Lejzerowicz et al., 2013; Pawłowska et al., 2014; Singer et al., 2023) and potentially improving the predictive power of RSL reconstructions (Liu et al., 2025).

Despite these advantages, eDNA-derived assemblages may be subject to similar challenges to morphologically based approaches, as they can be influenced by spatial and temporal variability. eDNA comprises genetic material from multiple sources, including DNA originating locally from deceased populations (autochthonous) or transported from external sources (allochthonous) (Ellegaard et al., 2020; Nagler et al., 2022), as well as propagules (Singer et al., 2023; Brinkmann et al., 2023; Fouet et al., 2024) and living individuals (Angeles et al., 2020), all of which may vary seasonally. Furthermore, eDNA is especially sensitive to temporal or small-scale spatial variations (Matsuoka et al., 2021; Bista et al., 2017; Urabe et al., 2025; Singer et al., 2023), making it a powerful tool for biomonitoring (Brinkmann et al., 2023; He et al., 2019; Pawlowski et al., 2016) and biodiversity assessment (Bakker et al., 2019; Gu et al., 2023; Lejzerowicz et al., 2014). However, this sensitivity to temporal and small-scale spatial variations may also pose challenges for paleoenvironmental applications. The extent and implications of this variability for sea-level reconstructions remain poorly understood and have not been systematically assessed.

To address this knowledge gap, we investigated the temporal and spatial variability in foraminiferal eDNA assemblages and their effects on transfer-function-based RSL reconstructions in mangrove and mudflat environments in Hong Kong. Over 2 years, we conducted semi-annual monitoring to (1) assess seasonal dynamics and spatial variation in eDNA assemblages across tidal elevations, (2) identify environmental drivers of temporal variation, and (3) determine how seasonality influences the accuracy of eDNA-based RSL reconstructions. Our findings reveal clear seasonality in foraminiferal eDNA assemblages in a subtropical coastal wetland and provide recommendations for improving eDNA sampling strategies for RSL reconstruction.

2 Material and methods

2.1 Study area

Our study was conducted in the Mai Po Nature Reserve, situated on the eastern fringe of the Pearl River Delta (PRD), South China (Fig. 1b). The PRD is one of the world's largest and most dynamic river deltas, supporting high biodiversity within its extensive networks of rivers, estuaries, and intertidal wetlands (WWF, 2020). Rapid urbanization and industrial development, particularly in neighboring cities Shenzhen and Hong Kong, have influenced regional hydrology, sediment dynamics, and water quality throughout the Deep Bay and the Mai Po wetlands (Sun et al., 2017). The region experiences semidiurnal tides with a great diurnal amplitude of 1.96 m (from mean lower low water [MLLW] to mean higher high water [MHHW]), recorded by the nearest tide gauge station at Tsim Bei Tsui.

Figure 1 Study area and sampling locations. (a) Location of the study area in the Pearl River Delta. (b) The location of Deep Bay within the Pearl River Delta, showing climatic monitoring stations at Lau Fau Shan (LFS) and Hong Kong International Airport (HKA) (red squares). (c) Monitoring stations set up at interior mangrove forest (upper-mangrove station), mangrove forest (mid-mangrove station), and mudflat–mangrove transitional zone (upper-mudflat station) are shown, along with the nearest tide gauge at Tsim Bei Tsui. The position of the surface transect of the established modern training set in Liu et al. (2025) is also shown. Land surface elevation was derived from the lidar digital elevation model supplied by the Hong Kong SAR government for the Mai Po areas.

The climate is humid subtropical, characterized by distinct wet and dry seasons driven by seasonal monsoons (Yu et al., 2023; Ding and Chan, 2005). The southwest monsoon brings the wet season (April–September), which accounts for 77 % of annual precipitation, with maximum monthly totals exceeding 1000 mm. The northeast monsoon corresponds with the dry season (October–March), when the lowest monthly total precipitation drops below 1 mm (Hong Kong Observatory, 2024) (Fig. 2, Table S1 in the Supplement). The precipitation shift between dry and wet seasons typically lags about 1 month behind monsoon onset, with heavy rainfall beginning in late April to early May and ending from late October to early November (Hong Kong Observatory, 2024). The average monthly air temperature during the study period (from November 2022 to May 2024) ranged from 15–30 °C and was lowest in February (15 °C) and highest in July (30 °C). The seasonal contrast between wet and dry conditions affects porewater pH, salinity, and other environmental factors (Li et al., 2011; Wang et al., 2013; Chen et al., 2022), likely influencing foraminiferal community composition and eDNA detection.

Figure 2 Climatic conditions at the study area during the monitoring period (November 2022–May 2024), with the sampling period for modern training set samples from Liu et al. (2025) also indicated. Time series of monthly averaged air temperature (top panel), sea surface temperature (middle panel), and total precipitation (bottom panel) show characteristic dry (red) and wet (blue) seasonal climate patterns of the region. Monthly averaged air temperature and total precipitation data were obtained from the monitoring station at Lau Fau Shan (the nearest station to our study area), while sea surface temperature data were sourced from the monitoring station at Hong Kong International Airport, both maintained by the Hong Kong Observatory (Hong Kong Observatory, 2024).

Within Mai Po, a shore-normal transect was established spanning the mudflat–mangrove transitional zone to the interior mangrove forest, coinciding with a previous foraminiferal eDNA modern training set (Fig. 1c) (Liu et al., 2025). This transect captures the ecological gradient from barren mudflat to mature mangrove forest, enabling assessment of how environmental factors shape foraminiferal eDNA assemblages. Three monitoring stations were positioned along this transect in different environmental zones – (1) mudflat–mangrove transitional zone (upper-mudflat station), (2) mangrove forest (mid-mangrove station), and (3) interior mangrove forest (upper-mangrove station) (Fig. 1c) – to capture the temporal and small-scale variation in eDNA assemblages along the seaward–landward gradient. The flora in the mudflat–mangrove transitional zone is dominated by Sporobolus alterniflorus, with occasional non-native Sonneratia apetala (Yu et al., 2025), while the mangrove stands are primarily occupied by Kandelia candel, Aegiceras corniculatum, and Avicennia marina (Li et al., 2019; Lee, 2000).

2.2 Sampling strategy

Given the study area's distinct dry and wet seasonal pattern, we conducted semi-annual sampling during both seasons. Samples were collected at each monitoring station in late November (dry season) and late May (wet season), after seasonal climatic factors had stabilized for at least 2 weeks. For example, the dry-season sampling in November 2022 occurred during a period of consistently low precipitation (daily total precipitation 0–9.6 mm), while wet-season sampling in May 2023 followed the onset of increased rainfall that characterized the transition from the dry winter months (Hong Kong Observatory, 2024). This timing ensured that the foraminiferal eDNA assemblages captured reflected true seasonal variation related to the respective collection periods. Within a 1 × 1 m plot at each monitoring station, three replicate sediment samples (approximately 20 cm³ each) for eDNA analysis were collected randomly from the top 1 cm of the surface (Liu et al., 2025). This replication was designed to assess microscale spatial variability (patchiness) in foraminiferal eDNA assemblage composition within each monitoring station (Lejzerowicz et al., 2014; Walker et al., 2020). The position of each collected sample was recorded within 10 × 10 cm cells using a quadrat (Smith et al., 2021) to avoid repeatedly sampling from the same location across seasons. The sampling procedure followed Liu et al. (2025). Briefly, all sampling equipment was sanitized with 10 % bleach between each collection to prevent cross-contamination. All eDNA samples were kept on ice during field collection and promptly transferred to a −20 °C freezer upon return to the laboratory on the sampling day.

To examine the potential environmental drivers of temporal variation in foraminiferal eDNA, an additional 20 cm³ of sediment was sampled simultaneously at each station to measure environmental variables, including porewater salinity and pH, and stable carbon isotope geochemistry (total organic carbon (TOC), total nitrogen (TN) and δ¹³C). These samples were also kept on ice once collected and stored in a 4 °C refrigerator in the lab.

The elevation of each monitoring station was measured following the procedure introduced in Liu et al. (2025), relative to the Hong Kong Principal Datum (PD). Specifically, we took three elevation measurements at each monitoring station due to the slightly uneven surface topography. Because of the flat topography of the mangrove forest in the study area, the elevations of the upper- and mid-mangrove stations fell within similar ranges (2.17–2.22 mPD). We converted the elevation to standardized water level index (SWLI) units (Kemp and Telford, 2015), where a value of 100 corresponds to mean tide level (MTL) and 200 to mean higher high water (MHHW), based on tidal data obtained from the Tsim Bei Tsui tide gauge (Fig. 1c).

2.3 eDNA analysis

2.3.1 eDNA extraction, amplification, and sequencing

We extracted eDNA from sediment samples using the 0.5 g PowerSoil^® DNA isolation kit (QIAGEN) according to manufacturer instructions. A two-step polymerase chain reaction (PCR) amplification protocol was conducted to prepare sequencing libraries. To minimize PCR bias and to cover total genetic and taxonomic diversity (Nichols et al., 2018; Singer et al., 2023), DNA extraction was performed in two replicates for each sample, followed by duplicate PCR amplifications of each extraction replicate. The PCR amplification procedure follows Liu et al. (2025). Briefly, the first PCR was conducted with primers s14F1–s15 to amplify the 135–190 bp fragment within the hypervariable region 37f of the foraminiferal SSU rDNA gene (Lejzerowicz et al., 2014). Both primers included overhang adapter sequences compatible with Illumina indices following the Illumina MiSeq system sequencing protocol. Replicate PCR products from each sample were pooled and purified using AMPpure XP beads (Beckman Coulter, Singapore) following the Illumina MiSeq system sequencing protocol. Subsequently, an indexing PCR was performed using the Nextera XT Index Kit (Illumina) to append dual indices to the purified amplicons. Indexed PCR products were purified again using AMPpure XP beads. Purified samples were pooled to a final concentration of 20 µM in 10 mM Tris (pH 8.5) and sent to Novogene Co., Ltd for sequencing on the NovaSeq platform using paired-end 150 bp reads (PE 150).

2.3.2 Bioinformatics

Bioinformatics processing followed the pipeline introduced in Liu et al. (2025). Briefly, raw sequence data were demultiplexed and primers removed by Cutadapt v2.8 (Martin, 2011); the produced pair-end reads were further merged and trimmed by PEAR v0.9.11 with a quality threshold of 26 (Zhang et al., 2014). Concatenated sequences were denoised and chimeras removed using VSEARCH (Rognes et al., 2016). Sequences were clustered into operational taxonomic units (OTUs) at 97 % similarity using VSEARCH. Quality control was applied to discard OTUs that were (1) either shorter than 90 bp or longer than 230 bp, (3) had fewer than 10 reads, or (3) appeared in negative controls with read counts equal to or exceeding those in any sample. Taxonomy was assigned using BLASTn v2.12.0 (Altschul et al., 1990) with an e value of 0.001 at species, family, and phylum levels based on identity thresholds of 98 %, 90 %, and 80 %, respectively, referring to the GenBank nucleotide database and an in-house foraminiferal DNA database (Table S2). OTUs with identities between 80 % and 90 % (assigned to phylum level) and OTUs with unknown taxonomy were classified as “undetermined OTUs” and treated individually as taxa in subsequent analyses. OTUs with identity thresholds below 80 % were not considered foraminifera and were removed.

2.4 Environmental variables

To assess the influence of environmental variables on foraminiferal eDNA assemblage composition over time, porewater salinity, pH, sediment total organic carbon (TOC), total nitrogen (TN), and δ¹³C were measured following the procedures introduced in Yu et al. (2025) and Liu et al. (2025). Briefly, salinity and pH were measured using a Thermo Scientific^™ A3255 pH/conductivity multimeter from porewater obtained by centrifuging sediment samples (Sawai et al., 2016). Sediment samples for δ¹³C, TOC, and $\mathrm{C}/\mathrm{N}$ analyses were pretreated with 5 % hydrochloric acid (HCl) prior to analysis on a EuroVector EA3028 Nu Horizon isotope-ratio mass spectrometer (IRMS) and an EA Isolink elemental analyzer. Additionally, climate data (air temperature, sea surface temperature, and precipitation) during our study period (November 2022–May 2024) were obtained from the monitoring station at Lau Fau Shan (nearest to our study area), with additional sea surface temperature measurements obtained from the Hong Kong International Airport station maintained by the Hong Kong Observatory (Hong Kong Observatory, 2024) to evaluate seasonal climatic influences on environmental variables and foraminiferal eDNA assemblages (Table S1).

2.5 Statistical analysis

To assess seasonal variation in the abundance of dominant foraminiferal taxa at each sampling station, we performed a one-way analysis of variance (ANOVA) with 1000 permutations (p-value threshold = 0.05) (Troth et al., 2021; Chung et al., 2024) (Table S3). Identical seasons from years 1 and 2 were merged as a single factor in this analysis. Dominant taxa were defined as those with > 5 % relative abundance in at least one sample. Prior to ANOVA, we assessed homogeneity of variances for each taxon's abundance at each station and applied log transformation as needed to meet ANOVA assumptions. The relative abundances of dominant taxa were plotted, and a heatmap was generated using merged replicate samples for each station and season to visualize patterns in taxonomic composition and diversity across environments and seasons. The heatmap was produced using the ggheatmap function from the heatmaply package (Galili et al., 2018) in R (v4.4.2).

To examine overall patterns of temporal (seasonal) and spatial variation in foraminiferal eDNA assemblages, we conducted non-metric multidimensional scaling (NMDS) based on Bray–Curtis dissimilarity of taxa relative abundance data, using the vegan package (v2.6-4) in R (Oksanen et al., 2007). The NMDS ordination was constructed using dominant taxa (> 5 % relative abundance in at least one sample), and the first two NMDS dimensions were used to visualize compositional dissimilarities among stations and seasons. To test the significance of observed patterns, we used permutational multivariate analysis of variance (PERMANOVA) with 999 permutations (Anderson, 2014), also based on Bray–Curtis distances (Table S4). PERMANOVA was used to evaluate whether eDNA assemblages differed significantly between seasons (with identical seasons from years 1 and 2 merged as a single factor) and among stations (p-value threshold = 0.05).

To evaluate microscale spatial variability (patchiness) within each monitoring station, we calculated pairwise Bray–Curtis dissimilarity among the three replicate samples collected within each 1 × 1 m plot at every station and sampling event. The magnitude of within-station dissimilarity (differences between replicate samples collected simultaneously from the same monitoring station representing microscale spatial variability) was then compared to among-station dissimilarity (differences among samples collected during the same season from different monitoring stations) using a permutation test (10 000 permutations, p-value threshold = 0.05) (Welch, 1990). This analysis tested whether microscale patchiness within stations was significant relative to broader environmental differences. Additionally, ANOVA (1000 permutations, p-value threshold = 0.05) was used to assess whether the dispersion (i.e., within-station variability) of eDNA assemblages differed significantly among the three monitoring stations during the same sampling period.

To determine the influence of environmental variables on the composition of foraminiferal eDNA assemblages, we performed redundancy analysis (RDA) and partial RDA (pRDA) using Hellinger-transformed relative abundance data (Table S5), implemented in the vegan package (v2.6-4) (Legendre and Gallagher, 2001) in R. Only dominant taxa (> 5 % relative abundance in at least one sample) were included in the RDA. The environmental variables considered included salinity, pH, TOC, carbon–nitrogen ratio ($\mathrm{C}/\mathrm{N}$), and δ¹³C. The significance (p<0.05) of environmental variables and the first two ordination axes was determined using ANOVA with 1000 permutations. To examine whether the control of environmental variables differs between seasons, separate RDA analyses were also performed for the wet- and dry-season samples. pRDA was conducted with tidal elevation (SWLI) of each monitoring station included to further assess the relative contributions of tidal elevation and the tested environmental variables to variance in the foraminiferal eDNA assemblage. Pairwise associations between foraminiferal taxa and environmental variables (salinity, pH, TOC, TN, $\mathrm{C}/\mathrm{N}$ ratio) were further evaluated using Spearman's rank correlation coefficient implemented in R (Spearman, 2010). Additionally, linear mixed-effects models (LMMs) were run to examine the significance of variation (p-value threshold = 0.05) in the environmental variables across seasons at each monitoring station (Luke, 2017). To further quantify the proportion of variance in assemblage composition explained by physicochemical (tidal elevation and environmental variables) and climatic (temperature, sea surface temperature, and precipitation) variables, we conducted variance partitioning using the varpart function in vegan (Peres-Neto et al., 2006).

To investigate how seasonal variations in eDNA assemblages at each monitoring station potentially influenced sea-level reconstruction, we applied a Bayesian transfer function (BTF) to estimate the elevation for each season at each station (Cahill et al., 2016; Liu et al., 2025; Yu et al., 2025). The BTF was constructed using a modern training set from the same study site (Liu et al., 2025), comprising 59 surface samples and spanning the elevation range from the lowest to the highest astronomical tide. The BTF was constructed following the same procedures described in Liu et al. (2025). Briefly, eDNA sequence counts were rarefied to the lowest sequencing depth in the dataset – after excluding one modern sample with low read numbers – to account for differences in sequencing depth among samples. To minimize bias from low-abundance taxa, only those exceeding 5 % abundance in at least one sample were included (Patterson and Fishbein, 1989; Fatela and Taborda, 2002; Liu et al., 2025), consistent with the criteria for taxa applied in RDA. Although OTUs assigned to the planktonic family Globigerinitidae exceeded this threshold, they were excluded from the BTF to reduce uncertainty, as they are not considered in situ taxa. The distribution of estimated elevations (SWLI) for each station and season was compared to observed elevations, and uncertainties were reported as the maximum value of the 2σ interval across replicates.

3 Results

3.1 Seasonal variation in dominant taxa

Several dominant foraminiferal taxa (> 5 % relative abundance in at least one sample) exhibited clear and, in some cases, statistically significant seasonal patterns in abundance across the three sampling stations (Figs. 3 and 4, Table S2). For example, Saccamminidae were significantly more abundant during the dry season at the upper-mangrove and upper-mudflat stations (ANOVA, p<0.05), while Ammoniidae showed significantly higher abundance in the wet season at the mid-mangrove station. Detailed relative abundances for all taxa at each station are provided in Table S6.

Figure 3 Relative abundance (%) of dominant foraminiferal taxa (> 5 % relative abundance in at least one sample) at three sampling stations during the wet (blue) and dry (red) seasons. Taxa that exhibited a statistically significant seasonal pattern at least at one of the monitoring stations determined by one-way ANOVA (Table S3) are marked with an asterisk.

3.2 Community structure across stations and seasons

The seasonal patterns in dominant taxa corresponded to distinct shifts in overall community composition across monitoring stations. During the dry season, the upper-mangrove, mid-mangrove, and upper-mudflat stations were all dominated by Saccamminidae, while Miliamminidae were also abundant in the mid-mangrove (Fig. 4). In the wet season, the upper mangrove was dominated by OTU23 and OTU59, while the mid-mangrove was characterized by high abundances of Ammoniidae and OTU3. During the wet season, the upper mudflat exhibited decreased abundance of Saccamminidae but increased abundances of Miliamminidae, Ammoniidae, and OTU3 compared to the dry season (Fig. 4).

Figure 4 Heatmap of foraminiferal eDNA assemblages showing the relative abundance (%) of dominant foraminiferal taxa at each sampling station across wet and dry seasons. Replicate samples from each station were pooled by season to represent seasonal shifts in overall community structure.

Both environment and season influenced community composition (Fig. 5a). Samples from the upper mangrove and mid-mangrove clustered together, indicating that their foraminiferal community compositions were more similar to each other than to the upper-mudflat samples, although there was still a significant difference between these two stations (PERMANOVA; F= 2.21, R²= 0.09, p= 0.014). In contrast, samples from the upper mudflat formed a more significantly separated cluster from the upper mangrove (PERMANOVA; F= 5.18, R²= 0.19, p= 0.001) and mid-mangrove (PERMANOVA; F= 5.65, R²= 0.20, p= 0.001), reflecting a distinct community composition (Fig. 5a). In the upper mudflat, dry- and wet-season samples were clearly separated in ordination space, also showing a significant difference between seasons (PERMANOVA; F= 4.42, R²= 0.31, p= 0.005; Table S4). Within the upper mangrove, dry- and wet-season samples largely overlapped in ordination space, yet PERMANOVA indicated a significant seasonal difference in community composition (F= 2.34, R²= 0.19, p= 0.01). In contrast, no significant seasonal difference was detected in the mid-mangrove samples (PERMANOVA; F= 1.93, R²= 0.16, p = 0.079). In addition, most taxa were centrally located in the NMDS ordination, suggesting they were common to both the upper-mangrove and mid-mangrove environments. However, monothalamids such as Saccamminidae, Vanhoeffenella sp., and calcareous-walled Nummulitidae were more closely associated with the upper-mudflat environment in the NMDS analysis (Fig. 5a).

3.3 Microscale (within-station) spatial variability

Samples showed significantly higher among-station heterogeneity than within-station (microscale) variability (Fig. 5b). Although the upper-mangrove station sampled in November 2022 exhibited the highest within-station dissimilarity (Table 1), within-station variability was significantly lower than among-station differences (permutation test, p<0.01). Furthermore, overall variation in within-station dispersion across all monitoring stations was not statistically significant (ANOVA, p>0.05), indicating limited patchiness at the scale sampled.

Figure 5 Dissimilarity of foraminiferal eDNA assemblages in each sample. (a) Non-metric multi-dimensional scaling (NMDS) ordination of foraminiferal eDNA assemblages. NMDS ordination based on Bray–Curtis dissimilarities of foraminiferal eDNA assemblages at each sampling station during dry and wet seasons (stress value = 0.17). The NMDS was constructed using dominant taxa with > 5 % relative abundance in at least one sample. Samples from different environments and seasons are indicated by distinct colors and shapes. (b) Box plot comparing average Bray–Curtis dissimilarity values among stations (1: among-station variation) and within stations (2: within-station variation).

Table 1 The average within-station Bray–Curtis (BC) dissimilarity of each monitoring station per collection.

3.4 Influence of environmental variables

Both static (elevation) and seasonally variable environmental factors (salinity, pH, TOC, TN, and δ¹³C) were measured to assess their influence on foraminiferal eDNA assemblages. The upper- (2.17–2.18 m PD, 189–190 SWLI) and mid-mangrove (2.17–2.22 m PD, 190–195 SWLI) monitoring stations shared similar elevations due to the flat topography of the mangrove forest compared to the upper-mudflat station (1.27–1.31 m PD, 93–99 SWLI). The upper-mangrove and mid-mangrove environments are characterized by lower salinity and higher pH, while the upper mudflat exhibits lower total organic carbon (TOC) and higher δ¹³C values (Fig. 6c, Table 2).

Figure 6 Environmental drivers of foraminiferal community structure. (a) Redundancy analysis (RDA) bioplot of dominant foraminiferal taxa (> 5 % relative abundance in at least one sample) and examined environmental variables. (b) Variance partitioning of community composition among physicochemical environmental variables and climatic data. Residual values are also shown. (c) Vertical bars represent the observed minimum–maximum range of environmental variables measured at three monitoring stations during dry and wet seasons. (1: upper-mangrove station; 2: mid-mangrove station; 3: upper-mudflat station).

Several measured environmental variables exhibited seasonal variation across the three environments (Fig. 6c, Table 2). Nevertheless, results from the LMMs indicated that only salinity at the upper-mudflat station exhibited a significant difference, being much higher in the dry season (18.0–20.2) than in the wet season (8.1–10.9, LMMs, p<0.05). Taxa distribution was influenced by examined environmental conditions and climatic factors (temperature, precipitation, sea surface temperature) (Fig. 6). Salinity, pH, and TOC (ANOVA, p<0.01) significantly influenced foraminiferal eDNA composition across environments and seasons. Nevertheless, tidal elevation of each monitoring station remained the dominant driver shown by pRDA, explaining 13 % of the variance in eDNA assemblages, followed by salinity (11 %) and pH (7 %). Most taxa clustered near the center of the RDA bioplot, suggesting broad environmental tolerance, but some (e.g., Nummulitidae; Spearman's ρ with salinity = 0.20, ρ with pH = 0.19) exhibited positive correlations with salinity and pH (Fig. 6a). Many undetermined OTUs (e.g., OTU59; Spearman's ρ with TOC = −0.61) showed negative correlations with TOC (Fig. 6a). RDA analyses performed separately for the dry and wet seasons indicated that salinity and pH were the only significant drivers of community variation in both seasons (ANOVA, p<0.01). However, their relationship shifted seasonally, showing a positive association in the wet season but an inverse relationship in the dry season (Fig. S1). Variance partitioning analysis revealed that physicochemical environmental properties explained 36.2 % of the variance in eDNA assemblages, while climatic variables (including temperature, sea surface temperature, and precipitation) accounted for 5.5 % (Fig. 6b).

Table 2 Environmental variables measured at each monitoring station per season.

3.5 Relative sea level (RSL) estimation of each monitoring station

A foraminiferal eDNA BTF developed from a previously established modern training set (Liu et al., 2025) was used with the present eDNA dataset to estimate elevation for all replicate samples collected from each monitoring station across seasons and years. For the mangrove stations, the estimated elevations fell within the observed elevation range (within the 95 % uncertainty interval) for both dry and wet seasons (Fig. 7). In contrast, the BTF consistently overpredicted elevations in both seasons at the upper-mudflat station. The dry-season samples from the upper-mudflat station exhibited the narrowest prediction range (121–132 SWLI; Table S7), while the wet-season samples from the upper-mudflat station had the lowest average 1σ uncertainty (7.5 SWLI, 0.38 m) (Table S7).

Figure 7 Elevation estimates (standardized water level index unit, SWLI) for each monitoring station during dry and wet seasons based on foraminiferal eDNA assemblage. Distributions represent six replicate samples per station and season. Error bars indicate the maximum 2σ uncertainty across replicates. The range of observed elevation at each station is shown by the dashed gray bar.

4 Discussion

4.1 Temporal and spatial variation in foraminiferal eDNA assemblage

Foraminiferal eDNA assemblages varied temporally across seasons within specific habitats but were spatially stable at small scales (i.e., low patchiness) across our three intertidal monitoring stations. Temporal variation patterns likely reflected fundamental differences in the composition and sources of eDNA across environments. At the mid-mangrove station, community composition remained relatively stable throughout the year (Figs. 4, 5a). This stability suggests an environment which may buffer short-term seasonal changes due to reduced variability in environmental conditions, as observed in previous mangrove studies (Camp et al., 2016; Alongi, 2022). Consequently, the eDNA assemblage tends to reflect long-term environmental gradients, such as tidal elevation, which are shaped by persistent factors like the duration and frequency of inundation or subaerial exposure (Horton and Culver, 2008; Woodroffe et al., 2005). Additionally, the eDNA pool could be dominated by accumulated extracellular DNA from local populations across multiple generations, which can contribute > 40 % of total sediment eDNA (Caro et al., 2023; Carini et al., 2016; Ascher et al., 2009). This interpretation – that extracellular DNA may potentially dominate the assemblage – is consistent with traditional morphological studies, where dead foraminiferal tests also show limited seasonal variability compared to living (rose Bengal-stained) assemblages (Buzas et al., 2015; Hippensteel et al., 2002; Horton and Edwards, 2003; Horton and Murray, 2006; Walker et al., 2020). However, validating this hypothesis will require further research employing methods such as environmental RNA (eRNA) (Pearman et al., 2022; Pochon et al., 2017) or eDNA fractionation (Nagler et al., 2018, 2021).

In contrast to the mid-mangrove station stability, the upper-mangrove and upper-mudflat stations exhibited marked seasonal changes in the relative abundance of dominant taxa such as Saccamminidae and OTU3 (Figs. 3, 4). This indicates a greater response of the foraminiferal community to seasonal drivers in these environments, suggesting that the eDNA assemblages there could be more reflective of short-term environmental fluctuations (Goldberg et al., 2011; Ellegaard et al., 2020; Singer et al., 2023; Torti et al., 2018). Ecologically, the observed seasonal variation in eDNA assemblages among stations appears to be driven by the life histories and environmental sensitivities of key foraminiferal taxa. Hard-shelled foraminifera, such as those in the family Ammoniidae, exhibited significant seasonal variation (Table S3), with higher abundance in the wet season at the mid-mangrove station (Fig. 3). This is consistent with morphological surveys, which report that calcareous taxa typically reach peak densities during spring and summer, coinciding with warmer temperatures and increased phytoplankton and zooplankton productivity (Alve and Murray, 1999; Horton and Edwards, 2003; Berkeley et al., 2008; Murray and Alve, 2000; Richirt et al., 2020). These seasonal population cycles likely lead to increased eDNA production rates during reproductive periods (Buxton et al., 2017; Spear et al., 2015; Fukumoto et al., 2015), resulting in changes in eDNA assemblages that respond more rapidly than morphological methods and influence the taxonomic composition detected in our samples.

Conversely, monothalamous taxa – particularly Saccamminidae – were especially abundant during the dry season and contributed disproportionately to the observed seasonality. These taxa are known to be more sensitive to environmental change (He et al., 2019), and their seasonal dynamics are not captured in traditional morphological surveys using routine methodology due to poor preservation of their delicate tests (Schönfeld et al., 2012). Recent studies suggest that Saccamminidae may decline during summer, possibly due to colonization of their preferred substrate by algal blooms (Henderson, 2023), while eDNA surveys in French estuarine mudflats have reported their highest abundances in the fall (Singer et al., 2023). This pronounced seasonality in Saccamminidae abundance helps explain our earlier RSL reconstruction efforts at this site (Liu et al., 2025), where the modern training set – based on wet-season samples (primarily collected between June and August 2022) – showed comparatively lower Saccamminidae representation. In that study, Saccamminidae abundances reached up to 23 % (from MTL to MHHW), similar to the 9 %–36 % observed in the wet-season samples of the present study (Fig. 3), but much lower than the abundances observed during the dry season (4 %–60 %; Fig. 3). Although including monothalamids improves BTF performance (Liu et al., 2025), our results suggest that the timing of the collection of training set samples is an important consideration. In our study area, sample collection during the wet season – when Saccamminidae are less dominant – yields more balanced representation of indicator taxa and prevents compositional bias in eDNA data (Gloor et al., 2017). This approach is preferable for new studies aiming to explore the relationship between eDNA assemblages and tidal elevation, as it reduces the risk of dominant taxa masking underlying ecological patterns.

Our study also provides insights into the spatial variability in eDNA assemblages. In RSL reconstruction, small-scale spatial variability in foraminiferal assemblages may introduce bias into modern training sets, reducing their accuracy in capturing tidal elevation signals (Kemp et al., 2011; Yu et al., 2025). The upper-mangrove station exhibited the greatest heterogeneity among replicate samples (Table 1), but overall community structure and composition of dominant taxa remained consistent across all monitoring stations, indicating strong spatial stability at microscales (within 1 × 1 m). Although small-scale spatial variability in intertidal foraminiferal communities can be driven by the adaptation of living populations to local microhabitats and biotic factors such as predation pressure (Buzas, 1970, 1982) or food resource availability (Alve and Murray, 2001; Fontanier et al., 2003), the environments and conditions at our study site seemed to be relatively homogeneous. Consequently, the overall eDNA assemblages were primarily shaped by larger-scale environmental differences among the different stations (Fig. 6c), including salinity, δ¹³C, TOC, and $\mathrm{C}/\mathrm{N}$ ratio, rather than by fine-scale spatial heterogeneity. Accordingly, they did not display the pronounced spatial patchiness observed in deep-sea (Lejzerowicz et al., 2014) or temperate mudflat environments (Singer et al., 2023). This spatial stability has important implications for both paleoenvironmental reconstruction and contemporary monitoring. It supports the sampling design used in modern training sets for tidal elevation reconstruction, where single samples from narrow elevation intervals (typically < 10 cm) are assumed to be representative of community structure at that elevation. It also suggests that eDNA-based monitoring in mangrove and mudflat systems can reliably capture the dominant community structure without intensive spatial replication under relatively uniform environmental conditions. This supports the use of single or pooled samples to represent local assemblages in coastal monitoring programs.

4.2 Factors influencing seasonality in foraminiferal eDNA assemblages

A range of abiotic and biotic environmental factors have long been recognized as important drivers of foraminiferal community composition in intertidal settings (Debenay et al., 2006; Armynot Du Châtelet et al., 2009b; Lei et al., 2017; Murray, 2006; Armynot Du Châtelet et al., 2009a). Our results further demonstrate that temporal fluctuations in environmental variables, shaped by seasonal and tidal cycles, are closely linked to the observed seasonality in foraminiferal eDNA assemblages. Salinity, pH, and total organic carbon (TOC) are significantly correlated with seasonal variation in foraminiferal eDNA community structure shown by RDA. These variables, in turn, are influenced by broader climatic drivers such as temperature, precipitation, and river discharge. Understanding these relationships is critical for the development and application of foraminiferal eDNA as an ecological indicator, as it highlights which environmental variables must be monitored or controlled to ensure reliable interpretation.

Among the short-term environmental variables examined here (Table 2), salinity was the primary driver of seasonal changes in foraminiferal eDNA assemblages (Fig. 6a). During the wet season, increased precipitation and river discharge lowered salinity, favoring taxa adapted to reduced salinity (e.g., Trochammina inflata (Montagu, 1808), which increased from 1 %–2 % to 1 %–18 % at the upper-mangrove station; Fig. 4) while constraining more stenohaline taxa (e.g., Saccaminidae sp. A279, which decreased from 1 %–17 % to < 2 % at the upper-mudflat station; Figs. 3, 4) that are intolerant of broader salinity ranges (Singer et al., 2023; Jorissen et al., 2022; Li et al., 2023). In parallel, lower pH during the wet season likely influenced the abundance of calcareous foraminifera by affecting calcification and metabolic processes (Dong et al., 2020; Le Cadre et al., 2003). Likewise, fluctuations in TOC content indicated changes in food availability important for juvenile recruitment (Murray, 2006), while fluctuation in $\mathrm{C}/\mathrm{N}$ ratio and δ¹³C is an indication of potential changes in the food resource composition (Bouillon et al., 2002). Collectively, these examined physicochemical environmental variables shaped the distinct seasonal community patterns observed across the intertidal zone (Figs. 3, 4). Although grain size variation was relatively minor within our study area (Liu et al., 2025; Yu et al., 2025) and was therefore not examined, sediment texture and local hydrodynamics may also influence foraminiferal eDNA distributions and should be considered in future studies spanning broader spatial scales.

In addition to environmental factors like salinity and pH, climatic variables such as temperature also influence foraminiferal assemblages in intertidal zones, as shown by both morphological (Kaminski et al., 2021; Lei et al., 2017; Murray and Alve, 2000) and eDNA-based studies (Singer et al., 2023). However, variance partitioning in our study revealed that physicochemical environmental variables collectively explained 36.2 % of variation in eDNA assemblages, while climatic factors collectively explained only 5.5 % of the variation in eDNA assemblages (Fig. 6b). Although increasing temperature can enhance organisms' metabolism and, consequently, foraminiferal growth (Dong et al., 2019; Prazeres and Pandolfi, 2016; Lombard et al., 2009), the seasonal temperature variation in our study area per sampling was relatively small (up to 4.6 °C; Fig. 2, Table S1). Therefore, it is likely that climatic factors primarily exert indirect effects on foraminiferal communities by modulating local short-term environmental conditions. The seasonal shift in the relationship between salinity and pH provides an example of this interaction (Fig. S1). During the wet season, salinity and pH showed a positive relationship, likely driven by freshwater input that simultaneously reduced both variables across the intertidal gradient. In contrast, during the dry season, salinity and pH were negatively correlated, indicating that site-specific processes may exert stronger control over climatic influences.

Environmental variables not only affect living foraminiferal communities but may also influence the properties and preservation of eDNA itself. Higher temperatures during the wet season may accelerate extracellular and dead organisms' eDNA degradation via increased physical stress and bacterial nuclease activity (Eichmiller et al., 2016; Lance et al., 2017; Strickler et al., 2015; Wood et al., 2020), potentially increasing the proportion of eDNA derived from living populations. On the other hand, increased river discharge and stronger hydrodynamics during the wet season may elevate transportation of exogenous eDNA (e.g., propagules) into the intertidal zone (Harrison et al., 2019; Ji et al., 2020; Nevers et al., 2020), potentially introducing uncertainty when correlating tidal elevation with vertical zonation of foraminiferal eDNA assemblages.

Meanwhile, the taphonomy of eDNA itself may partly shape the seasonal patterns observed in the foraminiferal eDNA assemblages. Although the upper-mangrove and upper-mudflat stations were enriched in Saccamminidae signals, these signals declined during the wet season (Figs. 3, 4). Because DNA can degrade rapidly in intertidal sediments (Schweizer, 2015), eDNA amplified from surface sediment is likely to be biased toward fresh intracellular DNA from living or recently active individuals (Angeles et al., 2020). When Saccamminidae no longer thrive during the wet season, their relative contribution to the total eDNA pool at the monitoring stations consequently decreases. Over longer timescales, however, eDNA preserved in the sediment archive is likely to represent a more integrated signal, derived mainly from extracellular or relict DNA and reflecting cumulative population flux rather than short-term living assemblages alone (Liu et al., 2025; Angeles et al., 2020).

Despite these seasonal effects, tidal elevation explained the largest proportion of variation in foraminiferal eDNA assemblages, highlighting its overriding and persistent influence that is consistent with previous morphological studies (Armynot Du Châtelet et al., 2018; Lei et al., 2017). Tidal elevation shapes distinct physicochemical gradients across the study area, as reflected in δ¹³C values and TOC content (Table 2), which both directly influence foraminiferal communities and mediate the effects of seasonal environmental changes at different elevations. While seasonal variations in salinity, pH, and TOC shape temporal patterns in foraminiferal assemblages, tidal elevation remains the dominant factor, emphasizing the key role of local environmental context in structuring intertidal communities.

4.3 Implications of temporal variation in foraminiferal eDNA for sea-level reconstruction

Our results demonstrate that foraminiferal eDNA assemblages serve as reliable proxies for tidal elevation in mangrove environments, with elevation estimates derived from eDNA BTFs for both dry and wet seasons generally falling within the observed elevation ranges of the sampled mangrove stations (Fig. 7). This reliability supports the use of mangrove-derived eDNA samples for constructing modern training sets throughout the year, although broader seasonal and interannual replication is required to confirm the stability of this relationship. In contrast, the accuracy of eDNA-based reconstructions was reduced at the upper-mudflat station situated at a transitional zone, where elevation estimates were consistently overpredicted (Fig. 7). This limitation is potentially driven by increased contributions of allochthonous and propagule-derived eDNA in these dynamic environments (Angeles et al., 2020; Goldstein and Alve, 2011; Singer et al., 2023), especially during the wet season. The potential influx of allochthonous eDNA, particularly from higher-elevation mangrove environments, expands the apparent environmental tolerances of some taxa (e.g., Vanhoeffenella sp. isolate 5180) and increases uncertainty in elevation estimates of the transitional zone. Additionally, anthropogenic disturbance in the transitional zone (Yu et al., 2025) may exacerbate habitat instability and further confound the eDNA signal. These findings mirror morphological studies, which report greater community turnover and weaker tidal elevation gradients at lower-elevation mudflats (Lei et al., 2017), resulting in larger offsets in RSL estimates.

Our findings support several recommendations for optimizing foraminiferal eDNA as a sea-level indicator for reconstructing RSL change: (1) prioritizing stable mangrove environments for constructing modern training sets, as these yield signals most closely linked to local tidal elevation; (2) avoiding environmental transitional zones, such as the upper mudflat, where possible to minimize allochthonous influences; and (3) utilizing single samples per homogeneous environment, as this is sufficient to capture local foraminiferal eDNA diversity in line with patterns observed in morphological studies (Horton and Edwards, 2006; Kemp et al., 2009; Kemp et al., 2012), which our results show are sufficient for spatial representation. By refining both sampling strategies and analytical approaches, foraminiferal eDNA has considerable potential to contribute to more accurate and robust reconstructions of past RSL change in coastal wetlands.

5 Conclusion

Our study provides new insights into the temporal and spatial variation in foraminiferal eDNA assemblages and their implications for RSL reconstruction. Foraminiferal eDNA assemblages at the mid-mangrove monitoring station remained relatively stable across seasons, demonstrating a possible buffering of seasonal variations by the environment or a potential extracellular eDNA dominance, whereas the upper-mangrove station showed significant seasonal variation. The mudflat environment displayed distinctly different foraminiferal eDNA assemblages between the dry and wet seasons, showing a potentially greater contribution from propagules and exogenous eDNA than in mangrove sites that are further from the open sea. Despite the influence of seasonally variable environmental and climatic factors, tidal elevation remained the dominant factor shaping eDNA assemblages. Elevation estimates produced by the eDNA BTF for the mangrove monitoring stations are accurate for both dry and wet seasons, while predictions for the upper-mudflat station during the wet season failed to be accurate. Our results reveal that foraminiferal eDNA assemblages display environment-specific seasonal variability, necessitating cautious interpretation of mudflat samples given their pronounced sensitivity to seasonal changes. Nevertheless, we demonstrate that eDNA remains a robust and broadly applicable proxy for RSL reconstruction, ensuring its utility for future studies.