Late Holocene environmental change in Celestun Lagoon, Yucatan, Mexico

Epikarst estuary response to hydroclimate change remains poorly understood, despite the well-studied link between climate and karst groundwater aquifers. The influence of sea-level rise and coastal geomorphic change on these estuaries obscures climate signals, thus requiring careful development of paleoenvironmental histories to interpret the paleoclimate archives. We used foraminifera assemblages, carbon stable isotope ratios (δ13C) and carbon:nitrogen (C:N) mass ratios of organic matter in sediment cores to infer environmental changes over the past 5300 years in Celestun Lagoon, Yucatan, Mexico. Specimens (> 125 µm) from modern core top sediments revealed three assemblages: (1) a brackish mangrove assemblage of agglutinated Miliammina and Ammotium taxa and hyaline Haynesina (2) an inner-shelf marine assemblage of Bolivina, Hanzawaia, and Rosalina, and (3) a brackish assemblage dominated by Ammonia and Elphidium. Assemblages changed along the lagoon channel in response to changes in salinity and vegetation, i.e. seagrass and mangrove. In addition to these three foraminifera assemblages, lagoon sediments deposited since 5300 cal yr BP are comprised of two more assemblages, defined by Archaias and Laevipeneroplis, which indicate marine Thalassia seagrasses, and Trichohyalus, which indicates restricted inland mangrove ponds. Our data suggest that Celestun Lagoon displayed four phases of development: (1) an inland mangrove pond (5300 BP) (2) a shallow unprotected coastline with marine seagrass and barrier island initiation (4900 BP) (3) a protected brackish lagoon (3000 BP), and (4) a protected lagoon surrounded by mangroves (1700 BP). Stratigraphic (temporal) changes in core assemblages resemble spatial differences in communities across the modern lagoon, from the southern marine sector to the northern brackish region. Similar temporal patterns have been reported from other Yucatan Peninsula lagoons and from cenotes (Nichupte, Aktun Ha), suggesting a regional coastal response to sea level rise and climate change, including geomorphic controls (longshore drift) on lagoon salinity, as observed today. Holocene barrier island development progressively protected the northwest Yucatan Peninsula coastline, reducing mixing between seawater and rain-fed submarine groundwater discharge. Superimposed on this geomorphic signal, assemblage changes that are observed reflect the most severe regional wet and dry climate episodes, which coincide with paleoclimate records from lowland lake archives (Chichancanab, Salpeten). Our results emphasize the need to consider coastal geomorphic evolution when using epikarst estuary and lagoon sediment archives for paleoclimate reconstruction and provide evidence of hydroclimate changes on the Yucatan Peninsula.

modern core top sediments revealed three assemblages: (1) a brackish mangrove assemblage of agglutinated Miliammina and Ammotium taxa and hyaline Haynesina (2) an inner-shelf marine assemblage of Bolivina, Hanzawaia, and Rosalina, and (3) a brackish assemblage dominated by Ammonia and Elphidium. Assemblages changed along the lagoon channel in response to changes in salinity and vegetation, i.e. seagrass and mangrove. In addition to these three foraminifera assemblages, lagoon sediments deposited since 5300 cal yr BP are comprised of two more assemblages, defined by Archaias and Laevipeneroplis, which indicate marine Thalassia seagrasses, and Trichohyalus, which indicates restricted inland mangrove ponds. Our data suggest that Celestun Lagoon displayed four phases of development: (1) an inland mangrove pond (5300 BP) (2) a shallow unprotected coastline with marine seagrass and barrier island initiation (4900 BP) (3) a protected brackish lagoon (3000 BP), and (4) a protected lagoon surrounded by mangroves (1700 BP). Stratigraphic (temporal) changes in core assemblages resemble spatial differences in communities across the modern lagoon, from the southern marine sector to the northern brackish region. Similar temporal patterns have been reported from other Yucatan Peninsula lagoons and from cenotes (Nichupte, Aktun Ha), suggesting a regional coastal response to sea level rise and climate change, including geomorphic controls (longshore drift) on lagoon salinity, as observed today. Holocene barrier

Introduction
Epikarst estuaries provide sediment archives that record sea level rise, storm histories, and coastal environmental change. Like siliciclastic estuary systems, epikarst estuaries back-step landward in response to Holocene sea-level rise and changes in coastal geomorphology such as barrier island progradation. However, epikarst estuaries are generally linked to a groundwater system, and hence modified by catchment hydrochemistry, and thus also respond to ocean-atmosphere forced changes in hydroclimate. Epikarst groundwater systems are important geologic features of the Caribbean region (van Hengstum et al. 2011(van Hengstum et al. , 2020Collins et al. 2015;Tamalavage et al. 2018), where environmental change affects humans and natural ecosystems dependent on karst groundwater resources (Gabriel et al. 2009;Perry et al. 2011;Peros et al. 2015;Jaijel et al. 2018). Coastal epikarst systems have recently been recognized as important paleoclimate archives (Gabriel et al. 2009;van Hengstum et al. 2018), but interpreting climate signals in these complex systems requires a robust understanding of the sedimentary processes that occur in epikarst estuaries and associated lagoons.
Here we present the longest paleoenvironmental record to date from Celestun Lagoon, a large epikarst estuary (sometimes locally called caletas) on the Yucatan Peninsula, Mexico. The microfossil and vegetation records span the last 5300 years and reflect changes in buildup of dune-ridge complexes and lagoon landward back-step progression, development of mangrove forest, storm and current events, and groundwater spring discharge and reveal an estuary's sensitivity to climate and coastal changes in the peninsula lowlands.
The Yucatan Peninsula is a carbonate platform that is dominated by karst geomorphology (Perry et al. 1995;Perry et al. 2003), where rainfall percolates rapidly into groundwater (within months), flows toward the coasts and discharges as submarine springs (Herrera-Silveira 1994;Perry et al. 2003). In the north Yucatan Peninsula, the discharge location of these submarine springs is strongly controlled by the Chicxulub impact crater, which focuses discharge into the Celestun and Dzilam Lagoons (Perry et al. 2002;Perry et al. 2003;Bauer-Gottwein et al. 2011). Climate on the Yucatan Peninsula is seasonal subtropical, with most rainfall occurring from June to October (Herrera-Silveira 1994;Bauer-Gottwein et al. 2011). Seasonality of precipitation is driven primarily by intra-annual migration of the Intertropical Convergence Zone about the equator (Hodell et al. 1991;Hughen et al. 1996;Haug et al. 2001), but inter-annual variability in rainfall is influenced by additional ocean-atmosphere circulation patterns, including the North Atlantic Subtropical High pressure zone in the Atlantic, Walker circulation across the Pacific, and solar insolation across both ocean basins (Bond et al. 2001;Laskar et al. 2004;Sutton and Hodson 2005;Douglas et al. 2016a;Bhattacharya et al. 2017). The peninsula is also characterized by steep gradients in precipitation, with rainfall up to 2000 mm y −1 in the southern highlands to less than 500 mm y −1 in the northern lowlands (Bauer-Gottwein et al. 2011;Metcalfe et al. 2015).
Regional climate variability, together with local hydrologic conditions, is expressed in paleoclimate records across the Yucatan Peninsula from lake and cenote sediment cores (Hodell et al. 1995(Hodell et al. , 2005Leyden et al. 1996;Douglas et al. 2016b), Caribbean coastal cave deposits (Gabriel et al. 2009;van Hengstum et al. 2010;Collins et al. 2015), and speleothems (Webster et al. 2007;Medina-Elizalde et al. 2010;Kennett et al. 2012;Douglas et al. 2016b;Akers et al. 2016). These records differ in the timing and severity of recorded drought events, making it difficult to reconstruct past atmospheric circulation and raising the need for an integrated spatio-temporal rainfall record. New records from coastal lagoons that represent integrated regional climate (via groundwater within a large catchment) could reconcile some of the discrepancies among other paleoclimate archives.
Foraminifera assemblages and bulk sediment organic matter δ 13 C and C:N elemental ratios are useful tools for exploring coastal paleoecology and environmental conditions that respond to climate changes. Modern foraminifera species are distributed along gradients of salinity, wave energy, water stratification, and substrate type, making them sensitive indicators of change in the physical and chemical environment (Debenay and Guillou 2002;Murray 2006;Schönfeld et al. 2012). Coastal paleoenvironmental reconstructions from foraminifera utilize groups or clusters of taxa to represent a particular environment (e.g., brackish estuary, protected lagoon, open coast) and/or individual species to indicate specific habitat (e.g., mangrove swamp, seagrass beds, suboxic sediment) (Sen Gupta 1999). Both methods have been used to characterize the coastal environment in the Yucatan Peninsula and broader Caribbean region. In particular, tropical estuaries are characterized by a landward (estuary's head) finesediment assemblage containing Ammotium, Miliammina, and Trochammina spp., which transitions to brackish mixed assemblages of Ammonia, Elphidium, and Quinqueloculina spp., followed by an ocean assemblage containing Bolivina and Rosalina spp., among other marine taxa toward the estuary opening to the ocean (estuary mouth) (Debenay and Guillou 2002).
In the carbonate-rich mangrove-dominated estuaries of the Caribbean, species richness increases with distance from sheltered coastal lagoons, where a few agglutinated taxa dominate, to tidal flats and a shallow inner shelf, where more diverse hyaline and porcelaneous taxa live. For example, Trochammina and other agglutinated taxa live in fringing mangroves of Puerto Rico whereas Rosalina, Peneroplis, and Archaias live far from protective mangroves in carbonate mud flats and shallow carbonate lagoons (Culver 1990). Changes in the coastal environment affect this distribution, as in Nichupte Lagoon, Quintana Roo (northeast Yucatan Peninsula coast), where barrier islands restricted circulation and changed the foraminifera assemblages from marine genera such as Triloculina and Quinqueloculina to brackish-tolerant Ammonia and Elphidium (Hart and Kaesler 1986). In addition to salinity change caused by reduced connectivity with seawater, restricted circulation results in diurnal or seasonal fluctuations in pH, temperature, and oxygen, with a corresponding decrease in foraminifera species diversity (Hart and Kaesler 1986;Sen Gupta 1999). Circulation restriction has been observed in multiple lagoons along the northeast Yucatan Peninsula coast, with significant changes in sediment deposition and faunal response (Brady 1972). Ecosystem changes based on specific indicator species have also been used around the Caribbean, for example Archaias angulatus and Peneroplis proteus are used to identify coral reef and seagrass habitat because these species require marine salinities and seagrass beds of Thalassia (Sen Gupta 1999; Gischler and Hudson 2004). The agglutinated Miliamina fusca requires hyposaline muddy environments, and assemblages composed entirely of agglutinated specimens often mark inner mangrove swamps (Sen Gupta 1999).
Carbon isotopes (δ 13 C org ) and C:N elemental ratios of sediment organic matter (OM) are robust tracers of the sources of OM to the lacustrine and marine sediments because of the strong difference in δ 13 C org and C:N of terrestrial plants relative to seagrass and marine phytoplankton (Gonneea et al. 2004;Lamb et al. 2006). These variables have been used in circum-Caribbean carbonate environments to infer changes in vegetation and to determine the causes for such changes, including Maya agricultural practices (C3 vs. C4 plants), climate conditions (moist vs. arid), and sea level rise (mangroves) (Gonneea et al. 2004;Polk et al. 2007;van Hengstum et al. 2010). Specifically, δ 13 C, C:N, and microfossil assemblage records were used to elucidate changes in climate variability and its effect on dominant lowlying wetland vegetation and karst sinkhole and cave flooding events during wet climate periods (Polk et al. 2007;van Hengstum et al. 2010van Hengstum et al. , 2020Tamalavage et al. 2018). These studies demonstrated the utility of combining OM variables with other hydroclimate and sea level proxies (e.g., δ 18 O, microfossils, and peat) for documenting ecologic and hydroclimate changes.
In this study we present data on modern benthic foraminifera assemblages and a 5300-year record of foraminifera assemblage changes, δ 13 C, and C:N in sediment from Celestun Lagoon, Yucatan, Mexico. We inferred the paleoenvironmental history of the lagoon with a particular focus on factors that control lagoon salinity. We assessed the causes of the changes in salinity (e.g., sea level rise and changes in coastal geomorphology) and then compared our record to other paleoenvironmental and paleoclimate records from the peninsula.

Site description
Celestun Lagoon is a narrow (0.5-2.4 km) and long (21.2 km) estuarine lagoon on the northwest coast of the Yucatan Peninsula in southern Mexico (Fig. 1a). The lagoon is currently protected from ocean waves by a series of coquina dune ridges about 30 km long (Fig. 1b, Lowery and Rankey 2017). The lagoon interior is characterized by autochthonous carbonate sediments that range in size from fine carbonate muds in the protected northern lagoon to coarse mollusk and gastropod shells and fragments near the southern lagoon, which opens to the ocean. Benthic foraminifera are found throughout the lagoon. The lagoon depth ranges from 0.5 m in the north to 3 m in the south, and bottom sediments in the center of the lagoon remain submerged at low tide (Young et al. 2008;Lowery and Rankey 2017). Red and black mangroves, Rhizophora mangle and Avicennia germinans, respectively, fringe all sides of the lagoon, and abundant, brackish-water-tolerant seagrass beds (Halodule wrightii, Ruppia mexicana) dominate the central lagoon channel (Herrera-Silveira et al. 1998). These polyhaline flora and fauna are supported by submarine groundwater springs that discharge brackish (salinity of 3-9 psu) groundwater predominantly in the northern lagoon, giving the lagoon its estuarine characteristics, with a persistent, but seasonally variable salinity gradient (Herrera-Silveira 1996; Gonneea et al. 2004). Peak groundwater discharge into Celestun Lagoon, which occurs in December, lags peak rainfall by one to two months (Perry et al. 2002), and prolonged wet periods accompany prolonged low-salinity conditions in the northern lagoon (Herrera-Silveira 1994;Herrera-Silveira et al. 1998). The interior of Celestun Lagoon is a low-energy settling basin, and the exterior (ocean side, Fig. 1c) is high-energy, characterized by ocean currents and winter storm waves that cause longshore transport of bioclastic carbonate sediments. Longshore drift transports sediment southward along the outer west bank of the lagoon, resulting in net barrier island  (Hodell et al. 2005), 4 Lake Chichancanab (Hodell et al. 1995 Fig. 2). Core 4A was collected in three drives, and core 4B, 3 m from core 4A, was also collected in three drives (Table 1), but used only for additional radiocarbon dating. Bulk wet density was measured with a gamma densitometer for cores 4A and 4B (Fig. 2) and used in Corelyzer software v.2.0.4 to match radiocarbon ages from core 4B to appropriate depths in core 4A. Bedrock is estimated to be about 3 m deep in the northern lagoon based on clinking sounds heard from coring (Livingstone, 298 cm total recovery, not used in this study) 0.7 km west of 4A and 4B. Auxiliary core UL-5 was collected in 2 drives (200 cm total length, Fig. S1) and was used to corroborate age-depth models and general trends observed in core 4A. Sediment cores 4A and UL-5 were subsampled in 1-cm sections every~5 cm, split into aliquots, and separated for OM δ 13 C and C:N mass ratio analyses and benthic foraminifera assemblages.

Chronology
Basal peat and terrestrial macro-organics (roots, leaves, and seeds) were picked from wet-sieved sediments from cores at sites 4A and 4B and used for radiocarbon dating. Autochthonous organic carbon (e.g., seagrass) was avoided where identifiable because of potential carbon reservoir effects that occur in karst systems (Fletcher et al. 2017). All radiocarbon dates (n=15, Table 2) were calibrated to calendar years using IntCal13 as detailed in Blaauw and Christen (2011). Site 5A was chosen for measuring 210 Pb activity because of its distance from the ocean to the south and from active flamingo flocks to the north, which can enhance bioturbation and sediment mixing. 210 Pb activity (n=7) was determined using 10-mL samples of homogenized bulk sediment sealed in vials for 24 h and measured in a high-purity germanium (HPGe) well detector, using established methods and calibration reference materials (Swarzenski et al. 2006). Sediment accretion rates were calculated from 210 Pb excess . We excluded 210 Pb excess values from the upper mixed 5 cm of the sediment and used values below 6 cm, which show a clear log-linear profile of 210 Pb excess with depth. We set 0 cm depth to −59 years before present, where 1950 is 0 years before present. The age-depth model (Fig. 3) was developed using the R package 'rbacon' v. 2.4.0, which uses Bayesian inference models of the radiocarbon and 210 Pbderived sediment accumulation rate (Blaauw and Christen 2011). Radiocarbon data from sites 4A and 4B and 210 Pb excess values from site 5A were combined to create the most probable age-depth model for core 4A. One radiocarbon age of 1485 at depth of 30.5 cm (Table 2) was excluded from the model because the 210 Pb excess accretion rate above 30 cm and core stratigraphy at 35 cm indicate this radiocarbon sample was not in its original stratigraphic position.

Foraminifera assemblages
Foraminifera were visually identified from~30-mg sediment aliquots, dry sieved to 125 µm (Table 3). Foraminifera between 63 and 125 µm were commonly covered in calcareous precipitates, which could not be dislodged with sieving or sonicating, rendering identification inconsistent or impossible, and hence were excluded. Each species was counted until at least 300 individuals (living or dead) were obtained. Counts were used to calculate relative abundance, and sediment aliquots were weighed to determine absolute abundance (individuals mg −1 of sediment) to address the different volumes of sample used (Schönfeld et al. 2012). This method was applied to coretop samples and core 4A downcore samples. We define "coretop" as the upper 5 cm intervals of sediment, which we assumed were mixed, as inferred from 210 Pb profiles (Fig. 3). In core UL-5, foraminifera counts were performed on the[250-µm fraction. This size fraction was chosen because of the common occurrence of large benthic foraminifera along the Yucatan Peninsula coast (Martinez et al. 2018) and at this site. Data from core UL-5 (Electronic Supplementary Material [ESM]) were not collected at high resolution and were only used to corroborate the trends observed in core 4A. Carbon isotope analysis Organic matter δ 13 C and C:N mass ratios were analyzed simultaneously using a Carlo Erba 1108 elemental analyzer connected to a ThermoFinnigan Delta Plus XP isotope ratio mass spectrometer. A 2-mg sample from each sediment interval was acidified repeatedly using chilled 6% sulphurous acid, following Verardo et al. (1990), with an added sonication between acidifications to promote complete decarbonation. This method minimized effervescence and sample loss in sediments contain-ing~95% carbonate mud (i.e., all samples in this study, excluding peat). Values are reported in reference to VPDB with instrumental reproducibility of δ 13 C=0.2‰ (UC Santa Cruz Stable Isotope Lab). The C:N mass ratio is reported as dry weight % organic C to dry weight % total N, with a reproducibility of 0.1.

Statistical analyses
All statistical analyses were performed in R software (v. 3.5.0) (R Core Team 2020) with packages vegan (v. 2.5.2) (Oksanen et al. 2019) and rioja (v. 0.9.21) (Juggins 2017). Species richness (S), Shannon diversity (H), and evenness (E) were determined from equations presented in Buzas and Hayek (1998). Wisconsin double standardization was applied to raw abundance counts to maximize influence of rare species, and Q-mode and R-mode agglomerative clustering were performed using the Bray-Curtis dissimilarity matrix and Ward's linkage method (Sen Gupta 1999). Foraminifera assemblage names were assigned based on R-mode results, and depth intervals of assemblage occurrence were assigned based on Q-mode results. Cluster cophenetic correlation coefficient, which describes data distortion in cluster dendrograms, was determined using Pearson's correlation between the Bray-Curtis dissimilarity matrix and cluster cophenetic matrix (Davis 2002). Nonmetric multidimensional scaling (NMDS) was performed using the Bray-Curtis dissimilarity matrix of Age-depth model for site 4A cores using radiocarbon dates from 4A (blue) and 4B (light yellow) and 210 Pb accretion-derived age from site 5A core (green, accretion rate from inset). Data is presented in Table 2. Gray dotted lines are the 95% confidence interval for ages. Offset at depth 35 cm represents an inferred erosional event ( Fig. 2) removing 800 years of sediment for site 4A. Inset 210 Pb-derived sediment accumulation rate calculated as -0.0311/-0.21=0.146 cm yr −11 .5 mm yr −1 which is consistent with accretion rates determined for the northern lagoon (Gonneea et al. 2004). More detail on the age-depth models can be found in ESM General salinity tolerance (psu) and environment Reference
Likely similar to T trigonula, 35-81psu Murray (1991Murray ( , 2006 Triloculina linneiana carbonate-rich mud, nearshore marine lagoons and mangroves Javaux and Scott (2003), Murray (2006) J Paleolimnol (2022) 67:131-162  139 the Wisconsin standardized data with k=3 dimensions. Continuous variables were natural logtransformed for plotting as vectors. PERMANOVA and beta-dispersion statistical tests were performed to determine significant differences between cluster centroids and between cluster spreads, respectively (Anderson 2001;Anderson et al. 2006). To check significant difference between centroid positions determined from PERMANOVA -which gives a significant result for centroid position or variance-a randomized resampling was performed to simulate a p value for centroid position alone. Clusters were pooled, randomly sampled without replacement, plotted with NMDS ordination, and distance was measured between NMDS centroids of the randomized clusters. A total of 10,000 iterations were performed and compared to the original centroid distance. Kendall correlation was performed for non-normally distributed continuous variables.

Modern salinity gradient
Lagoon water salinity increases from a minimum of 14 psu in the northern lagoon (with variability near springs) to 33 psu in the southern lagoon, close to the ocean (Table 1). Based on observations in the field, wave energy is lower in the northern interior region (more protected) than in the southern region, close to the coastline.

Sedimentology
Core 4A provided the longest sediment record (262 cm, Fig. 2), hence most analyses were performed on core 4A. The stratigraphy consists of a basal peat layer from 262 to 230 cm, with a transition from 230-228 cm to massive carbonate mud that persist from 228 to 0 cm, with increasing light brown pigmentation upcore (Fig. 2). Fossil bivalves, gastropods, foraminifera, and fragments thereof are rare in peat and common in carbonate mud. Pieces of charcoal and organic matter\1 mm occur occasionally in the carbonate mud. Interval 210-170 cm contains more gastropods than other calcareous microfossils, and interval 170-140 cm contains more carbonate mud with light brown pigmentation and frequent occurrence of gastropods, ostracods, and foraminifera. From 140 to 112 cm, coarse-grained material (1-cm clams) and carbonate mud predominate, and bivalves disappear above 112 cm. Interval 112-0 cm is characterized by carbonate mud, frequent gastropods, and abundant foraminifera. From 36 to 34 cm a coarse (0.5-1.0 mm) layer of foraminifera and gastropods occurs, with a notable color transition from light brown to reddish brown (Fig. 2b). Core 4B exhibits nearly identical stratigraphy (Fig. 2a) and consists of a basal peat layer overlain by massive carbonate mud and fossil bivalves, gastropods, and foraminifera from 230 to 0 cm, with a coarse layer of gastropods and foraminifera from 36 to 34 cm (Fig. 2b). Wet bulk density values are similar between the two cores, with 1.3 g/mL in the peat layer, a maximum of 1.9 g/mL just above the peat-mud transition (229 cm), and a trend to 1.4 g/mL at the surface (Fig. 2a). Based on stratigraphy, wet sediment density, and core General salinity tolerance (psu) and environment Reference
The Bayesian age-depth model for core 4A (Fig. 3) implies a long-term mean linear sediment accumulation rate between 0.4 and 0.6 mm y −1 . Assuming the lagoon-wide coarse sediment layer at~35 cm depth ( Fig. 2b) is an erosional event, the age model for core 4A indicates a missing time span of 800 years, from 1300 to 500 calendar years Before Present (hereafter BP). A comparison of age models with and without the erosional event and 210 Pb excess data can be found in ESM. The model with an erosional event and with 210 Pb excess accretion rates, though complex, better reconciles radioisotope and sedimentological observations and more accurately reflects age model uncertainty for this dynamic coastal site than does a simpler radiocarbon model.

Benthic foraminifera identification and distribution
More than 27,000 foraminifera were counted in 90 samples, and 26 species were identified (Table 3). Foraminifera appeared pristine, with no pitting or dissolution noted in coretop or downcore specimens.
Diversity, clustering, and statistical relationships R-mode cluster analysis of fossil and modern foraminifera indicates five clusters of commonly cooccurring taxa with a cophenetic correlation coefficient of 0.72 (Fig. 5a). We refer to these clusters by simplified assemblage names, which are Bolivina-Hanzawaia-Rosalina, Trichohyalus, Archaias-Laevipeneroplis, Ammonia-Elphidium, and Agglutinated-Haynesina. All five assemblages were observed in core 4A, whereas assemblages Archaias-Laevipeneroplis and Trichohyalus were absent from modern coretop sediments.
Diversity indices S, H, and E change with the faunal succession upcore in 4A (Fig. 4a). From Phases 1 to 4 (dashed lines in Fig. 4a), respective mean S changed from 8.8-11.8 to 8.9-7.6, mean H from 1.05-1.31 to 0.44-0.52, and mean E from 0.36-0.34 to 0.44-0.52. A maximum S of 18 occurs at 229 cm (4900 BP) and a minimum of 6 occurs at 170 cm (3400 BP). Diversity H follows a similar trend, but remains high for the last 2900 BP, due to an increase in absolute abundance, and E increases over time. In modern samples, the northern lagoon has values of S=9, H=1.64, and E=0.57, and the southern lagoon has values of S=11, H=1.71, and E=0.53 (Fig. 4b).
NMDS of core 4A clusters and environmental vectors are presented in Fig. 6 and represent statistical relationships given in Table 4. All four clusters plot without overlap in ordination space (Fig. 6A), and foraminifera composition of each cluster is distinct and significantly different ( Fig. 6C and Table 4, PERMANOVA Source: Cluster, p=0.003; Bootstrapping Source: Cluster, p\ \0.01) with no significant difference in cluster dispersion (Beta dispersion Source: Cluster, p=0.111). Bulk sediment δ 13 C org and C:N ratio significantly correlate with cluster composition (p=0.001 for both). Significant correlation between δ 13 C org and E (Kendall τ= −0.291, p=0.014) and S (Kendall τ=0.370, p= 0.003) also show a significant, but weak relationship. S and H showed no significant correlations, but plot as long vectors in the general direction of δ 13 C org (Fig. 6a).
NMDS of modern coretop Clusters 1 and 2 show no overlap (Fig. 6b) and are significantly different in assemblage composition ( Fig. 6d and Table 4, PERMANOVA Source: Cluster, p=0.001; Bootstrapping Source: Cluster, p\ \0.01), but are not significantly different in the variance of species abundance (Beta dispersion Source: Cluster, p= 0.128). The largest gradients (longest vectors) along the lagoon transect are δ 13 C org , C:N ratio, and salinity (Fig. 6b). PERMANOVA results are significant for salinity, δ 13 C org , and cluster composition (p=0.001, Table 4), and temperature is significant (p=0.014). Species richness (S), δ 13 C org , and C:N ratios are significantly correlated (Kendall correlation, τ=0.− 460 and 0.46, respectively, p=0.002 for both), but no significant correlation is observed between the δ 13 C org and assemblage diversity (H) or evenness (E). S is significantly correlated with distance along the transect from the northern to southern lagoon (Table 4), coincident with the observed decrease in lagoon wave energy from north to south. Results and figures related to auxiliary core UL-5 are described in the ESM.

Distributions of flora and fauna
Foraminiferal and geochemical analyses in surface sediments and water along Celestun Lagoon reveal differences between the northern and southern lagoon regions. The Agglutinated-Haynesina assemblage occurs only in the northern lagoon, whereas the Bolivina-Hanzawaia-Rosalina assemblage occurs only in the southern lagoon (Fig. 4b). The Ammonia-Elphidium assemblage is common in both regions. The Trichohyalus and Archaias-Laevipeneroplis assemblages were not observed in the modern lagoon. Maximum species richness (15) and diversity (1.96) both occur in the middle of the lagoon, about 13 km from the southern opening, where assemblages overlap (Fig. 4b) and the habitat transitions from a protected low-salinity lagoon to a higher-energy marine-influenced lagoon. The assemblage distribution aligns with gradients in salinity, temperature, δ 13 C, and C:N of OM (Fig. 6b). Strong association (clustering) of foraminifera assemblages with specific salinity and carbon isotope ranges suggests that the distribution of benthic foraminifera in Celestun Lagoon is controlled in part by salinity, vegetation, and, to a lesser extent, temperature. PERMANOVA A Fig. 5 Stratigraphic diagrams of benthic foraminifera colored by cluster results (from Fig. 5) and paired with diversity indices and bulk sediment organic carbon isotope data for A site 4A cores, with an age scale obtained from Fig. 3, and B coretop samples from north to south in the lagoon (Fig. 1). Dashed lines represent stratigraphic transitions between assemblages. Note change in scale on the bottom axis of each panel and the age jump at 35 cm reflecting a missing time span from 1300 to 500 cal yr BP results support significant relationships among these variables (Table 4). Clusters dominated by the Agglutinated-Haynesina assemblage represent brackish salinity and shallow lagoon mangroves, and clusters dominated by the Bolivina-Hanzawaia-Rosalina assemblage represent inner marine shelf and seagrass habitat, which are the environments in which species of these assemblages have been found around the Caribbean (Table 3). Low wave energy and muddy substrate in the enclosed northern lagoon also promote the abundance of agglutinated foraminifera like Miliammina fusca (Debenay and Guillou 2002). Lagoon water temperature is weakly correlated with cluster composition (Table 4, Fig. 6b), but does not exhibit large changes in mean temperature among sites sampled (29.2± 2.6°C, Table 1). Salinity exhibits the greatest variability in the lagoon. Although Trichohyalus and Archaias-Laevipeneroplis were not observed in modern Celestun, we infer their habitats based on modern occurrence at other sites in the Gulf of Mexico and Caribbean. Trichohyalus aguayoi occurs in inland mangrove ponds in Bermuda (Javaux and Scott 2003) and groundwater-flooded blue holes in the Bahamas (van Hengstum et al. 2020), and thus represents lowsalinity, isolated mangrove ponds with minimal ocean influence (Table 3). Archaias angulatus and Laevipeneroplis proteus (also known as Peneroplis proteus) inhabit marine seagrass of the genus Thalassia off the coasts of Belize and Florida (Hallock and Peebles 1993;Sen Gupta 1999). Thus, the Archaias-Laevipeneroplis assemblage represents marine-salinity, shallow-shelf environments with Thalassia seagrass vegetation (Table 3).
In the marine-brackish conditions of southern Celestun Lagoon, euryhaline Ammonia and Elphidium were found alongside epiphytic Quinqueloculina and Rosalina taxa, suggesting seagrass vegetation may contribute to assemblage diversity, but alone is insufficient in determining modern assemblage distribution. The importance of salinity is further demonstrated by benthic foraminifera observed 10 km south of Celestun Lagoon, in Isla Arena, which has similar sedimentology, geometry, and vegetation as Celestun, but higher salinity because of greater seawater input (Lowery and Rankey 2017).
As a result of this higher salinity, the Isla Arena assemblages are dominated by miliolid and peneroplid foraminifera that prefer near-marine salinity.   Variable salinity in coastal environments restricts foraminifera occurrence to only those that can withstand a large and fluctuating salinity range (Sen Gupta 1999). Such conditions occur in the northern lagoon and result in low species richness (low S values,Figs. 4b,6b). The pattern of assemblage change along a salinity gradient is found in other Caribbean and Atlantic coastal systems with large salinity ranges. On the mangrove forest coast of French Guiana,  noted that seasonal salinity oscillations from 20 to 35 psu are accompanied by assemblage compositional changes from Miliammina-dominated (agglutinated) to Ammonia-and Elphidium-dominated. A narrow estuary of Brazil exhibited shallow protected brackish-water assemblages of agglutinated Miliammina and Ammotium taxa and seaward marine-influenced assemblages dominated by A. tepida, E. gunteri, and Bolivina taxa (Debenay et al. 1998), similar to the southern opening of Celestun Lagoon. Collectively, these studies point to salinity as being the most important environmental variable controlling foraminifera distribution in Celestun Lagoon. Based on these associations, we interpret changes in downcore assemblages as representing changes in salinity in the northern lagoon over time.
We acknowledge that analyzing the[125-µm size fraction may exclude some foraminifera species or smaller adult specimens of the species we identified (Table 3, Murray 2006). However, benthic foraminifera at our site and other coastal Yucatan Peninsula sites are characterized by abundant large ([250 µm) specimens (Martinez et al. 2018), suggesting we observed most taxa present. Based on the number of taxa identified (n=26 , Table 3), the number of specimens counted (n=27,000), and the striking differences in assemblage composition observed in the[125-µm fraction (Fig. 4), we consider the statistical relationships noted between environmental parameters and assemblage composition (Table 4) to be real patterns for the modern system, despite exclusion of small foraminifera. Furthermore, δ 13 C and C:N of bulk sediment OM are independent of foraminifera size fractions and are proxies for mangrove and seagrass vegetation that are also distributed along the lagoon salinity gradient (Gonneea et al. 2004;Khan et al. 2017;Carnero-Bravo et al. 2018), supporting interpretation of environments inferred from assemblages. It is possible that drying the samples prior to counting may affect preservation of agglutinated taxa (Murray and Alve 1999) and bias assemblage composition. However, the hyaline foraminifera Haynesina germanica observed in the low-energy muddy substrate preferred by agglutinated taxa provides an additional A B C D Fig. 7 Conceptual diagram of coastal development around Celestun Lagoon. The star represents site 4A in the present day northern lagoon. A Phase 1, during low sea level the core site is located within a inland low-salinity mangrove pond with Trichohyalus assemblage. B Phase 2, sea level rises beginning 5100 BP, the site is inundated, and coastal marine Thalassia seagrass and Archaias-Laevipeneroplis assemblages colonize by 4900 BP. Long-shore drift continues, initiating spits and barrier islands which begin a prolonged period of assemblage change. C Phase 3, decelerating sea-level rise allows for spits and barrier islands to accumulate stably. Bolivina-Hanzawaia-Rosalina dominates the site but is replaced by Ammonia-Elphidium as groundwater lowers salinity the nascent lagoon. D Phase 4, sea level is stable, barrier island coastline expands, and mangroves colonize new shoreline in the lagoon which increases in length. Agglutinated-Haynesina assemblage develops alongside Ammonia-Elphidium. Barrier islands continue to accumulate to present. Precipitation (climate) recharges groundwater throughout all phases, and final assemblage reflects reduction in mixing between the groundwater and seawater as well as absolute decrease in groundwater (and increase in sweater) during droughts constraint on the distribution of the Agglutinated-Haynesina assemblage in space and time. We thus are confident in using assemblages defined in the modern environment to infer downcore, paleoenvironmental conditions, specifically changes in salinity.

Lagoon salinity changes-geomorphic processes and climate changes
The salinity in a restricted lagoon such as Celestun Lagoon can vary as a consequence of two distinct processes that affect the relative contribution of seawater and fresh/brackish groundwater inputs to the lagoon. The first is physical (geomorphic) changes to the coastline (i.e., lagoon restriction or sea level changes) that determine the amount of seawater input to the lagoon, and the second is climate changes that impact regional and local precipitation, which control the amount of groundwater input.
Coastal lagoon morphology in the Yucatan Peninsula was studied by Brady (1972), which remains the most thorough compilation of epikarst sedimentological processes for the peninsula. Brady indicates that modern barrier islands evolved from Pleistocene aeolian dunes on the continental shelf. Once formed, dunes were submerged by rising seas and provide shallow sites to initiate submarine barrier island formation and subaerial aeolian accumulation. Island growth may be episodic, particularly when hurricane overwash adds or removes sediments, but mangrove encroachment rapidly stabilizes island shorelines and facilitates additional shoreline sediment accumulation. Importantly, sediment cores of Nichupte and Holbox Lagoons (Fig. 1c) reflect sharp contacts between basal mangrove peat and overlying carbonate mud, indicating sea level rise and lagoon formation, but are subsequently characterized by faunal shifts from inner-shelf miliolids and clams to restricted lagoon Elphidium and encrusting algae, along with a fining grain size (Brady 1972). Faunal change was attributed to barrier island growth, but constraints on accumulation rate were poor.
More recent work indicates that longshore sediment transport along the Yucatan Peninsula coastline can build barrier island complexes within decades. Based on the observed net 200 m per 33 years of barrier island progradation offshore near Celestun Lagoon (Fig. 1b, Lowery and Rankey 2017), a 600 m barrier can develop within a century, and around 30 km can accumulate in 5000 years. Volumetric rates of sediment transport (35,000 m 3 yr −1 ) along the northwest Yucatan Peninsula coast are sufficient to create the estimated volume of sediment of the Celestun west bank (5.3910 7 m 2 93 m elevation) on this timescale (Appendini et al. 2012). These rates of barrier island development are consistent with gradual accretion of the west Celestun Lagoon bank, resulting in an increase in lagoon length over time and a parallel gradual decrease in the mixing and input of seawater into the northern lagoon. Specifically, sediment deposition at the Celestun coast has resulted from both a Pleistocene dune high point, as indicated by Brady (1972), and a change of shoreline orientation relative to the west-flowing Yucatan Current. Cuevas et al. (2013) noted that offshore dunes tend to accumulate at abrupt changes in shoreline orientation on the north coast of the Yucatan Peninsula, and both the northwest (Celestun) and northeast (Holbox) corners of the Peninsula exhibit extensive barrier-dune complexes.
We note similar barrier-dune complexes from Google Earth™ near lagoons of Holbox, Ria Lagartos, Dzilam de Bravo, and Progreso, on the north coast of the Yucatan Peninsula, and near Isla Arena and Terminos on the west coast (Fig. 1c). In Terminos Lagoon, the dune-barrier complex called Isla de Carmen consists of 38 km of shelly coquina (Phleger and Ayala-Castanares 1971). Based on the progradation rate determined by Lowery and Rankey (2017), Isla de Carmen took 6300 years to accumulate, consistent with the island's estimated age of 7000 years based on Holocene sea level terraces (Phleger and Ayala-Castanares 1971). These studies all indicate that the rate and scale of longshore transport and barrier island formation observed today were similar over the Late Holocene.
Although the gradual restriction of seawater into Celestun Lagoon will result in a gradual decrease in salinity, variability superimposed on this trend may be linked to climate-induced changes in groundwater input. Specifically, an increase in the salinity would suggest lower precipitation and possible drought conditions. J Paleolimnol (2022) Fig. 1 for site locations. A Relative sea-level and 95% confidence shading (Khan et al. 2017). B Summer (August) insolation at 20°N (Laskar et al., 2004). C Hydrogen isotopes of plant wax from Lake Salpeten (Douglas et al. 2015). D NMDS of microfossils and Cenote Aktun Ha (van Hengstrum et al. 2010). E δ 18 O of Pyrgophorus sp. in Lake Punta Laguna (Curtis et al. 1998). F δ 18 O of Pyrgophorus coronatus in Lake Chichancanab (Hodell et al. 1995). L NMDS scores of foraminiferal assemblages (see Fig. 6a). M H diversity (reflecting total species and absolute abundance) of assemblages. N Phase of lagoon development for Celestun Lagoon (see conceptual diagram, Fig. 7). 14 C data points are for this study (Table 2). Shaded columns represent documented dry periods for the north-central Yucatan Peninsula. Dotted box represents possible storm or tsunami deposit in coastal records from the west, south, and east of Celestun Lagoon (Fig. 1). Note that NMDS scores in D and l cannot be directly compared and are presented to show relative trends

Reconstruction of the northern lagoon environmental setting
Based on comparison between observations from core-top sediments along the present-day lagoon transect and downcore data from core 4A in the northern lagoon, we suggest that the northern lagoon environment has experienced four phases of development, characterized by distinct environmental settings controlled by climate, sea-level rise, and coastal geomorphic change (Fig. 7). Based on general agreement between cores 4A and UL-5 (ESM), these phases are representative of the entire northern lagoon.
Phase 1 (5300-5100 BP) Between 5300 and 5100 BP, the northern area of Celestun Lagoon was within a mangrove forest that was isolated from ocean influence. The peat δ 13 C org value of − 28.26‰ indicates mangroves and other terrestrial vegetation (Gonneea et al. 2004;Young et al. 2005;Lamb et al. 2006;Tamalavage et al. 2018) (Fig. 2a), and peat accumulation suggests the absence of oxidants like dissolved oxygen and marine sulfate, consistent with a water-logged freshwater swamp setting. The oligohaline assemblage Trichohyalus observed during this interval is found today in isolated, brackish, groundwater-fed pond environments in Bermuda, supporting interpretation of minimal marine influence (Javaux and Scott 2003).
Regional sea level at that time is estimated to have been 2-3 m lower than today (Khan et al. 2017), and a bathymetric slope of − 0.4 m km −1 off the northwest Campeche Bank near Celestun (Appendini et al. 2012) would place the coastline 5 km to the northwest of the current location of Celestun Lagoon, consistent with reduced marine influence. Groundwater was likely provided via springs in the permeable Ring of Cenotes because this Cretaceous-age feature existed prior to Holocene coastal development of the Yucatan Peninsula (Perry et al. 2002). The environment may have resembled the "Petenes" that are present south of Celestun Lagoon today, a periodically flooded black mangrove forest characterized by other marsh and subtropical vegetation, diffuse groundwater inputs, organic-rich sediments, and gastropods (Lowery and Rankey 2017). Based on the low-salinity habitat of Trichohyalus in Bermuda (Javaux and Scott 2003) and modern groundwater discharge salinity values (Herrera-Silveira 1996), salinity was likely less than 10 psu in these mangrove ponds and probably closer to 3 psu, the value of modern groundwater at this site (Young et al. 2008). We know of no published climate records that span the period between 5300-5100 BP in the northwest Yucatan Peninsula. However, pollen records suggest mesic (semi-moist) tropical forests in the southwest and northeast areas of the Yucatan Peninsula and in Haiti ca. 5200-4800 BP (Hodell et al. 1991;Higuera-Gundy et al. 1999;Aragón-Moreno et al. 2012;Vela-Pelaez et al. 2018). Regionally wet conditions would provide groundwater to inland mangrove ponds to support Trichohyalus assemblages whose presence in this environment distinguishes this phase from subsequent, more saline conditions in the lagoon.

Phase 2 (5100-3000 BP)
This phase of over 2000 years is characterized by both changes in foraminifera assemblages and sources of organic matter to the lagoon. Between 5100 and 4900 BP, the peat in the core transitions to carbonate mud (Fig. 2a) as δ 13 C org increases from − 28.26 to − 20.21‰ and C:N decreases from 25 to 11 (Fig. 4a). The Trichohyalus assemblage disappears from Celestun Lagoon, and at 4900 BP the Archaias-Laevipeneroplis assemblage appears alongside the Ammonia-Elphidium assemblage (Fig. 4a). The δ 13 C org shift to less negative values indicates an influx of marine particulate matter and seagrass vegetation (Gonneea et al. 2004), corroborated by the Archaias-Laevipeneroplis assemblage which inhabits marine Thalassia seagrass beds (Table 3, Sen Gupta 1999). These data suggest that coastal marine conditions at Celestun Lagoon commenced at this time, probably caused by sea-level rise ca. 5100 BP (Fig. 8). On the Caribbean side of the Yucatan Peninsula, in Puerto Morelos, mangrove peat layers also date to 5000 years ago (Islebe and Sánchez 2002), and coastal caves in Belize and Quintana Roo (Yucatan Peninsula east coast) exhibit basal peats dated to 6000-4000 years ago, with mixed pollen and marine microfossil assemblages suggesting a transition from terrestrial to marine-influenced coastal zones (Brady 1972;Polk et al. 2007;Gabriel et al. 2009;van Hengstum et al. 2010;Collins et al. 2015). A composite sea-level curve indicates that at 5000 BP, sea level was 2.5 m lower than today along the north coast of the Yucatan Peninsula (Khan et al. 2017), consistent with the 2.6 m depth at which peat is found at Celestun (Fig. 2a).
Rising sea level may have been coincident with drought in the northwest Yucatan Peninsula, as seen from records in Lakes Chichancanab and Punta Laguna ( Fig. 8; Hodell et al. 1995;Curtis et al. 1998). The Ammonia-Elphidium assemblage at Celestun suggests brackish groundwater discharge was still present, but limited, otherwise the brackish-intolerant Archaias angulatus would not be observed (Hart and Kaesler 1986). This is also consistent with palynological and sedimentary changes of riverine soils in the south Yucatan Peninsula that suggest the onset of droughts ca. 4800 BP, which would have impacted the catchment area draining to Celestun Lagoon (Bauer-Gottwein et al. 2011;Aragón-Moreno et al. 2018;Vela-Pelaez et al. 2018). These events occur within chronological uncertainty of the appearance of Archaias-Laevipeneroplis (4900 BP) and may also explain the short duration (century-scale drought event) of this assemblage in the lagoon, which would have disappeared once average precipitation had resumed.
By 4500 BP, the Archaias-Laevipeneroplis assemblage was replaced by the Bolivina-Hanzawaia-Rosalina assemblage with increased abundance of the Ammonia-Elphidium assemblage, while δ 13 C remained between − 19 and − 21‰ (Fig. 4a), indicating persistent seagrass and marine vegetation from 4500 to 3500 BP. The foraminifera community of that time interval greatly resembles the modern assemblage in the southern part of the lagoon (Fig. 4b), suggesting that salinity was likely\30 psu, which is inhospitable to Archaias (Hallock and Peebles 1993), but ideal for the other taxa (Table 3). This represents a decrease from a marine salinity supporting Thalassia seagrass and Archaias-Laevipeneroplis at 4800 BP, to salinities of\30 psu, likely caused by increased groundwater discharge or reduced contribution of seawater. Dry conditions persisted from 4600 to 4400 BP in Belize (Akers et al. 2016), at 4500 BP in the south-central Yucatan Peninsula (Curtis et al. 1998;Wahl et al. 2014), from 4700 to 3600 BP in the southwest Yucatan Peninsula (Islebe et al. 2019), though conditions in the north-central Yucatan Peninsula at that time are unclear. Hence, precipitation-induced increase in groundwater discharge in Celestun Lagoon at 4500 BP is at odds with available climate data (Fig. 8).
It is likely that formation of barrier islands and spits around the northern lagoon, where springs are located (Fig. 1a, b), reduced mixing between groundwater and seawater, resulting in the observed decrease in salinity (Fig. 7). Highest species richness and diversity occurred between 4500 and 3800 BP (S =17, H=1.97, Fig. 4a), similar in magnitude to, and comprised of the same taxa, as the modern central lagoon assemblages (Fig. 4b). Indeed, high diversity is expected of foraminifera that colonize the diverse habitats typical of such an environment (fringing mangrove, brackish-tolerant seagrass, semi-protected lagoon, and unrestricted seawater circulation). West of Celestun Lagoon, across the Gulf of Mexico (Fig. 1c), La Mancha Lagoon in Veracruz, Mexico, also exhibits a rapid change in foraminifera composition between 6500 and 4500 BP, attributed to longshore transport and barrier island formation (Arellano-Torres et al. 2019). Initiation of coastal dunes and barrier island formation across the Gulf of Mexico, as sea-level rose, is consistent with apparent freshening and maximum diversity in Celestun Lagoon, despite regionally dry conditions (i.e., the restriction of seawater mixing is more effective in controlling the salinity than the reduction in groundwater input caused by drought). A more recent example (1000 BP) of this process was observed in sediment cores from Nichupte Lagoon on the eastern Yucatan Peninsula (Fig. 1c), where assemblages shifted from an Archaias assemblage to an Ammonia-Elphidium assemblage as barrier islands isolated the lagoon and reduced mixing with seawater (Hart and Kaesler 1986).
From 4500 to 3600 BP, diversity decreased as marine taxa progressively disappeared, while mean δ 13 C (− 19.5‰) and C:N (8.2) changed little (Fig. 4a). the Ammonia-Elphidium assemblage became dominant by 3400 BP, while the Bolivina-Hanzawaia-Rosalina assemblage decreased in abundance, indicating salinity ranged from 20-25 psu based on modern lagoon distributions (Table 3, Fig. 4b). We note, however, that sea level rise during that time (Fig. 8) should have resulted in greater marine influence, not less. Hence, we interpret the shift towards brackish taxa as indicative of greater influence of groundwater discharge, which may have been caused by increasingly moist climate between 4500 and 3400 BP.
The lower salinity could also be attributed to reduced contribution of seawater related to sand-bar buildup, however rapid fluctuation between dominant assemblages suggests an unstable environment, more likely reflecting climate change rather than more gradual and unidirectional barrier island accumulation. Between 3600 and 3400 BP (Fig. 4a), assemblages oscillated between Bolivina-Hanzawaia-Rosalina and Ammonia-Elphidium, accompanied by large changes in S and H. Specifically, because salinity correlates best with assemblage distribution (Table 4), shifting relative abundances between these two assemblages suggest fluctuating lagoon salinity, and is not consistent with the generally unidirectional formation of barrier islands. Instead, the lagoon geometry at that time was such that relatively small changes in absolute groundwater discharge greatly affected the salinity in the lagoon, thereby impacting foraminifera assemblages (Fig. 7). These relatively short-term salinity fluctuations were likely related to climate changes, such as variability in the ITCZ position and regional precipitation, as manifested by high variability in other climate records during this period ( Fig. 8; Hughen et al. 1996;Douglas et al. 2016b;Aragón-Moreno et al. 2018). A wet period ca. 3500 BP is supported by negative oxygen isotope excursions in lake records from Lake Salpeten and mesic forest pollen from Lake Silvituc in the southwest peninsula (Wahl et al. 2006;Douglas et al. 2015;Vela-Pelaez et al. 2018). However, dry conditions are indicated at the same time by the oxygen records from Lakes Chichancanab and Puerto Arturo in the central peninsula (Hodell et al. 1991;Wahl et al. 2014). Indeed, at 3500 BP our data best match the Fe/Ca data (a runoff proxy) of La Mancha Lagoon, west of Celestun Lagoon (Figs. 1c, 8), across the Gulf of Mexico, demonstrating the regional link between climate and assemblage variability in the Gulf of Mexico southern coasts.
In Celestun Lagoon, from 3400 to 3000 BP, the Ammonia-Elphidium assemblage became dominant and Bolivina-Hanzawaia-Rosalina decreased in prevalence, indicating a reduction in salinity, while a decrease in δ 13 C suggests more mangrove-derived carbon in sediments (Fig. 4a). It is difficult to say whether these salinity and carbon trends are the result of barrier island growth or climate. Mangrove forest cover in particular may respond to both wetter climate and/or to additional shoreline to colonize. K/Ca ratios in nearby Los Petenes and Fe/Ca ratios in La Mancha Lagoon suggest reduced erosion and hence dry conditions across the southern Gulf of Mexico (Roy et al. 2017;Arellano-Torres et al. 2019), whereas δ 18 O records from Punta Laguna and δD records from Lake Salpeten imply wet conditions over the east and south Yucatan Peninsula around 3000 BP ( Fig. 8; Curtis et al. 1998;Douglas et al. 2015). The general divide between dry northwest and wet southeast among records in the region may suggest that a specific atmospheric circulation pattern drove precipitation at that time. Indeed, previous studies indicate changes in the positions of the ITCZ, North Atlantic Subtropical High, and descending limb of Walker circulation (Metcalfe et al. 2015;Bhattacharya et al. 2017).
Overall, Phase 2 is distinguished by the apparent faunal succession from coastal marine taxa to brackish lagoon taxa (Fig. 4a), despite a sea-level rise of 1 m (Khan et al. 2017), with occasional oscillations between dominant assemblages during periods of climate variability (Fig. 8). Close agreement among observed sediment transport rates (Appendini et al. 2012), calculated barrier island aggradation time (Lowery and Rankey 2017), and transition between dominant foraminifera assemblages (Fig. 4a), suggests coastal geomorphology continued to influence the hydrogeology and salinity of Celestun Lagoon during Phase 2, with assemblages exhibiting variable sensitivity to changes in groundwater discharge induced by climate at that time (Fig. 7).

Phase 3 (3000-1700 BP)
Between 3000 and 1700 BP, species richness decreased and assemblages were dominated by the near-exclusive presence of Ammonia-Elphidium (Fig. 4a). Bulk sediment δ 13 C org decreased from − 19.51 to − 23.27‰ by 1700 BP, indicating vegetation shifted from seagrass/phytoplankton to a greater proportion of mangroves. The decrease in species richness indicates a more homogenous setting (fewer habitats) compared to Phase 2 and was probably similar to the modern northern lagoon (Fig. 4b).
Celestun Lagoon salinity during this phase was likely15-25 psu, indicated by the Ammonia-Elphidium assemblage and decrease in species richness (Fig. 4b). This trend is consistent with the ongoing barrier island buildup that reduced seawater input. However, at 2500 BP (105-106 cm), Quinqueloculina bicarinata and Hanzawaia concentrica appear briefly and Ammonia-Elphidium abundance decrease, suggesting a reduction in groundwater and more marine influence in the lagoon (Fig. 4a). This change corresponds to drought conditions noted in Lake Chichancanab (Central Yucatan Peninsula), and Macal Chasm (Belize) speleothem records ca. 2600 BP (Hodell et al. 1991;Akers et al. 2016). The generally drier conditions on the Yucatan Peninsula from 3000 to 2000 BP (Douglas et al. 2015) also match the timing of increased marine influence in Celestun Lagoon, as suggested by slight carbon isotope changes from -21.7 to -21.5‰ (indicating greater input of marine particulate matter). After 2000 BP, Ammonia-Elphidium returned, and the conclusion of this phase is the onset of the brackish, semi-protected lagoon environmental conditions that characterize the present-day northern lagoon.

Phase 4 (1700 BP-present)
The last 1700 years of northern Celestun Lagoon history are characterized by the appearance of the Agglutinated-Haynesina assemblage, low species richness (S=7), δ 13 C org decrease from − 22.97 to − 24.22‰, and variable C:N (mean=7.5). Elphidium galvestonense disappeared, but other taxa of the Ammonia-Elphidium assemblage (Fig. 5c) remained, alongside increased abundance of Q. laevigata. These conditions persist to the present, thus representing the onset of modern conditions in the northern lagoon beginning at 1700 BP. The Agglutinated-Haynesina assemblage indicates increased muddy substrate, which is observed in the protected mangrove lagoon today and is consistent with mangrove forest expansion, indicated by decreasing carbon isotope values (Fig. 4a). Mangrove expansion around Celestun Lagoon is consistent with expansion of mangrove forest documented elsewhere in the Caribbean (Peros et al. 2007). Though barrier islands provide more land for mangrove colonization, the effect is local, so a regional expansion of mangroves suggests regional climate influence. Onset of moister climate ca. 1700 BP is consistent with observation in the Chichancanab and Salpeten lake records ( Fig. 8; Hodell et al. 1995;Douglas et al. 2015).
The age-depth model indicates that the coarse bioclastic layer at 35 cm depth (Fig. 2b) begins at 1300 BP, or about 700 Common Era (CE). In core 4A this layer consists of an increase in absolute abundance of foraminifera in the Ammonia-Elphidium assemblage and a pronounced decrease in agglutinated taxa and H and E indices (Fig. 4b). In the southern lagoon, abundant bivalves and fragments are noted. This interval coincides with the well-documented Maya droughts (Douglas et al. 2016b), which implies decreased groundwater discharge in Celestun Lagoon. Although an increase in the polyhaline assemblage Ammonia-Elphidium could be a response to increased salinity, the lack of fine mud across the cores would imply that authigenic carbonate precipitation was reduced in the lagoon at the same time. If included in the age-depth model, radiocarbon sampled at 31 cm (Table 2, 1485 14 C yr) causes a rapid change in slope in the age-depth regression line that results in abrupt decrease in sedimentation to 0.2 mm y −1 . This slowed sedimentation rate is not consistent with the 1.5 mm y −1 rate determined by 210 Pb excess (Fig. 3) or rising sea level at this coastal setting (Fig. 8, Khan et al. 2017). We thus interpret this lagoon-wide layer to be caused by storm erosion or a tsunami. Hurricanes are known to sort and transport shelf foraminifera into coastal systems (Murray 2006), and identified features of this coarse layer, including the lack of fine sediments or agglutinated taxa, are consistent with such an event. The rapid increase in absolute abundance of Ammonia and Elphidium, which occurs in hurricane overwash deposits on the Atlantic coast of the United States (Collins et al. 1999), also supports this interpretation. The presence of an erosive shell layer representing the missing interval, 1300 to 500 BP, provides some paleoclimate information, as other studies suggest increased hurricane activity from 800 to 1350 AD (1150-600 BP) in the Caribbean (Peros 2015), although our age model uncertainty is high near the shell layer at 35 cm (Fig. 3).
Alternatively, the bioclastic layer could represent a tsunami deposit, likely from submarine continental slope failure (Horrillo et al. 2013). Holocene sealevel rise has been suggested as a cause of slope failure from water overburden pressure (Smith et al. 2013), and in the Gulf of Mexico, sea level has risen 9 m during the last 7.8 ka (Blum et al. 2002). Continental slopes may also fail because of additional sediment load, which would be expected from increased continental weathering and riverine transport to coasts. Fe/Ca and K/Rb ratios in sediments of La Mancha Lagoon (Fig. 1c), indicate increased detrital input and coastal weathering at 1000 BP (Arellano-Torres et al. 2019). The timing is coincident with the high-energy shell deposit in Celestun Lagoon (Figs. 2, 3), and evidence of a tsunami deposit might be recognized in sediment cores from other lagoons around the southern Gulf of Mexico. Brady (1972) noted an erosional contact with increased grain size at 20 cm depth in Holbox Lagoon (Fig. 1c). South of Celestun in Los Petenes (Fig. 1c), sediments transition from peat at 38 cm to silt between 35 and 15 cm and are dated around 700 BP (1300 CE) with a spike in K/Ca that indicates rapid erosion inland (Roy et al. 2017;Fig. 8). Rapid erosion from rainfall is unlikely during the Terminal Classic Maya mega-droughts of this time (Hodell et al. 1995;Douglas et al. 2016b), further implying a catastrophic depositional event such as a tsunami inundation and subsequent outwash. Lastly, an anomalous berm of unsorted boulders and sand on the Caribbean coast of the Yucatan Peninsula has been interpreted as a potential tsunami deposit with a radiocarbon age of 1500 BP (Shaw and Benson 2015). At our site, the occurrence of a single (anomalous) shell layer with an upper age limit of 1300 BP suggests that the causal event was relatively uncommon, further supporting a tsunami as opposed to hurricanes, which are common in the region. Additional research is required to better characterize the tsunami history of the Gulf of Mexico, but a tsunami remains a plausible cause of the shell layer at 35 cm.
From 500 BP to present, the Ammonia-Elphidium and Agglutinated-Haynesina assemblages dominate, with a notable increase in A. tepida at 150 BP, or about 1800 CE (Fig. 4a). The only other occurrence of abundant A. tepida was during Phase 1, at 5200 BP. Thin-walled A. tepida thrive in low-salinity, protected conditions (Table 3, Murray 2006). Barrier island aggradation is unlikely to have altered the entire lagoon geometry substantially over just the last few hundred years, so the increase in A. tepida may be a response to low-salinity conditions, induced by wet climate. The nearest speleothem-derived oxygen isotope record, from Tecoh Cave (100 km west of Celestun Lagoon), shows three century-scale pulses of relatively wet conditions occurred during the Medieval Warm Period and Little Ice Age (Medina-Elizalde et al. 2010), and a wet interval is also noted in Belize speleothems (Kennett et al. 2012). A diatom record from Cenote San Jose Chulchaca, 22 km east of Celestun Lagoon indicates a decrease in salinity over the last 800 years Hodell et al. 2005). However, Whitmore et al. (1996) did not note the same trend in Cenote Sayaucil (north-central Yucatan lowlands). Furthermore, the pollen record from Lake San Jose suggests a decline of forest taxa that depend on high rainfall at the same time the cenote was freshening (Hodell et al. 2005), a trend corroborated by pollen and microfossil observations from Aguada X'caamal (75 km southeast of Celestun Lagoon) that imply drying in the northwest region during the Little Ice Age (600 to 100 BP) (Hodell et al. 2005). Together these records suggest both increased precipitation and increased evaporation for the northwest Yucatan Peninsula, which today has the highest evapotranspiration potential of the entire region (Bauer-Gottwein et al. 2011). Phase 4 of Celestun Lagoon development is thus characterized by expanding mangrove forests during a period of variable precipitation. Ammonia-Elphidium and Agglutinated-Haynesina assemblages dominate foraminifera composition of the northern semi-protected habitat with a wide range of salinity, features observed in the modern northern lagoon today.

Conclusions
This study presents the oldest paleoenvironmental record currently available for the northwest corner of the Yucatan Peninsula. Over the past 5300 years, Celestun Lagoon, on the northwest coast of the Yucatan Peninsula, has evolved from an inland mixed mangrove marsh saturated by groundwater, to a protected brackish lagoon surrounded by mangrove forest. We identified three primary controls for this change: sea level rise, coastal barrier island buildup and climate change. Geomorphic control of lagoon salinity, via reduced contribution of seawater, is a first-order control on the composition of benthic foraminifera assemblages that trend from marine to brackish taxa over time. Importantly, the 5000-year age of our record suggests barrier island aggradation, and thus lagoon formation, occurred on Holocene timescales, and led to development of these important features on the modern coastline of the Yucatan Peninsula. Thus, studies seeking to interpret causes of paleoenvironmental change on the coastal Yucatan Peninsula will need to account for changes in lagoon basin size and shape, which impact the proxy environmental records. When geomorphologic change is considered, excursions from brackish to marine taxa during lagoon development signify decreases in groundwater discharge and hence reduced precipitation over the catchment area, which is responsible for recharging the groundwater. Agreement between the Celestun sediment record and paleoclimate records from the far western and southern regions of the Gulf of Mexico suggest that Celestun Lagoon sediments reflect regional climatic events and thus regionally integrated precipitation. The inferred paleoecologic variability presented in this study, derived from processes of barrier island aggradation, sea-level rise, and climate-induced change in groundwater discharge, will contribute to future interpretations of paleoenvironmental and paleoclimate archives found along the coast of the Yucatan Peninsula.

Data Repository
All data presented in this study are available in the National Oceanic and Atmospheric Administration (NOAA) Paleoclimatology Database at the URL ncdc.noaa.gov/paleo-search/study/29171. Sediment cores are stored at the National Lacustrine Core Facility (LacCore), University of Minnesota, Minneapolis, MN, USA), and cores, core images, metadata, and bulk physical properties can be retrieved using the project code FLAM-CEL09.