Changes in Organic Matter Deposition Can Impact Benthic Marine Meiofauna in Karst Subterranean Estuaries

Subsurface mixing of seawater and terrestrial-borne meteoric waters on carbonate landscapes creates karst subterranean estuaries, an area of the coastal aquifer with poorly understood carbon cycling, ecosystem functioning, and impact on submarine groundwater discharge. Caves in karst platforms facilitate water and material exchange between the marine and terrestrial environments, and their internal sedimentation patterns document long-term environmental change. Sediment records from a flooded coastal cave in Cozumel Island (Mexico) document decreasing terrestrial organic matter (OM) deposition within the karst subterranean estuary over the last ∼1,000 years, with older sediment likely exported out of the cave by intense storm events. While stable carbon isotopic values (δ¹³C_org ranging from −22.5 to −27.1‰) and C:N ratios (ranging from 9.9 to 18.9) indicate that mangrove and other terrestrial detritus surrounding an inland sinkhole are the primarily sedimentary OM supply, an upcore decrease in bulk OM and enrichment of δ¹³C_org values are observed. These patterns suggest that a reduction in the local mangrove habitat decreased the terrestrial particulate OM input to the cave over time. The benthic foraminiferal community in basal core sediment have higher proportions of infaunal taxa (i.e., Bolivina) and Ammonia, and assemblages shift to increased miliolids and less infaunal taxa at the core-top sediment. The combined results suggest that a decrease in terrestrial OM through time had a concomitant impact on benthic meiofaunal habitats, potentially by impacting dissolved oxygen availability at the microhabitat scale or resource partitioning by foraminifera. The evidence presented here indicates that landscape and watershed level changes can impact ecosystem functioning within adjacent subterranean estuaries.

Introduction

Karst subterranean estuaries (KSEs), areas in coastal carbonate aquifers where seawater, saline groundwater, and fresh meteoric waters mix, are biogeochemical reaction zones between the land and the sea that provide habitat to characteristic biological communities. Global coastlines are ∼25% karst landscapes (Ford and Williams, 2013), which account for ∼12% of global submarine groundwater discharge (Beck et al., 2013). Dissolution and diagenesis of eogenetic carbonate landscapes create abundant cave networks and sinkholes that can become subsequently inundated through sea-level rise. This geologic setting simultaneously enhances seawater intrusion that impacts groundwater resources inland on carbonate landscapes (Perry et al., 2002; Coronado-Álvarez et al., 2011) and land-to-sea material transport through conduits and submarine springs (e.g., Fleury et al., 2007). The terms “anchialine cave” and “marine cave” have been widely applied to describe the habitat of cave-adapted biological communities living in the freshwater, mixing zone, brackish water, and also in the marine sector of karst and volcanic subterranean estuaries in the Mediterranean (e.g., Sket, 1996; Bergamin et al., 2020; Chen et al., 2020), the Tropical North Atlantic Ocean (e.g., Rubio et al., 2015; van Hengstum et al., 2019), and elsewhere (Holthuis, 1973; Stock et al., 1986; Iliffe and Kornicker, 2009; Bishop et al., 2015; Riera et al., 2018). Given modern challenges of coastal urbanization, seawater nutrient loading, and porosity of karst landscapes, there is a growing need to better understand organic matter delivery and carbon exchange between the ocean and subsurface within KSEs.

Bulk organic matter (OM) can be an important sedimentary constituent and an essential carbon source in KSEs. Direct openings to the surface facilitate particulate organic matter (POM) transport from eroding terrestrial soils and aquatic/marine primary producers into flooded coastal caves and sinkholes that are flooded by the subterranean estuary (van Hengstum et al., 2010; Gregory et al., 2017). As a result, changing terrestrial vegetation, such as wetland colonization, can impact local POM inputs to KSEs (Collins et al., 2015; Tamalavage et al., 2018). Refractory POM is either exported from the coastal aquifer or accumulates along the land-sea environmental continuum in submerged sinkholes and caves (e.g., Gabriel et al., 2009; Collins et al., 2015; Tamalavage et al., 2018; van Hengstum et al., 2020). POM distribution and accumulation in the sediment is also affected by horizontal groundwater flow through the caves (Gabriel et al., 2009; van Hengstum et al., 2010; Winkler et al., 2016). In addition, chemoautotrophic primary producers can supply in-situ POM to the benthos in the form of sinking material in marine caves (Airoldi and Cinelli, 1996). The bioavailable fraction of OM fuels biogeochemical cycling and food webs in coastal caves and sinkholes in KSEs (Pohlman et al., 1997; Socki et al., 2002; Seymour et al., 2007; Brankovits et al., 2017). However, an important question remains: Do changes in the source of sediment OM actually influence benthic meiofaunal communities in KSEs?

This study presents a 1,000-year stratigraphic record of organic matter and benthic meiofaunal community change from a flooded coastal cave on Cozumel in Mexico. Bulk OM and benthic foraminiferal community shifts were used to investigate how changes related to the sources and quantity of sedimentary OM impacted benthic meiofauna over time. The results indicate that 1) a sinkhole allows terrestrial OM to enter a subsurface stream of brackish to saline groundwater before its discharge to the coastal sea, 2) the influx of OM from mangroves in the adjacent sinkhole have decreased over the last millennium, and 3) changes in benthic foraminiferal assemblages through time are most likely related to changes in the source and decreased supply of OM.

Materials and Methods

Study Site

The island of Cozumel is located approximately 20 km from the northeastern coast of the Yucatán Peninsula in Mexico (Figure 1) and separated by the Cozumel Channel (>400 m in water depth). Cozumel and the northeastern Yucatán coast have similar Pleistocene sea-level histories and carbonate lithology (Spaw, 1978), and the flooded caves filled with speleothems (e.g., stalactites, stalagmites, columns) indicate that regional caves were in the vadose zone during previous glacial episodes (Spaw, 1978; Richards et al., 1994). The region’s tropical climate has dry (from December to April/May) and wet seasons (from May to November) (Curtis et al., 1996; Kottek et al., 2006). Like the nearby mainland, the carbonate island of Cozumel is topographically flat (mean 5 m in elevation) and the absence of fluvial systems means rainfall either directly evapotranspires or infiltrates into the unconfined aquifers (Perry et al., 2002; Beddows et al., 2007; Gondwe et al., 2010). Abundant sinkholes and well-developed cave networks impact subsurface hydrologic transport depending on seasonal rainfall and storms (Frausto-Martínez et al., 2018; Frausto-Martinez et al., 2021), which may change groundwater discharge rates and associated constituent transport to the coastal sea like on the Yucatán Peninsula (Young et al., 2008; Gonneea et al., 2014; Null et al., 2014). The region experiences microtides that span 0.1–0.2 m (Kjerfve, 1981), which also impacts daily groundwater table elevation. Tourism is driving rapid regional development and urbanization, which threatens groundwater resources (Coronado-Álvarez et al., 2011), groundwater-dependent habitats (Bauer-Gottwein et al., 2011), and other native ecosystems, like coral reefs and coastal vegetation (e.g., mangroves, seagrass) (Muckelbauer, 1990).

FIGURE 1

FIGURE 1. Cozumel Island located east of the Yucatán Peninsula, Mexico. Aerolito Cave (red star) is located on the west coast of Cozumel. Other common research sites in karst subterranean estuaries are also shown (green triangles) in Bermuda (Cresswell and van Hengstum, 2021) and the Yucatán region, such as Hoyo Negro (Kovacs et al., 2017), Carwash (Gabriel et al., 2009), Yax Chen (Collins et al., 2015), Maya Blue (Pohlman et al., 1997), and Cenote Bang (Brankovits et al., 2017).

The sinkhole (cenote) entrance to Aerolito Cave (N 20.47°, W 86.98°; 71 m by 21 m) (Figure 2) provides access to 6.1 km of mapped cave passages traversing both inland (upstream) and seaward (downstream), and red mangroves (Rhizophora mangle) are present along its northern periphery (Frontana-Uribe and Solis-Weiss, 2011). There is a ∼240 m long subsurface connection from the sinkhole to a marine bay. These marine and terrestrial openings facilitate direct detritus inputs to the cave (Mejía-Ortíz et al., 2007). These conduits hydrographically connect the sea and the aquifer as meteoric freshwater flows to the coast mixing with oceanic seawater intruding into the aquifer. Downstream passages are well-mixed and resemble typical marine conditions, whereas the water column in the open pool of the sinkhole and the upstream passages is stratified. Here, a sharp mixing zone (halocline) between brackish (∼15 psu) and saline (∼37 psu) water masses occurs typically at ∼7 m water depth (Mejía-Ortíz et al., 2007; Calderón-Gutiérrez et al., 2018). Shifts in the depth and thickness of the mixing zone in response to high-intensity rain events have been documented in Aerolito (Calderón-Gutiérrez et al., 2018), similar to the mainland Yucatán Peninsula (Coutino et al., 2017; Brankovits et al., 2018). Previous local research includes documenting endemic cave fauna (Mejía et al., 2008; Mejía-Ortíz et al., 2017; Olesen et al., 2017) and marine invertebrate ecology (Mejía-Ortíz et al., 2007; Pisanty et al., 2010; Frontana-Uribe and Solis-Weiss, 2011; Ortiz and Winfield, 2015; Márquez-Borrás et al., 2016; Calderón-Gutiérrez et al., 2018; Gómez and Calderón-Guitérrez, 2020).

FIGURE 2

FIGURE 2. Map of Aerolito Cave in Cozumel Island. Located below sea-level, the cave network promotes seaward material transport and a hydrologically interconnected environment between the sea and the aquifer, providing access to the subterranean estuary. The terrestrial sinkhole opening, colloquially known as Cenote El Aerolito, is partially surrounded by mangroves. Sediment cores (orange circles) were obtained from downstream (Core 5, Core 6, Core 7, and Core 8) and upstream (Core 1, Core 2, Core 3, and Core 4) passages. Map is adapted after Mejía-Ortíz et al. (2007).

Sediment Core Collection

Sediment cores were obtained from Aerolito using cave diving procedures in July 2014 following safety protocols outlined by the American Academy of Underwater Sciences (AAUS). Sediment cores were collected using 5 or 8 cm diameter clear polycarbonate pipes along a land-sea transect. Four cores were collected proximal to the inland cenote pool (Core 1 to Core 4, collectively referred to as the upstream sites) and four cores between the terrestrial opening and the ocean (Core 5 to Core 8, referred to as the downstream sites) (Figure 2). All coring drives ceased at a hardground that is most likely the host carbonate platform. This observation suggests that the cores sampled the complete sedimentary accumulation at each location. The water table elevation is assumed to be sea-level, given the core sites and oceanic proximity. Core collection depths are measured with respect to the modern groundwater table (∼sea-level) at the time of sampling using scuba diver depth gauges (vertical accuracy ±0.3 m). Collection depths were not corrected for local microtidal variability. The sediment water interface was located at 10.5 and 6 m below sea-level in the upstream and downstream passages, respectively (Figure 2; Supplementary Table S1). Sediment cores were then transported back to the laboratory to be split lengthwise, photographed, stratigraphically described, x-radiographed to image downcore density, and continuously stored at 4°C until further analysis.

Lithological Analysis and Radiocarbon Dating

To characterize changes in the coarse carbonate sedimentary fraction, the volumetric mass of sedimentary particles greater than 63 µm was measured in all cores using the Sieve-first Loss-on-Ignition procedure (Sieve-first LOI) (van Hengstum et al., 2016). This approach was chosen because of the known challenges associated with measuring particle size distributions of carbonate sediments using different methods (Murray, 2002). Contiguous 1-cm sediment sub-samples, with a standardized initial volume of 2.5 cm³ (1.25 cm³ for smaller diameter cores), were first sieved wet over a 63-μm mesh, dried for 12 h in an oven at 80°C, and weighed to determine the final coarse sediment residue that included both coarse organic particles and mineral matter. Thereafter, desiccated coarse sediment residues were ignited in a muffle furnace at 550°C for 4.5 h to re-mineralize coarse organic matter particles. Remaining sediment residues were then re-weighed to determine a final concentration in milligrams of coarse sedimentary particles exceeding 63 µm in diameter per unit cm³ (expressed throughout as D_>63µm mg cm⁻³). Finally, new sediment sub-samples were obtained at contiguous 1-cm intervals downcore and subjected to a Classic LOI procedure to determine variation in bulk OM content, as per standard methods (550°C for 4.5 h) and reported as weight percent (%) of the original sediment mass (Dean, 1974; Heiri et al., 2001). Uncertainty on laboratory replicates using the Classic LOI approach is typically less than ±1%.

Age constraint was provided by radiocarbon dating bulk sedimentary OM particles and, when available, terrestrial plant macrofossils (e.g., twigs). Biogenic carbonates were deemed inappropriate for age purposes because hardwater and marine reservoir effects create problems for calibrating conventional radiocarbon results from marine invertebrates into sidereal years. Six radiocarbon dates were obtained across four cores: Core 1 (2 bulk sedimentary OM), Core 5 (2 terrestrial plant macrofossils), and Core 7 (2 terrestrial plant macrofossils; Supplementary Table S2). Material sampled for radiocarbon dating was rinsed with deionized water, dried, weighed, and then submitted for radiocarbon determination at the National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) Facility at Woods Hole Oceanographic Institution. Conventional radiocarbon dates were calibrated to calendar years before present (cal yrs BP) using IntCal20-NH (Reimer et al., 2020) calibration curves using CALIB 8.2 (Stuiver et al., 2021) (Supplementary Table S2).

Stable Carbon Isotope and C:N Ratio Analysis of Organic Matter

Seven cores (Core 1 to Core 7) were sampled contiguously downcore to determine the stable carbon isotopic ratio of bulk organic matter (δ¹³C_org) and the C:N ratio. In coastal settings, δ¹³C_org and C:N measured from sediment cores can diagnose the organic matter sources deposited through time (e.g., Gonneea et al., 2004; Lamb et al., 2006; Tamalavage et al., 2018; Gonneea et al., 2019). For our cores, a 5 mm sub-sample was obtained every 1 or 2 cm downcore for the entire length (n = 222). Approximately 1.25 cm³ of sediment was directly acidified with ∼8 ml of 1 M hydrochloric acid for 24 h, or until effervescence ceased, then desiccated at 50°C and homogenized into a powder. Determination of δ¹³C_org and C:N ratio occurred at the Texas A&M University Stable Isotope Geoscience Facility using an Erba NA 1500 Elemental Analyzer that is coupled to a Thermo Scientific DELTA^plus isotope ratio mass spectrometer. Final isotopic ratios are reported relative to the standard Vienna Pee Dee Belemnite (VPDB) for carbon in per mil (‰). Analytical precision for two international standards, USGS 40 and 41a which represent a range of δ¹³C values for a two-point calibration, are within 0.2‰ (1σ) δ¹³C. Additional internal standards were used as an internal check for reproducibility. The residual organic carbon %C (acidified residue) was divided by %N (acidified residue) to calculate C:N ratio.

A simple two-source mixing model using δ¹³C_org values (Fry, 2007) calculated an estimate for the relative terrigenous versus marine supply to sedimentary OM through the terrestrial versus marine openings. This calculation was based on the assumptions that: 1) OM in the cave sediments is derived from two main sources (marine-derived OM transported landward through the marine cave opening and terrestrial-derived OM eroding into the cave through interconnected sinkholes) and 2) the endmember values are preserved with no or minimal changes after burial. We recognize that a two-source mixing model grossly oversimplifies the organic carbon sources delivered to coastal settings because it does not account for the complexity of variable OM sources and diagenesis. For example, potential sedimentary OM sources in coastal settings include primary producers, C3 and C4 plants of coastal vegetation, terrestrial vegetation, and marine and freshwater algae (e.g., Fourqurean et al., 2005; Lamb et al., 2006). Moreover, diagenesis (i.e., biological, chemical, and physical processes) can affect sedimentary OM composition, including δ¹³C_org values, after burial (e.g., Fogel et al., 1989; Benner et al., 1991; Gonneea et al., 2004). Despite the complexities, this approach can provide a first-order estimate for the relative contribution of terrigenous versus marine sources to the sediment OM in a dark cave setting that facilitates carbon exchange between the land and the sea.

Considering the habitats adjacent to Aerolito, we assume that terrestrially derived OM primarily includes mangrove-derived detritus, other terrestrial and wetland plants, and, in this case, even some biomass from aquatic primary producers in the sinkhole. The assumption that sediment OM was a mix of marine and terrigenous sources can be expressed mathematically by:

${\mathrm{f}}_{terrestrial}+{\mathrm{f}}_{marine}=1,(1)$

where f_terrestrial and f_marine are the fractions of OM from the two different sources. To calculate f_terrestrial, the following equation (Eq. 2) was used:

${\mathrm{f}}_{terrestrial}=\frac{{\delta }_{se\mathit{dime}nt}-{\delta }_{marine}}{{\delta }_{terrestrial}-{\delta }_{marine}},(2)$

where δ_sediment is the measured δ¹³C_org value of the sediment sample, δ_terrestrial is the average δ¹³C_org value measured across a variety of mangrove plants from multiple habitats in the Yucatán Peninsula (−28.24‰; Gonneea et al., 2004), δ_marine is the average δ¹³C_org value of POM measured from the sea along the coastline in the nearby Cozumel Channel (−20.1‰; Brankovits et al., 2017) which is assumed to represent the isotopic content of marine algae and other components of the coastal marine POM pool. The determined fraction contribution of a source will only reflect the component of the OM that is preserved and not the original flux to the sediment.

Benthic Foraminifera

Subfossil benthic foraminifera can provide a physical record of changing benthic environmental conditions, such as changes in groundwater salinity, dissolved oxygen concentrations, or organic matter supply (Caralp, 1989; Kaiho, 1994; Gooday, 2003; Jorissen et al., 2007). As such, benthic foraminifera were the model organism used to provide a proxy for any long-term changes in the benthic meiofaunal community in the passages proximal to the ocean of Aerolito Cave. After an assessment of the radiocarbon results (discussed further below), our research design included: 1) sampling both the oldest (basal core samples) and youngest (core-tops) parts of the stratigraphic record (n = 13), and 2) comparing Aerolito foraminiferal assemblages to the dataset of recent benthic foraminifera in Bermudian underwater caves (Cresswell and van Hengstum, 2021; n = 125).

Foraminifera were sampled from Core 1, Core 2, Core 4, Core 5, and Core 7 (n = 13, Figure 3). Core-top samples (n = 5) were collected from 5 to 7 cm core depths to ensure that both epifaunal (living at or near the sediment water interface) and infaunal (living deeper in the sediment) foraminiferal taxa would be represented, along with any taxonomic smoothing imparted by local taphonomic effects (e.g., currents or any bioturbation in oxygenated marine settings) or patchiness in local benthic foraminiferal communities. Samples were also collected near the base of the selected cores (n = 8). For each sample, a 1-cm wide interval of 0.63 cm³ of sediment was obtained from the core and rinsed over a >63-μm sieve to concentrate benthic foraminiferal tests. Benthic foraminifera were counted with a minimum census of 300 individuals per sample (Patterson and Fishbein, 1989; Fatela and Taborda, 2002; Schönfeld et al., 2012). In one sample, 0.63 cm³ of sediment was insufficient to reach our minimal census count goal and additional sediment was sieved from the same stratigraphic level, resulting in sample volumes spanning 0.63–3.13 cm³. To promote counting efficiency, four samples were split into representative aliquots using a wet-splitter prior to census counting (Scott and Hermelin, 1993). Foraminifera were wet counted and identified to the species level, with taxonomic identification based on taxonomic literature (Cushman, 1933; Loeblich and Tappan, 1987) and previous studies on cave foraminifera in Bermuda (van Hengstum and Scott, 2011; Little and van Hengstum, 2019; Cresswell and van Hengstum, 2021). The raw foraminiferal database of census counts from the 13 samples from Aerolito included 40 taxonomic units identified to the species level from 60 genera. The resulting dataset was subjected to both unconstrained Q-mode and R-mode cluster analysis to observe any natural groupings between the samples using R Studio software (version 1.1.423; R Core Team, 2020). Using a similarity profile analysis (SIMPROF), 1000 permutations of the cluster analysis were performed to yield a dendrogram with significant clusters (p < 0.05). First, the proportional abundance of each taxa was calculated using the “decostand” function in the “vegan” package (Oksanen et al., 2019) and bubble plots were generated using “ggplot2” (Wickham, 2016). Cluster diagrams and SIMPROF were produced in the “vegan” and “clustsig” (Whitaker and Christman, 2014) packages, using the proportional abundance values, Bray-Curtis dissimilarity index, and UPGMA cluster algorithms. Taxa were excluded from cluster analyses if the calculated standard error was greater than the proportional abundance for a taxonomic unit in all samples.

FIGURE 3

FIGURE 3. Detailed photographs and x-radiographs of sediment cores are shown with downcore coarse-grained particle deposition (Sieve-first LOI, D_>63 μm), bulk organic matter (Classic LOI, weight %), stable carbon isotopic values (δ¹³C_org), and C:N ratios. Radiocarbon and microfossil sample positions are marked on relevant cores.

The foraminiferal results from Aerolito (n = 13) were then compared to dataset of recent foraminiferal distributions in surface samples from Bermudian marine caves (n = 125, Cresswell and van Hengstum, 2021) to create a new composite dataset of 138 samples. Taxonomic differences were evaluated at the genus level in the composite dataset to 1) promote equitable comparison between regions, and 2) taxonomic bias imparted by separate microscopists. While the genera Patellina and Patellinoides were originally differentiated in the samples from Aerolito Cave, they were subsequently grouped into a pseudo-Patellina genus for comparison to Bermuda (see discussion in Cresswell and van Hengstum, 2021). Biodiversity indices [Shannon-Weiner Index (H’) and Fisher’s Alpha (α)] for each sample was calculated from the raw counts in the final database using the “vegan” package. Species richness (S_rare) was calculated using rarefaction with the “rarefy” function in “vegan” because of the spread in the number of individuals counted between samples (ranging from 377 to 1,151). The sub-sample size was set to 377 individuals, which is the minimum sample size observed. Finally, ternary diagrams were generated using the “ternary” package in R (Smith, 2017) with respect to: 1) the proportions of wall structure and life mode (epifaunal versus infaunal), and 2) the δ¹³C_org value of bulk OM. These ternary plots compared the total percent miliolid taxa (epifaunal and infaunal) to the remaining epifaunal (hyaline and agglutinated) and infaunal taxa (hyaline and agglutinated). Since agglutinated taxa were not common in the dataset from these shallow, tropical carbonate settings, this ternary division best highlighted structural differences in the miliolid-dominated dataset. Designation of hyaline and agglutinated taxa as epifaunal or infaunal primarily followed: Corliss and Chen (1988), Corliss and Emerson (1990), Corliss (1991), Jorissen and Wittling (1999), Murray (2003), with additional references to Rathburn and Corliss (1994), Alve and Bernhard (1995), Gooday (1996), Bernhard and Bowser (1999), Fontanier et al. (2002), and van Hengstum and Scott (2012).

Results

Lithology and Chronology

Sediment in Aerolito is dominated by carbonate materials and organic matter, with negligible contributions from siliciclastic minerals and mean sediment accumulation is just 107 ± 19 cm (eight cores, range: 77–135 cm, Supplementary Table S1). Qualitatively, the sediment presents as fine carbonate matrix (e.g., micrite), with inclusions of coarse-grained shell material (bivalves and gastropods), carbonate fecal pellets, and sponge spicules that could be readily observed through stereomicroscopy of wet-sieved sediment (over a >63-μm mesh). In general, basal core sediment was dark brown in color versus the core-tops, and orange-colored iron-rich sediment appears at the base of Core 2 and Core 4, and from 25–17 cm in Core 1 (Figure 3). This sedimentary unit is separated with a sharp interface from overlying carbonate sediments and has decreased bulk organic matter content. Lastly, sedimentary bedding in the carbonate mud was not immediately visible in photographs, but linear bedding could be observed in x-radiographs (Figure 3). These observations indicate minimal (<3 cm) vertical post-depositional sediment mixing from bioturbation. However, bedding structures cannot be assumed to indicate long periods of continuous sedimentation because similar sedimentary structures can be rapidly generated following current dissipation and re-dissipation of suspended sediment load after an intense storm (e.g., hurricane).

The downstream passages have generally more bulk OM than the upstream passages (Figure 3). Downstream sites have highest bulk OM content toward the base of the core (up to 17.7%) with a gradual upcore decrease at an average rate of −0.11% OM cm⁻¹. This upcore decrease in bulk OM is also present in upstream core sites (Core 2, Core 4, above 17 cm in Core 1), but with a lower rate of −0.04% OM cm⁻¹. The lightest-shaded sedimentary layers in the greyscale radiographs document the densest layers, which typically correspond to iron-rich sediments or relatively higher fractional abundance of coarse sediment grains (D_>63µm mg cm⁻³) that have the lowest measured bulk OM content (as low as 7.8% OM in Core 4). For instance, the coarser-grained iron oxide rich layer from 16 to 23 cm in Core 1 has a mean D_>63µm of 100 mg cm⁻³ ± 40, whereas the rest of the 22 cm core has a mean D_>63µm is 28 mg cm⁻³ ± 15. Coarse sediment mass generally increases upcore in the downstream cores, but there are no consistent trends in the upstream section (Figure 3). Spatially, the downstream cave areas between the sinkhole and the ocean are typically coarser than the upstream, with mean D_>63µm values of 116 mg cm⁻³ ± 79 (n = 163 cm) versus 67 mg cm⁻³ ± 47 (n = 178 cm) upstream. Sediment with higher bulk OM that is typically near the base of the cores is hereafter referred to as OM-enriched deposits, whereas sediment with lower bulk OM is hereon referred to as OM-limited deposits.

Six radiocarbon dates broadly constrain the onset of sedimentation (Figure 3; Supplementary Table S2). The basal date from Core 1 has a probable 1σ age of 940 ± 10 cal yrs BP at 27–28 cm, which represents the oldest reliable date for the minimum onset of preserved sediment accumulation within Aerolito. The second date from Core 1 that is higher in the succession at 19–20 cm dates older to 2325 ± 20 cal yrs BP. However, this date was pulled from an interval with relatively high coarse-grained sediment content (>63 µm), which is perhaps dominated by older sediments reworked by a high energy hydrologic event and were subsequently incorporated into the radiocarbon dated bulk OM. In the downstream passage at Core 5, a date at 36–37 cm is aged to 2450 ± 50 cal yrs BP, but a date just 2 cm higher at 34–35 cm is 660 ± 10 cal yrs BP when bulk OM content increases to ∼18%. If the deeper date from Core 5 is accurate, it would suggest ∼1,800 years of non-deposition at this core site. Sediments from the basal 35–43 cm interval in Core 5 have relatively low bulk OM (12.9% ± 1.1), which is not represented in other downstream cores. Alternatively, the date at 36–37 cm in Core 5 was obtained on an older plant fragment that became deposited in the cave long after the plant died. The dated plant fragment deeper in Core 7 at 37–38.5 cm is 580 ± 20 cal yrs BP, but the dated plant fragment higher in the succession at 33.5–34 cm has an age of 740 ± 10 cal yrs BP. Again, this shallower and older date is likely a reworked plant fragment from the terrestrial surface that was deposited long after the original plant died. Based on the youngest possible dates in Core 1 (940 ± 10 cal yrs BP, 27–28 cm), Core 5 (660 ± 10 cal yrs BP, 34–35 cm), and Core 7 (580 ± 20 cal yrs BP, 37–38.5 cm), the sediment sampled by the cores accumulated in the last millennium, which is long after the cave was first flooded by Holocene sea-level rise.

Stable Carbon Isotopic Values and C:N Ratios

The δ¹³C_org values across all core samples (n = 222) range from −22.5 to −27.1‰, and the C:N values vary from 9.9 to 18.9. There is a general upcore decrease in C:N ratios, and δ¹³C_org values become more enriched (i.e., positive) (Figure 3). As a result, the OM-limited sedimentary sections near the top are relatively ¹³C-enriched and lower in their C:N values compared to the core bottoms. Spatially, sediment archived in downstream passages are characterized with relatively depleted δ¹³C_org values (mean: −25.0 ± 0.8‰, ranging from −23.6 to −27.1‰) and high C:N ratios (mean: 14.1 ± 1.6, ranging from 10.2 to 18.9). In contrast, sediments located upstream from the terrestrial opening show more enriched δ¹³C_org values (mean: −23.8 ± 0.9‰, ranging from −22.5 to −26.2‰) and lower C:N ratios (mean: 12.0 ± 1.0, ranging from 9.9 to 14.8).

In general, the two-source mixing model based on δ¹³C_org values indicates that ¹³C-depleted OM sources (e.g., terrestrial plants, mangroves), deliver a considerable portion of the organic carbon to the sediment (dataset range: 29–85%). However, there is generally greater terrigenous contribution at the base rather than at the core-tops (Figure 4). These patterns suggest that terrestrial organic carbon delivery to the sediment was greater in the past (basal core sediment) than in the present (core-top). Comparing the measured δ¹³C_org and C:N values with previously reported values of potential OM sources (Lamb et al., 2006) further indicates that organic carbon preserved in the Aerolito sediment is derived from two primary sources: vascular C3 plants and marine POM/algae (Figure 5A). In comparison to similar coastal caves in Bermuda (Cresswell and van Hengstum, 2021), sediment in Aerolito preserved less marine-sourced organic carbon, like from marine algae (Figure 5A). Interestingly, sedimentary OM closer to the marine opening (downstream) received greater contributions from vascular plant detritus than those near the terrestrial opening (upstream) (Figure 5). Moreover, higher bulk OM accumulation (up to 17.7% OM) was associated with greater contribution of percent terrestrial organic carbon (up to 85%) (Figure 5B).

FIGURE 4

FIGURE 4. Isotope mixing model results. Downcore trends in the percent terrestrial organic carbon (OC) as calculated from the δ¹³C_org values using a simple two-source mixing model. The terrestrial endmember value is δ¹³C_org = −28.2‰, as the mangrove vegetation reported by Gonneea et al. (2004), while the marine endmember is δ¹³C_org = −20.1‰, as the marine POM along the coast of the Cozumel Channel reported by Brankovits et al. (2017).

FIGURE 5

FIGURE 5. (A) Cross-plot between stable carbon isotopic values (δ¹³C_org) and C:N ratio of bulk organic matter in Aerolito Cave (this work) and Bermuda caves (Cresswell and van Hengstum, 2021). (B) Cross-plot between bulk organic matter content (Classic LOI, weight %) and the percent marine/terrestrial organic carbon as calculated from the two-endmember mixing model based on δ¹³C_org values.

Foraminiferal Communities

Benthic foraminifera were abundant in sediment from sampled areas in Aerolito (Figure 2). Two groupings can be identified in the dendrogram produced from unconstrained Q-mode cluster analysis (Figure 6), which generally separates core-top samples (i.e., less bulk OM, more depleted δ¹³C_org values) from the bottom of the cores (i.e., higher bulk OM, more enriched δ¹³C_org values). The two communities are hereafter referred to as the OM-limited Community (recent) and the OM-enriched Community (past) (Figure 6). There is little taxonomic difference between the two communities identified by Q-mode cluster analysis, which indicates negligible long-term changes in water column structure in the sampled areas.

FIGURE 6

FIGURE 6. Dendrogram produced from two-way hierarchical cluster analysis of the benthic foraminiferal data in Aerolito. The Q-mode cluster analysis identifies two main groups: OM-limited Community (purple field) and OM-enriched Community (brown). The R-mode cluster analysis highlights dominant (Group I), common (Group II), rare (Group IV), and very rare taxa (Group IV). Sample codes (bottom row) are based on the core ID followed by sample depth in cm (e.g., C7_5.5 was sampled from Core 7 at 5.5 cm core depth).

However, the relative abundance of infaunal (i.e., Bolivina) versus miliolid taxa differentiates the two foraminiferal communities, along with the proportions of shared common taxa (Figure 6). In the OM-limited Community (core-tops), the dominant fauna (>10%) are the porcelaneous miliolid taxa (33.1%; e.g., Quinqueloculina, Triloculina, Spirophthalmidium) and the infaunal hyaline genus Bolivina (14.1%). Other common taxa (>5%) include Rosalina spp. (mean: 10.0%) and Spirillina (mean: 5.7%) (Supplementary Table S3). In contrast, the OM-enriched Community (basal sediments) has an increased abundance of Bolivina (20.7%), decreased miliolids (21.2%), and common occurrences of Rosalina (9.0%) and Eponides (6.1%). Foraminiferal density is also higher in the OM-limited Community (mean: 4357 individuals cm⁻³) than the OM-enriched Community (mean: 915 individuals cm⁻³), though both communities have a similar diversity (H’_OML: 2.36; H’_OME: 2.62).

The R-mode cluster analysis identifies four groups of species (Figure 6). Group I contains the top four dominant taxa (mean: ≥5%) present in both communities from the unconstrained Q-mode cluster analysis: miliolid taxa, Bolivina, Rosalina, and S. vivipara. Group II comprises the greatest number of taxonomic units (n = 19), which are common to both communities and includes taxa such as the infaunal taxa Bulimina, Ammonia, and rare agglutinated fauna. Ammonia is rare in the OM-limited Community (mean: <1%), but it achieves a higher proportion in the OM-enriched Community (mean: 4.7%). Group III (n = 10) contains rare taxa (mean: <2%) that are present in nearly every sample, whereas Group IV consists of only two genera which occur in only two samples (mean: <1%).

The two communities also exhibit differences in foraminiferal test structure and morphology (Figure 7; Supplementary Table S4). The ternary diagrams depict an overall increase in miliolid taxa in the core-top (i.e., youngest) versus bottom (i.e., oldest), and differences in the proportions of infaunal taxa versus epifaunal taxa. The older OM-enriched Community (numbered brown data points on Figure 7A) has more infaunal taxa (mean: 46.5%; range: 42.5–52.1%), such as Bolivina, Buliminella, and Fursenkoina, with decreased miliolid abundance (20.7%; range: 17.3–23.7%) and epifaunal taxa (mean: 31%; range: 27.6–38.4%). The younger OM-limited Community (numbered purple data points on Figure 7A) has higher miliolid abundances (mean: 20.7%; range: 17.3–23.7%), epifaunal species (mean: 28.7%; range: 23.8–38.4%), and decreased infaunal taxa (mean: 34.4%; range: 26.1–43.5%).

FIGURE 7

FIGURE 7. Ternary diagrams depicting the proportion of benthic foraminiferal taxa in each sample organized into miliolid, non-miliolid epifaunal, and non-miliolid infaunal taxa. Results from this study (numbered circles and squares) are compared to benthic foraminifera from karst subterranean estuaries in Bermuda (smaller circles) (Cresswell and van Hengstum, 2021). Samples in the ternary diagram are color coded (A) by community designation (Aerolito Cave OM-limited Community in purple and OM-enriched Community in brown), and (B) by δ¹³C_org values of sediment OM.

Discussion

Integrity of Sedimentation

Based on the radiocarbon ages, the basal age result of ∼940 cal yrs BP from Core 1 most likely constrains the earliest possible onset of modern sedimentation. An overall young (i.e., last millennium at most) age for the recovered stratigraphy from this part of Aerolito Cave is also supported by the young basal age from Core 7 of 580 ± 20 cal yrs BP at 37–38.5 cm. Parts of the western North Atlantic margin experienced increased intense hurricane activity from 2600 to 1000 cal yrs BP (Donnelly and Woodruff, 2007; van Hengstum et al., 2016). Thus, it is possible that any older than ∼1,000 years old sedimentary deposits were washed out from the cave when more frequent intense hurricane activity repetitively increased the velocity of groundwater flow.

Although an assessment of temporal change is tempered by potential impacts to the stratigraphic integrity by more recent high-energy storm events, the recovered cores do preserve a record of long-term environmental change in Aerolito during the last millennium. Older radiocarbon dated material (e.g., ∼2450 and 2325 cal yrs BP in Core 5 and Core 7) could either have been transported and deposited into the cave after death of the original plant, or potentially indicate re-suspension of preexisting sedimentary material in the cave following a hurricane. Occasional increase in coarse-grained sediment content in the cores collectively suggest that modern sedimentation has likely been affected by episodic hydrologic events like hurricanes. However, distinctive coarse-grained tempestite beds were not visible in the stratigraphy despite the available sediment supply (e.g., shell material), and the upcore stratigraphic trends are not characterized by repetitive fining upward sequences that commonly develop from hurricane-induced suspension and redeposition of sediment in marine environments (Kuenen and Migliorini, 1950; Wanless, 1981; Bentley et al., 2002; Toomey et al., 2013). The original horizontal bedding imaged by the x-radiographs and the intact iron-rich sedimentary deposit in Core 1, with sharp contacts with overlying and underlying carbonate sediment, all lend further support that the sediments can be used to evaluate long-term environmental trends (Figure 3). While some bioturbation is expected, these sedimentary structures indicate that complete vertical post-depositional sediment mixing and homogenization from bioturbation or groundwater currents have not occurred. Taken together, these observations lend confidence that the recovered cores preserve, at least, a low-resolution record of long-term environmental change in Aerolito.

Iron Oxide Sediment

Sedimentary layers with increased iron oxide concentrations occur both intercalated within a core (Core 1) and at the base of cores (Core 2 and Core 4; Figure 3). The intercalated iron-rich sedimentary unit in Core 1, which has decreased bulk OM content, is separated with a sharp interface from overlying sediments. It is highly likely that Cozumel experiences some Saharan dust deposition from atmospheric transport (Prospero and Nees, 1986; Prospero and Lamb, 2003; Muhs et al., 2007). However, the recovered iron-rich sediments are likely not derived from Saharan dust deposition because these sediments are most likely less than 1,000 years old (see above), and Aerolito and other caves flanking the Cozumel Channel at similar elevations would have already been flooded by Holocene sea-level rise (Gabriel et al., 2009; Milne and Peros, 2013; Khan et al., 2019).

An alternative hypothesis is that these iron-rich sedimentary deposits have an in situ, subaquatic origin. In siliciclastic subterranean estuaries, the oxidative precipitation of dissolved Fe(II) from the mixing of seawater and groundwater causes iron oxide deposition in the subsurface (Charette and Sholkovitz, 2002). Similar iron-rich deposits have been observed in Palm Cave (van Hengstum et al., 2019) and Green Bay Cave (van Hengstum et al., 2011) in Bermuda, where it has been hypothesized that they originate from the oxidative precipitation of Fe(II)-rich minerals in an “iron curtain” along a redox gradient at the sediment-water interface as anoxic saline groundwater upwells into the oxygenated caves within carbonate platforms (van Hengstum et al., 2011; van Hengstum et al., 2019). We speculate that a similar mechanism is also occurring in Aerolito Cave, but further work is required to assess these geochemical processes in KSEs.

Evidence for Reduction in Organic Matter Inputs From a Mangrove Habitat

Based on the simplistic two-endmember mixing model, the terrestrial-sourced OM contributions to the sediment OM pool vary from 29 to 85% (Figure 5B). Although the proportion of ¹³C-depleted OM has decreased over time, its sources continued to contribute to the sediments throughout the record (Figure 4). These patterns are most likely a result of mixed OM supply from the mangroves or other vegetation surrounding the terrestrial opening, similar to observations in the Yucatán Peninsula (van Hengstum et al., 2010; Collins et al., 2015), The Bahamas (Tamalavage et al., 2018; Björnerås et al., 2020), and Bermuda (van Hengstum et al., 2011; Little and van Hengstum, 2019). Globally, mangrove ecosystems are important carbon reservoirs (Donato et al., 2011) and a major source of organic carbon to the oceans (Dittmar et al., 2006; Maher et al., 2013), and they clearly contribute considerable organic carbon to tropical KSEs. POM with low C:N ratios and ¹³C-depleted isotopic values, which is indicative of in-situ biomass production by freshwater algae in the flooded sinkhole, is another likely source of OM that fluxes into the cave environment (Tamalavage et al., 2018). It is possible that the higher C:N ratios in basal core sediment is partly from diagenetic alteration of OM from the preferential loss of nitrogen (e.g., Fogel et al., 1989; Benner et al., 1991; Gonneea et al., 2004; Lamb et al., 2006). However, increased C:N values in the basal sediments are associated with high bulk OM content and depleted δ¹³C_org values, which suggests that OM inputs from the terrestrial sinkhole were higher in the past than in the present.

An upcore decrease in bulk OM is present in all cores and is associated with an overall decline in C:N ratios, increasing δ¹³C_org values, and consequent decrease in terrestrially derived OM contribution as estimated by the two-source mixing model (Figure 4). Because the joint shifts in bulk OM, C:N ratios, and δ¹³C_org values are most likely a result of lessening OM supply by mangrove vegetation through the sinkhole pool, it is possible that the mangrove habitat on the adjacent surface have shrunk in the last ∼1,000 years (Figure 8). One may consider that increased POM inputs from the marine algae (with more enriched δ¹³C_org values) could also explain the decrease in the fraction of terrestrially derived OM in the sediments. However, the OM-enriched sediment in Aerolito has depleted δ¹³C_org values and relatively high C:N ratios, which suggests that the observed shift in OM accumulation was primarily driven by changing supplies from wetland and/or terrestrial plants. Shifts in OM supply may affect biogeochemical cycling in KSEs, which could impact the composition of the discharged groundwater (Moore, 1999; Moore, 2010), but how this change in OM supply impacted elemental cycling in Aerolito through time remains unknown.

FIGURE 8

FIGURE 8. Landscape-level conceptual model summarizing organic matter (OM) dynamics in the karst subterranean estuary (KSE). Hydrologic forcing (blue arrows) through cave conduits accelerates the transport and exchange of material between the terrestrial and marine environments. These processes are especially amplified during episodic events (e.g., hurricanes). OM inputs from mangroves and terrestrial vegetation into the KSE were greater in the past (gray arrow) than in the present (black arrow), resulting in a decrease in the accumulated OM in the caves over time. Reduced OM inputs likely caused a shift in biogeochemical conditions which prompted a change in the benthic foraminifera communities. Higher terrestrial OM accumulation near the marine opening is from persistent seaward material transport during the last millennium.

Bulk OM content is generally higher near the marine opening that is downstream from the sinkhole, versus core sites that are proximal to and upstream from the sinkhole. Higher bulk OM has also been documented downstream from inland sinkholes (i.e., cenotes) on the Yucatán Peninsula (van Hengstum et al., 2009; Collins et al., 2015). Higher OM content downstream is associated with high C:N and depleted δ¹³C_org values (Figure 3) and, thus, likely receives larger contributions from terrestrial sources (Figures 4, 5B). Tidal-pumping of seawater between the ocean and karst aquifer facilitates the inland transport of marine-derived OM through the caves up to 400 m from the marine opening (Figure 4). At the same time, the general land-to-sea direction of fresh groundwater flow in the region promotes the seaward transport of non-marine POM delivered to the flooded cave environment through the terrestrial opening, rather than upstream from the terrestrial opening. This scenario is consistent with our observations in situ. During sample collections, sediment accumulation was negligible in the conduits further inland from the upstream sampling locations. Extreme hydrologic events (e.g., hurricanes) can amplify the seaward groundwater transport and subsurface currents, which may further influence POM distribution in the conduits. Hydrodynamic sorting may also contribute to the dissemination of terrigenous POM during hydrologic transport (Keil et al., 1994; Goñi et al., 1997).

In summary, evidence for continuous accumulation of terrestrially derived OM in the downstream passages confirms that the conduits have facilitated land-to-sea transport of OM during the last millennium (Figure 8). These observations support findings in other KSEs that the type of surface vegetation and the presence of terrestrial openings are key factors regulating POM inputs to coastal sinkholes (e.g., van Hengstum et al., 2010; van Hengstum et al., 2011; Collins et al., 2015; Gregory et al., 2017; Tamalavage et al., 2018; Little and van Hengstum, 2019). Such high POM inputs can promote eutrophic KSE sub-habitats that are in contrast with oligotrophic cave areas distant from openings (Pohlman, 2011; Brankovits et al., 2017) and the greater aquifer (groundwater-rock matrix) with limited connectivity to the surface. Open sinkholes, however, are known to exert a clear biogeochemical influence on the subterranean realm on tropical carbonate landscapes, and their influence is transmitted to oligotrophic caves and to the ocean by groundwater flow. For example, quantitative evidence from time-series data have shown the importance of rainfall-induced hydrological transport of OM from inland mangroves to the coastline through cave networks (Coutino et al., 2020). Although karst platforms account for a substantial percentage of certain trace element inputs to the ocean globally (Beck et al., 2013; Gonneea et al., 2014; Mayfield et al., 2021), point-sources of groundwater discharge only represents a (likely small) portion of these inputs. Point-source discharge, however, can influence coral reefs (Crook et al., 2012) and other local ecosystems (Moosdorf et al., 2015) while contributing to the heterogeneity of land-to-sea material fluxes along global coastlines (Luijendijk et al., 2020; Moore and Joye, 2021; Rocha et al., 2021).

Broader Comparison of Benthic Foraminifera Communities

Aerolito Cave contains benthic foraminiferal taxa that are similar to flooded coastal caves in Bermuda (Cresswell et al., 2017; Cresswell and van Hengstum, 2021). However, the foraminiferal communities of Aerolito have a lower proportion of epifaunal taxa and increased infaunal taxa. For example, non-miliolid epifaunal taxa reach proportions as high as >80% in Bermuda while the maximum proportion recorded in Aerolito is only 38.4% (Figure 7). The OM-enriched Community of Aerolito is most similar to foraminiferal communities from Bermuda that are also found in sediments with depleted δ¹³C_org values (likely because of increased terrestrial OM input), such as Community 4 from Green Bay Cave (GBC), as both contain a larger proportion of infaunal taxa (mean_OME: 46.5%; mean_GBC: 33.7%) and decreased miliolids (mean_OME: 20.7%; mean_GBC: 12%) (Figure 7). Aerolito’s OM-limited Community (core-top) is more similar to Community 5B in Green Bay Cave, as both foraminifera communities have similar abundances of miliolid taxa (mean_OML: 35.6%; mean_GBC: 36.2%). Green Bay Cave’s well-oxygenated water column is exposed to regular tidally forced groundwater-seawater exchange.

Overall sediment δ¹³C_org values in Aerolito are more depleted than observed in Bermuda marine caves (Figures 5A, 7B), which suggests that marine POM provides a larger contribution to the sedimentary OM pool in Bermuda KSEs than in Aerolito. Like Bermuda, sediments with a higher proportion of miliolid taxa in Aerolito tend to have more enriched δ¹³C_org values (Figure 7B), which suggests that marine-derived POM delivery is beneficial to miliolid taxa in these environments. Sediment δ¹³C_org values in Aerolito are most similar to the terrestrially influenced regions of Green Bay Cave (mean: −24.37‰) and Cow Cave (mean: −25.93‰) in Bermuda. More broadly, the presence of Spirophthalmidium emaciatum is also consistent with marine caves in Spain (Bergamin et al., 2020) and Bermuda (Cresswell and van Hengstum, 2021). Future comparisons to marine caves elsewhere in the tropical North Atlantic, such as Giant Cave in Caye Caulker (Fosshagen and Iliffe, 1991), could provide further insight into the differential impacts of OM delivery on cave benthic meiofaunal communities. It remains unknown how variation in the relative contributions from marine and terrigenous OM sources cascades and impacts the subterranean food webs in Bermuda versus Aerolito, but perhaps these differences contribute to regional biogeographical variances in higher-order cave fauna.

Benthic Foraminifera Respond to Organic Matter Changes

Foraminiferal community shifts suggest that decreasing OM inputs from a mangrove habitat over time impacted the benthic meiofaunal community in Aerolito Cave. Benthic foraminiferal communities are significantly controlled by OM quantity and source (e.g., terrestrial versus marine), both as a food source and by impacting bottom water and pore-water oxygenation (Corliss and Chen, 1988; Gooday and Turley, 1990; Bernhard, 1992; Jorissen et al., 1995; Belanger et al., 2020; Dimiza et al., 2020). The environmental control of OM supply on benthic foraminiferal communities has been documented in a variety of coastal environments such as the Mediterranean (Dimiza et al., 2020) and the Bay of Biscay (Fontanier et al., 2002), as well as the deep sea (Corliss and Chen, 1988; Jorissen et al., 1995). OM effects on benthic foraminiferal microhabitats have been previously explored within the KSEs of Bermuda (van Hengstum and Scott, 2011; Cresswell et al., 2017; Little and van Hengstum, 2019; Cresswell and van Hengstum, 2021), where habitats with more increased terrestrial-based OM (depleted δ¹³C_org values) are favored by certain taxa (e.g., Helenina, Bolivina, Rosalina).

In Aerolito, the basal part of the cores that have increased bulk OM and depleted δ¹³C_org values are associated with increased infaunal benthic foraminifera (e.g., Bolivina, Bulimina, Fursenkoina) and Ammonia beccarii (mean: 4.7% versus <1% in core-tops, Figures 6, 7). However, there was likely no widespread or long-term change in the physicochemical characteristics of the water column (e.g., dissolved oxygen and salinity). While increased Ammonia and infaunal taxa can be related to decreased dissolved oxygen concentration (Moodley and Hess, 1992; Jorissen et al., 1995; Sen Gupta et al., 1996), miliolid foraminifera are generally intolerant to long-term sub-oxic to anoxic conditions (<1.5 ml/L dissolved oxygen; Kaiho, 1994) and they are abundant in all Aerolito samples examined. In addition, the lack of agglutinated taxa (e.g., Trochammina, Entzia, Miliammina) associated with brackish water in caves and sinkholes (van Hengstum et al., 2008; Little and van Hengstum, 2019; Cresswell and van Hengstum, 2021) indicates that brackish conditions at the sediment-water interface were not persistent over the last 1,000 years in Aerolito. These results suggest that despite observations for short-term water column changes in Aerolito from increased freshwater discharge during hurricanes (Calderón-Gutiérrez et al., 2018), the conditions rebound back to the pre-storm conditions fairly quickly. Alternatively, subtle changes in OM sourcing and supply may alter: 1) the dissolved oxygen concentrations in the benthic foraminifera microhabitats, in turn impacting foraminiferal distributions, or 2) OM resource allocation and exploitation by different foraminifera genera in KSEs (Cresswell and van Hengstum, 2021). While in-situ microcosm experiments could help resolve the precise mechanisms, these findings indicate that changes in OM delivery into Aerolito over time have impacted benthic meiofaunal communities.

Conclusion

• Only young sediment from the last 1,000 years is preserved in Aerolito, despite the fact that the cave was flooded by the coastal aquifer thousands of years prior from deglacial sea-level rise. More frequent intense hurricane activity prior to 1,000 years ago may have exported older sediment out of the cave from constant increases in groundwater flow and discharge rates.

• The source and quantity of OM exported from the subterranean estuary to adjacent ocean is not constant through time, and its export is linked to: 1) changes in adjacent aquatic and terrestrial habitats and 2) potentially large storm events. In contrast, the coastal caves have allowed marine-derived OM to enter the aquifer in spite of the net land-to-sea hydrologic transport from the discharge of coastal groundwater.

• Over the last millennium, a decrease in terrestrial OM deposition (reduction in C:N ratios and increase in δ¹³C_org values) in the cave suggests that the magnitude of OM inputs from adjacent subaerial environments has decreased. We hypothesize that the contribution of ¹³C-depleted OM from mangroves in the adjacent sinkhole has changed through time, caused by a decrease in the areal extent of the mangroves in the area.

• Benthic foraminiferal community changes reveal an ecological response to decreased OM supply to the cave through time. In the past, when OM supply was higher and dominated by non-marine sources (i.e., depleted δ¹³C_org values), the benthic foraminiferal community in Aerolito was composed of taxa preferentially found in sediment with ¹³C-depleted bulk OM in Bermudian caves that have well-oxygenated water columns. The hypothesized cause of these community shifts is a change in food resource partitioning between foraminiferal genera through time.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

The project was conceived by PvH and LMM-O. PvH, TW, and GY-M collected the samples. DB, SL, TW, AT, and CM completed the laboratory activities. DB and SL conducted data analyses. DB, TW, and SL generated the artwork. All authors contributed to data interpretation, and DB wrote the manuscript with contributions from all.

Funding

Funding for DB was provided by the Texas A&M University at Galveston Postdoctoral Fellowship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

John W. Pohlman and Fernando Calderón-Gutiérrez provided valuable feedback during the writing. Two reviewers are thanked for comments to improve the final manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2021.670914/full#supplementary-material

References

Airoldi, L., and Cinelli, F. (1996). Variability of Fluxes of Particulate Material in a Submarine Cave with Chemolithoautotrophic Inputs of Organic Carbon. Mar. Ecol. Prog. Ser. 139, 205–217. doi:10.3354/meps139205

CrossRef Full Text | Google Scholar

Alve, E., and Bernhard, J. M. (1995). Vertical Migratory Response of Benthic Foraminifera to Controlled Oxygen Concentrations in an Experimental Mesocosm. Mar. Ecol. Prog. Ser. 137–151.

Google Scholar

Bauer-Gottwein, P., Gondwe, B. R. N., Charvet, G., Marín, L. E., Rebolledo-Vieyra, M., and Merediz-Alonso, G. (2011). Review: The Yucatán Peninsula Karst Aquifer, Mexico. Hydrogeol J. 19, 507–524. doi:10.1007/s10040-010-0699-5

CrossRef Full Text | Google Scholar

Beck, A. J., Charette, M. A., Cochran, J. K., Gonneea, M. E., and Peucker-Ehrenbrink, B. (2013). Dissolved Strontium in the Subterranean Estuary - Implications for the Marine Strontium Isotope Budget. Geochimica et Cosmochimica Acta 117, 33–52. doi:10.1016/j.gca.2013.03.021

CrossRef Full Text | Google Scholar

Beddows, P. A., Smart, P. L., Whitaker, F. F., and Smith, S. L. (2007). Decoupled Fresh-Saline Groundwater Circulation of a Coastal Carbonate Aquifer: Spatial Patterns of Temperature and Specific Electrical Conductivity. J. Hydrol. 346, 18–32. doi:10.1016/j.jhydrol.2007.08.013

CrossRef Full Text | Google Scholar

Belanger, C. L., Sharon, J., Du, J., Payne, C. R., and Mix, A. C. (2020). North Pacific Deep-Sea Ecosystem Responses Reflect Post-glacial Switch to Pulsed Export Productivity, Deoxygenation, and Destratification. Deep Sea Res. Oceanographic Res. Pap. 164, 103341. doi:10.1016/j.dsr.2020.103341

CrossRef Full Text | Google Scholar

Benner, R., Fogel, M. L., and Sprague, E. K. (1991). Diagenesis of Belowground Biomass of Spartina Alterniflora in Salt-Marsh Sediments. Limnol. Oceanogr. 36, 1358–1374. doi:10.4319/lo.1991.36.7.1358

CrossRef Full Text | Google Scholar

Bentley, S. J., Keen, T. R., Blain, C. A., and Vaughan, W. C. (2002). The Origin and Preservation of a Major Hurricane Event Bed in the Northern Gulf of Mexico: Hurricane Camille, 1969. Mar. Geol. 186, 423–446. doi:10.1016/s0025-3227(02)00297-9

CrossRef Full Text | Google Scholar

Bergamin, L., Taddei Ruggiero, E., Pierfranceschi, G., Andres, B., Constantino, R., Crovato, C., et al. (2020). Benthic Foraminifera and Brachiopods from a Marine Cave in Spain: Environmental Significance. Medit. Mar. Sci. 21, 506–518. doi:10.12681/mms.23482

CrossRef Full Text | Google Scholar

Bernhard, J. M. (1992). Benthic Foraminiferal Distribution and Biomass Related to Pore-Water Oxygen Content: Central California Continental Slope and Rise. Deep Sea Research Part A. Oceanographic Res. Pap. 39 (3-4), 585–605. doi:10.1016/0198-0149(92)90090-g

CrossRef Full Text | Google Scholar

Bernhard, J. M., and Bowser, S. S. (1999). Benthic Foraminifera of Dysoxic Sediments: Chloroplast Sequestration and Functional Morphology. Earth Sci. Rev. 46 (1-4), 149–165. doi:10.1016/s0012-8252(99)00017-3

CrossRef Full Text | Google Scholar

Bishop, R. E., Humphreys, W. F., Kršinić, F., Sket, B., Iliffe, T. M., Žic, V., et al. (2015). 'Anchialine' Redefined as a Subterranean Estuary in a Crevicular or Cavernous Geological Setting. J. Crustacean Biol. 35, 511–514. doi:10.1163/1937240x-00002335

CrossRef Full Text | Google Scholar

Björnerås, C., Škerlep, M., Gollnisch, R., Herzog, S. D., Ugge, G. E., Hegg, A., et al. (2020). Inland Blue Holes of the Bahamas–chemistry and Biology in a Unique Aquatic Environment. Fundam. Appl. Limnol. 95–106. doi:10.1127/fal/2020/1330

CrossRef Full Text | Google Scholar

Brankovits, D., Pohlman, J. W., Ganju, N. K., Iliffe, T. M., Lowell, N., Roth, E., et al. (2018). Hydrologic Controls of Methane Dynamics in Karst Subterranean Estuaries. Glob. Biogeochem. Cycles 32, 1759–1775. doi:10.1029/2018gb006026

CrossRef Full Text | Google Scholar

Brankovits, D., Pohlman, J. W., Niemann, H., Leigh, M. B., Leewis, M.-C., Becker, K. W., et al. (2017). Methane-and Dissolved Organic Carbon-Fueled Microbial Loop Supports a Tropical Subterranean Estuary Ecosystem. Nat. Commun. 8, 1835. doi:10.1038/s41467-017-01776-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Calderón-Gutiérrez, F., Sánchez-Ortiz, C., and Huato-Soberanis, L. (2018). Ecological Patterns in Anchialine Caves. PLoS ONE 13. doi:10.1371/journal.pone.0202909

CrossRef Full Text | Google Scholar

Caralp, M. H. (1989). Size and Morphology of the Benthic Foraminifer Melonis Barleeanum; Relationships with Marine Organic Matter. J. Foraminiferal Res. 19 (3), 235–245. doi:10.2113/gsjfr.19.3.235

CrossRef Full Text | Google Scholar

Charette, M. A., and Sholkovitz, E. R. (2002). Oxidative Precipitation of Groundwater-Derived Ferrous Iron in the Subterranean Estuary of a Coastal Bay. Geophys. Res. Lett. 29, 8585–9184. doi:10.1029/2001gl014512

CrossRef Full Text | Google Scholar

Chen, X., Cukrov, N., Santos, I. R., Rodellas, V., Cukrov, N., and Du, J. (2020). Karstic Submarine Groundwater Discharge into the Mediterranean: Radon-Based Nutrient Fluxes in an Anchialine Cave and a Basin-wide Upscaling. Geochimica et Cosmochimica Acta 268, 467–484. doi:10.1016/j.gca.2019.08.019

CrossRef Full Text | Google Scholar

Collins, S. V., Reinhardt, E. G., Werner, C. L., Le Maillot, C., Devos, F., and Rissolo, D. (2015). Late Holocene Mangrove Development and Onset of Sedimentation in the Yax Chen Cave System (Ox Bel Ha) Yucatan, Mexico: Implications for Using Cave Sediments as a Sea-Level Indicator. Palaeogeogr. Palaeoclimatol. Palaeoecol. 438, 124–134. doi:10.1016/j.palaeo.2015.07.042

CrossRef Full Text | Google Scholar

Corliss, B. H., and Chen, C. (1988). Morphotype Patterns of Norwegian Sea Deep-Sea Benthic Foraminifera and Ecological Implications. Geol 16 (8), 716–719. doi:10.1130/0091-7613(1988)016<0716:mponsd>2.3.co;2

CrossRef Full Text | Google Scholar

Corliss, B. H., and Emerson, S. (1990). Distribution of Rose Bengal Stained Deep-Sea Benthic Foraminifera from the Nova Scotian Continental Margin and Gulf of Maine. Deep Sea Res. A. Oceanographic Res. Pap. 37 (3), 381–400. doi:10.1016/0198-0149(90)90015-n

CrossRef Full Text | Google Scholar

Corliss, B. H. (1991). Morphology and Microhabitat Preferences of Benthic Foraminifera from the Northwest Atlantic Ocean. Mar. Micropaleontol. 17 (3-4), 195–236. doi:10.1016/0377-8398(91)90014-w

CrossRef Full Text | Google Scholar

Coronado-Álvarez, L., Gutiérrez-Aguirre, M. A., and Cervantes-Martínez, A. (2011). Water Quality in Wells from Cozumel Island, Mexico. Trop. Subtropical Agroecosystems 13, 233–241.

Google Scholar

Coutino, A., Stastna, M., Kovacs, S., and Reinhardt, E. (2017). Hurricanes Ingrid and Manuel (2013) and Their Impact on the Salinity of the Meteoric Water Mass, Quintana Roo, Mexico. J. Hydrol. 551, 715–729. doi:10.1016/j.jhydrol.2017.04.022

CrossRef Full Text | Google Scholar

Coutino, A., Stastna, M., and Reinhardt, E. G. (2020). Interaction of Mangrove Surface Coverage and Groundwater Inputs on the Temperature and Water Level Near Tulum, Quintana Roo, Mexico: Observations and Modelling. J. Hydrol. 583, 124566. doi:10.1016/j.jhydrol.2020.124566

CrossRef Full Text | Google Scholar

Cresswell, J. N., and van Hengstum, P. J. (2021). Habitat Partitioning in the Marine Sector of Karst Subterranean Estuaries and Bermuda's Marine Caves: Benthic Foraminiferal Evidence. Front. Environ. Sci. 8, 249. doi:10.3389/fenvs.2020.594554

CrossRef Full Text | Google Scholar

Cresswell, J. N., van Hengstum, P. J., Iliffe, T. M., Williams, B. E., and Nolan, G. (2017). Anthropogenic Infilling of a Bermudian Sinkhole and its Impact on Sedimentation and Benthic Foraminifera in the Adjacent Anchialine Cave Environment. Int. J. Speleol. 46 (3), 7. doi:10.5038/1827-806x.46.3.2128

CrossRef Full Text | Google Scholar

Crook, E. D., Potts, D., Rebolledo-Vieyra, M., Hernandez, L., and Paytan, A. (2012). Calcifying Coral Abundance Near Low-pH Springs: Implications for Future Ocean Acidification. Coral Reefs 31, 239–245. doi:10.1007/s00338-011-0839-y

CrossRef Full Text | Google Scholar

Curtis, J. H., Hodell, D. A., and Brenner, M. (1996). Climate Variability on the Yucatan Peninsula (Mexico) during the Past 3500 years, and Implications for Maya Cultural Evolution. Quat. Res. 46, 37–47. doi:10.1006/qres.1996.0042

CrossRef Full Text | Google Scholar

Cushman, J. A. (1933). Some New Recent Foraminifera from the Tropical Pacific. Contrib. Cushman Lab. Foraminiferal Res. 9, 77–95.

Google Scholar

Dean, W. E. (1974). Determination of Carbonate and Organic Matter in Calcareous Sediments and Sedimentary Rocks by Loss on Ignition; Comparison with Other Methods. J. Sediment. Res. 44, 242–248. doi:10.1306/74D729D2-2B21-11D7-8648000102C1865D

Google Scholar

Dimiza, M. D., Fatourou, M., Arabas, A., Panagiotopoulos, I., Gogou, A., Kouli, K., et al. (2020). Deep-sea Benthic Foraminifera Record of the Last 1500 Years in the North Aegean Trough (Northeastern Mediterranean): A Paleoclimatic Reconstruction Scenario. Deep Sea Res. Part Topical Stud. Oceanography 171, 104705. doi:10.1016/j.dsr2.2019.104705

CrossRef Full Text | Google Scholar

Dittmar, T., Hertkorn, N., Kattner, G., and Lara, R. J. (2006). Mangroves, a Major Source of Dissolved Organic Carbon to the Oceans. Glob. Biogeochem. Cycles 20, a–n. doi:10.1029/2005gb002570

CrossRef Full Text | Google Scholar

Donato, D. C., Kauffman, J. B., Murdiyarso, D., Kurnianto, S., Stidham, M., and Kanninen, M. (2011). Mangroves Among the Most Carbon-Rich Forests in the Tropics. Nat. Geosci 4, 293–297. doi:10.1038/ngeo1123

CrossRef Full Text | Google Scholar

Donnelly, J. P., and Woodruff, J. D. (2007). Intense Hurricane Activity over the Past 5,000 Years Controlled by El Niño and the West African Monsoon. Nature 447, 465–468. doi:10.1038/nature05834

PubMed Abstract | CrossRef Full Text | Google Scholar

Fatela, F., and Taborda, R. (2002). Confidence Limits of Species Proportions in Microfossil Assemblages. Mar. Micropaleontol. 45 (2), 169–174. doi:10.1016/s0377-8398(02)00021-x

CrossRef Full Text | Google Scholar

Fleury, P., Bakalowicz, M., and de Marsily, G. (2007). Submarine Springs and Coastal Karst Aquifers: a Review. J. Hydrol. 339, 79–92. doi:10.1016/j.jhydrol.2007.03.009

CrossRef Full Text | Google Scholar

Fogel, M. L., Kent Sprague, E., Gize, A. P., and Frey, R. W. (1989). Diagenesis of Organic Matter in georgia Salt Marshes. Estuarine, Coastal Shelf Sci. 28, 211–230. doi:10.1016/0272-7714(89)90067-x

CrossRef Full Text | Google Scholar

Fontanier, C., Jorissen, F. J., Licari, L., Alexandre, A., Anschutz, P., and Carbonel, P. (2002). Live Benthic Foraminiferal Faunas from the Bay of Biscay: Faunal Density, Composition, and Microhabitats. Deep Sea Res. Part Oceanographic Res. Pap. 49 (4), 751–785. doi:10.1016/s0967-0637(01)00078-4

CrossRef Full Text | Google Scholar

Ford, D., and Williams, P. D. (2013). Karst Hydrogeology and Geomorphology. West Sussex: John Wiley & Sons.

Fosshagen, A., and Iliffe, T. M. (1991). A New Genus of Calanoid Copepod from an Anchialine Cave in Belize. Bull. Plankton Soc. Jpn. 27. 339-350.

Google Scholar

Fourqurean, J. W., Escorcia, S. P., Anderson, W. T., and Zieman, J. C. (2005). Spatial and Seasonal Variability in Elemental Content, δ13C, and δ15N ofThalassia Testudinum from South Florida and its Implications for Ecosystem Studies. Estuaries 28, 447–461. doi:10.1007/bf02693926

CrossRef Full Text | Google Scholar

Frausto-Martínez, O., Castillo, J. F. R., and Olivares, O. C. (2021). Morphometry of Karst Depressions at Detailed Scale: El Cedral, Cozumel, Mexico. Trop. Subtropical Agroecosystems 24.

Google Scholar

Frausto-Martínez, O., Zapi-Salazar, N. A., and Colin-Olivares, O. (2018). Identification of Karst Forms Using LiDAR Technology: Cozumel Island, Mexico. Trends Geomatics-an Earth Sci. Perspective 11. doi:10.5772/intechopen.79196

CrossRef Full Text | Google Scholar

Frontana-Uribe, S., and Soli's-Weiss, V. (2011). First Records of Polychaetous Annelids from Cenote Aerolito (Sinkhole and Anchialine Cave) in Cozumel Island. Jcks 73, 1–10. doi:10.4311/jcks2009lsc0107

CrossRef Full Text | Google Scholar

Fry, B. (2007). Stable Isotope Ecology. New York: Springer Science & Business Media.

Gabriel, J. J., Reinhardt, E. G., Peros, M. C., Davidson, D. E., van Hengstum, P. J., and Beddows, P. A. (2009). Palaeoenvironmental Evolution of Cenote Aktun Ha (Carwash) on the Yucatan Peninsula, Mexico and its Response to Holocene Sea-Level Rise. J. Paleolimnol 42, 199–213. doi:10.1007/s10933-008-9271-x

CrossRef Full Text | Google Scholar

Gómez, P., and Calderón-gutiérrez, F. (2020). Anchialine Cave-Dwelling Sponge Fauna (Porifera) from La Quebrada, Mexico, with the Description of the First Mexican Stygobiont Sponges. Zootaxa 4803, 125–151. doi:10.11646/zootaxa.4803.1.7

CrossRef Full Text | Google Scholar

Gondwe, B. R. N., Lerer, S., Stisen, S., Marín, L., Rebolledo-Vieyra, M., Merediz-Alonso, G., et al. (2010). Hydrogeology of the South-Eastern Yucatan Peninsula: New Insights from Water Level Measurements, Geochemistry, Geophysics and Remote Sensing. J. Hydrol. 389, 1–17. doi:10.1016/j.jhydrol.2010.04.044

CrossRef Full Text | Google Scholar

Goñi, M. A., Ruttenberg, K. C., and Eglinton, T. I. (1997). Sources and Contribution of Terrigenous Organic Carbon to Surface Sediments in the Gulf of Mexico. Nature 389, 275–278. doi:10.1038/38477

CrossRef Full Text | Google Scholar

Gonneea, M. E., Charette, M. A., Liu, Q., Herrera-Silveira, J. A., and Morales-Ojeda, S. M. (2014). Trace Element Geochemistry of Groundwater in a Karst Subterranean Estuary (Yucatan Peninsula, Mexico). Geochimica et Cosmochimica Acta 132, 31–49. doi:10.1016/j.gca.2014.01.037

CrossRef Full Text | Google Scholar

Gonneea, M. E., Maio, C. V., Kroeger, K. D., Hawkes, A. D., Mora, J., Sullivan, R., et al. (2019). Salt Marsh Ecosystem Restructuring Enhances Elevation Resilience and Carbon Storage during Accelerating Relative Sea-Level Rise. Estuarine, Coastal Shelf Sci. 217, 56–68. doi:10.1016/j.ecss.2018.11.003

CrossRef Full Text | Google Scholar

Gonneea, M. E., Paytan, A., and Herrera-Silveira, J. A. (2004). Tracing Organic Matter Sources and Carbon Burial in Mangrove Sediments over the Past 160 Years. Estuarine, Coastal Shelf Sci. 61, 211–227. doi:10.1016/j.ecss.2004.04.015

CrossRef Full Text | Google Scholar

Gooday, A. J. (2003). Benthic Foraminifera (Protista) as Tools in Deep-Water Palaeoceanography: Environmental Influences on Faunal Characteristics. Adv. Mar. Biol. 46, 1–90. doi:10.1016/s0065-2881(03)46002-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gooday, A. J. (1996). Epifaunal and Shallow Infaunal Foraminiferal Communities at Three Abyssal NE Atlantic Sites Subject to Differing Phytodetritus Input Regimes. Deep Sea Res. Part Oceanographic Res. Pap. 43 (9), 1395–1421. doi:10.1016/s0967-0637(96)00072-6

CrossRef Full Text | Google Scholar

Gooday, A. J., and Turley, C. M. (1990). Responses by Benthic Organisms to Inputs of Organic Material to the Ocean Floor: a Review. Philosophical Trans. R. Soc. A 331 (1616), 119–138. doi:10.1098/rsta.1990.0060

CrossRef Full Text | Google Scholar

Gregory, B. R. B., Reinhardt, E. G., and Gifford, J. A. (2017). The Influence of Morphology on Sinkhole Sedimentation at Little Salt Spring, Florida. J. Coastal Res. 332, 359–371. doi:10.2112/jcoastres-d-15-00169.1

CrossRef Full Text | Google Scholar

Heiri, O., Lotter, A. F., and Lemcke, G. (2001). Loss on Ignition as a Method for Estimating Organic and Carbonate Content in Sediments: Reproducibility and Comparability of Results. J. Paleolimnol. 25, 101–110. doi:10.1023/a:1008119611481

CrossRef Full Text | Google Scholar

Holthuis, L. B. (1973). Caridean Shrimps found in Land-Locked Saltwater Pools at four Indo-West Pacific Localities (Sinai Peninsula, Funafuti Atoll, Maui and Hawaii Islands), with the description of one new genus and four new species. Zoologische Verhandelingen 128 (1), 1–48.

Iliffe, T. M., and Kornicker, L. S. (2009). Worldwide Diving Discoveries of Living Fossil Animals from the Depths of Anchialine and Marine Caves. Smithsonian Contrib. Mar. Sci. 38, 269–280.

Google Scholar

Jorissen, F. J., de Stigter, H. C., and Widmark, J. G. (1995). A Conceptual Model Explaining Benthic Foraminiferal Microhabitats. Mar. Micropaleontol. 26 (1-4), 3–15. doi:10.1016/0377-8398(95)00047-x

CrossRef Full Text | Google Scholar

Jorissen, F. J., Fontanier, C., and Thomas, E. (2007). Chapter Seven Paleoceanographical Proxies Based on Deep-Sea Benthic Foraminiferal Assemblage Characteristics. Dev. Mar. Geol. 1, 263–325. doi:10.1016/s1572-5480(07)01012-3

CrossRef Full Text | Google Scholar

Jorissen, F. J., and Wittling, I. (1999). Ecological Evidence from Live–Dead Comparisons of Benthic Foraminiferal Faunas off Cape Blanc (Northwest Africa). Palaeogeogr. Palaeoclimatol. Palaeoecol. 149 (1-4), 151–170. doi:10.1016/s0031-0182(98)00198-9

CrossRef Full Text | Google Scholar

Kaiho, K. (1994). Benthic Foraminiferal Dissolved-Oxygen Index and Dissolved-Oxygen Levels in the Modern Ocean. Geol 22 (8), 719–722. doi:10.1130/0091-7613(1994)022<0719:bfdoia>2.3.co;2

CrossRef Full Text | Google Scholar

Keil, R. G., Tsamakis, E., Fuh, C. B., Giddings, J. C., and Hedges, J. I. (1994). Mineralogical and Textural Controls on the Organic Composition of Coastal Marine Sediments: Hydrodynamic Separation Using SPLITT-Fractionation. Geochimica et Cosmochimica Acta 58, 879–893. doi:10.1016/0016-7037(94)90512-6

CrossRef Full Text | Google Scholar

Khan, N. S., Vane, C. H., Engelhart, S. E., Kendrick, C., and Horton, B. P. (2019). The Application of δ13C, TOC and C/N Geochemistry of Mangrove Sediments to Reconstruct Holocene Paleoenvironments and Relative Sea Levels, Puerto Rico. Mar. Geol. 415, 105963. doi:10.1016/j.margeo.2019.105963

CrossRef Full Text | Google Scholar

Kjerfve, B. (1981). Tides of the Caribbean Sea. J. Geophys. Res. 86, 4243–4247. doi:10.1029/jc086ic05p04243

CrossRef Full Text | Google Scholar

Kottek, M., Grieser, J., Beck, C., Rudolf, B., and Rubel, F. (2006). World Map of the Köppen-Geiger Climate Classification Updated. metz 15, 259–263. doi:10.1127/0941-2948/2006/0130

CrossRef Full Text | Google Scholar

Kovacs, S. E., Reinhardt, E. G., Chatters, J. C., Rissolo, D., Schwarcz, H. P., Collins, S. V., et al. (2017). Calcite Raft Geochemistry as a Hydrological Proxy for Holocene Aquifer Conditions in Hoyo Negro and Ich Balam (Sac Actun Cave System), Quintana Roo, Mexico. Quat. Sci. Rev. 175, 97–111. doi:10.1016/j.quascirev.2017.09.006

CrossRef Full Text | Google Scholar

Kuenen, P. H., and Migliorini, C. I. (1950). Turbidity Currents as a Cause of Graded Bedding. J. Geol. 58, 91–127. doi:10.1086/625710

CrossRef Full Text | Google Scholar

Lamb, A. L., Wilson, G. P., and Leng, M. J. (2006). A Review of Coastal Palaeoclimate and Relative Sea-Level Reconstructions Using δ13C and C/N Ratios in Organic Material. Earth-Science Rev. 75, 29–57. doi:10.1016/j.earscirev.2005.10.003

CrossRef Full Text | Google Scholar

Little, S. N., and van Hengstum, P. J. (2019). Intertidal and Subtidal Benthic Foraminifera in Flooded Caves: Implications for Reconstructing Coastal Karst Aquifers and Cave Paleoenvironments. Mar. Micropaleontol. 149, 19–34. doi:10.1016/j.marmicro.2019.03.005

CrossRef Full Text | Google Scholar

Loeblich, A. R., and Tappan, H. (1987). Foraminiferal Genera and Their Classification. New York: Springer.

Luijendijk, E., Gleeson, T., and Moosdorf, N. (2020). Fresh Groundwater Discharge Insignificant for the World’s Oceans but Important for Coastal Ecosystems. Nat. Commun. 11, 1–12. doi:10.1038/s41467-020-15064-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Maher, D. T., Santos, I. R., Golsby-Smith, L., Gleeson, J., and Eyre, B. D. (2013). Groundwater-derived Dissolved Inorganic and Organic Carbon Exports from a Mangrove Tidal Creek: The Missing Mangrove Carbon Sink?. Limnol. Oceanogr. 58, 475–488. doi:10.4319/lo.2013.58.2.0475

CrossRef Full Text | Google Scholar

Márquez-Borrás, F., Solís-Marín, F. A., Bribiesca-Contreras, G., and Laguarda-Figueras, A. (2016). First Record of Ophiura Ljungmani (Echinodermata: Ophiuroidea) from an Anchialine Cave in the Mexican Caribbean. Revista Mexicana de Biodiversidad 87, 1127–1130. doi:10.1016/j.rmb.2016.07.006

CrossRef Full Text | Google Scholar

Mayfield, K. K., Eisenhauer, A., Ramos, D. P. S., Higgins, J. A., Horner, T. J., Auro, M., et al. (2021). Groundwater Discharge Impacts Marine Isotope Budgets of Li, Mg, Ca, Sr, and Ba. Nat. Commun. 12, 1–9. doi:10.1038/s41467-020-20248-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mejía, L., Zarza, E., and López, M. (2008). Barbouria Yanezi Sp. nov., a New Species of Cave Shrimp (Decapoda, Barbouriidae) from Cozumel Island, Mexico. Crustaceana 81, 663–672. doi:10.1163/156854008784513474

CrossRef Full Text | Google Scholar

Mejía-ortíz, L. M., Yañez, G., and López-Mejía, M. (2017). Anchialocarididae, a New Family of Anchialine Decapods and a New Species of the Genus Agostocaris from Cozumel Island, Mexico. Crustac 90, 381–398. doi:10.1163/15685403-00003657

CrossRef Full Text | Google Scholar

Mejía-Ortíz, L. M., Yánez, G., and López-Mejía, M. (2007). Echinoderms in an Anchialine Cave in Mexico. Mar. Ecol. 28, 31–34. doi:10.1111/j.1439-0485.2007.00174.x

CrossRef Full Text | Google Scholar

Milne, G. A., and Peros, M. (2013). Data-model Comparison of Holocene Sea-Level Change in the Circum-Caribbean Region. Glob. Planet. Change 107, 119–131. doi:10.1016/j.gloplacha.2013.04.014

CrossRef Full Text | Google Scholar

Moodley, L., and Hess, C. (1992). Tolerance of Infaunal Benthic Foraminifera for Low and High Oxygen Concentrations. Biol. Bull. 183 (1), 94–98. doi:10.2307/1542410

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, W. S., and Joye, S. B. (2021). Saltwater Intrusion and Submarine Groundwater Discharge: Acceleration of Biogeochemical Reactions in Changing Coastal Aquifers. Front. Earth Sci. 9. doi:10.3389/feart.2021.600710

CrossRef Full Text | Google Scholar

Moore, W. S. (2010). The Effect of Submarine Groundwater Discharge on the Ocean. Annu. Rev. Mar. Sci. 2, 59–88. doi:10.1146/annurev-marine-120308-081019

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, W. S. (1999). The Subterranean Estuary: a Reaction Zone of Ground Water and Sea Water. Mar. Chem. 65, 111–125. doi:10.1016/s0304-4203(99)00014-6

CrossRef Full Text | Google Scholar

Moosdorf, N., Stieglitz, T., Waska, H., Dürr, H. H., and Hartmann, J. (2015). Submarine Groundwater Discharge from Tropical Islands: a Review. Grundwasser 20, 53–67. doi:10.1007/s00767-014-0275-3

CrossRef Full Text | Google Scholar

Muckelbauer, G. (1990). The Shelf of Cozumel, Mexico: Topography and Organisms. Facies 23, 201–239. doi:10.1007/bf02536714

CrossRef Full Text | Google Scholar

Muhs, D. R., Budahn, J. R., Prospero, J. M., and Carey, S. N. (2007). Geochemical Evidence for African Dust Inputs to Soils of Western Atlantic Islands: Barbados, the Bahamas, and Florida. J. Geophys. Res. Earth Surf. 112. doi:10.1029/2005jf000445

CrossRef Full Text | Google Scholar

Murray, J. W. (2003). An Illustrated Guide to the Benthic Foraminifera of the Hebridean Shelf, West of Scotland, with Notes on Their Mode of Life. Palaeontol. Electronica 5 (1), 31.

Google Scholar

Murray, M. R. (2002). Is Laser Particle Size Determination Possible for Carbonate-Rich Lake Sediments?. J. Paleolimnol. 27 (2), 173–183. doi:10.1023/a:1014281412035

CrossRef Full Text | Google Scholar

Null, K. A., Knee, K. L., Crook, E. D., de Sieyes, N. R., Rebolledo-Vieyra, M., Hernández-Terrones, L., et al. (2014). Composition and Fluxes of Submarine Groundwater along the Caribbean Coast of the Yucatan Peninsula. Continental Shelf Res. 77, 38–50. doi:10.1016/j.csr.2014.01.011

CrossRef Full Text | Google Scholar

Oksanen, J., Guillaume Blanchet, F., Kindt, R., Legendre, P., McGlinn, D., Minchin, P. R., et al. (2019). Vegan: Community Ecology Package. Retrieved from https://CRAN.R-project.org/package=vegan (Accessed January 24, 2021).

Google Scholar

Olesen, J., Meland, K., Glenner, H., van Hengstum, P. J., and Iliffe, T. M. (2017). Xibalbanus Cozumelensis, a New Species of Remipedia (Crustacea) from Cozumel, Mexico, and a Molecular Phylogeny of Xibalbanus on the Yucatán Peninsula. Eur. J. Taxonomy 316, 1–27. doi:10.5852/ejt.2017.316

CrossRef Full Text | Google Scholar

Ortiz, M., and Winfield, I. (2015). A New Amphipod Species (Peracarida: Amphipoda: Ampithoidae) Collected from Cenote Aerolito, Cozumel Island, Quintana Roo. Revista Mexicana de Biodiversidad 86, 332–336. doi:10.1016/j.rmb.2015.04.008

CrossRef Full Text | Google Scholar

Patterson, R. T., and Fishbein, E. (1989). Re-examination of the Statistical Methods Used to Determine the Number of Point Counts Needed for Micropaleontological Quantitative Research. J. Paleontol. 63 (02), 245–248. doi:10.1017/s0022336000019272

CrossRef Full Text | Google Scholar

Perry, E., Velazquez-Oliman, G., and Marin, L. (2002). The Hydrogeochemistry of the Karst Aquifer System of the Northern Yucatan Peninsula, Mexico. Int. Geol. Rev. 44, 191–221. doi:10.2747/0020-6814.44.3.191

CrossRef Full Text | Google Scholar

Pisanty, D. T., Solis, E. M. C., Weiss, V. S., and Salazar, M. H. (2010). Peracarids (Crustacea: Malacostraca) from Cenote Aerolito, Cozumel, Mexican Caribbean. Cahiers de Biologie Mar. 51, 177–180.

Google Scholar

Pohlman, J., Iliffe, T., and Cifuentes, L. (1997). A Stable Isotope Study of Organic Cycling and the Ecology of an Anchialine Cave Ecosystem. Mar. Ecol. Prog. Ser. 155, 17–27. doi:10.3354/meps155017

CrossRef Full Text | Google Scholar

Pohlman, J. W. (2011). The Biogeochemistry of Anchialine Caves: Progress and Possibilities. Hydrobiologia 677, 33–51. doi:10.1007/s10750-011-0624-5

CrossRef Full Text | Google Scholar

Prospero, J. M., and Lamb, P. J. (2003). African Droughts and Dust Transport to the Caribbean: Climate Change Implications. Science 302, 1024–1027. doi:10.1126/science.1089915

PubMed Abstract | CrossRef Full Text | Google Scholar

Prospero, J. M., and Nees, R. T. (1986). Impact of the North African Drought and El Niño on Mineral Dust in the Barbados Trade Winds. Nature 320, 735–738. doi:10.1038/320735a0

CrossRef Full Text | Google Scholar

R Core Team (2020). RStudio: Integrated Development for R. USA: RStudio, Inc. Boston. 423.

Rathburn, A. E., and Corliss, B. H. (1994). The Ecology of Living (Stained) Deep-Sea Benthic Foraminifera from the Sulu Sea. Paleoceanography 9 (1), 87–150. doi:10.1029/93pa02327

CrossRef Full Text | Google Scholar

Reimer, P. J., Austin, W. E. N., Bard, E., Bayliss, A., Blackwell, P. G., Bronk Ramsey, C., et al. (2020). The IntCal20 Northern Hemisphere Radiocarbon Age Calibration Curve (0-55 Cal kBP). Radiocarbon 62, 725–757. doi:10.1017/rdc.2020.41

CrossRef Full Text | Google Scholar

Richards, D. A., Smart, P. L., and Lawrence Edwards, R. (1994). Maximum Sea Levels for the Last Glacial Period from U-Series Ages of Submerged Speleothems. Nature 367, 357–360. doi:10.1038/367357a0

CrossRef Full Text | Google Scholar

Riera, R., Monterroso, Ó., Núñez, J., and Martínez, A. (2018). Distribution of Meiofaunal Abundances in a Marine Cave Complex with Secondary Openings and Freshwater Filtrations. Mar. Biodiv 48, 203–215. doi:10.1007/s12526-016-0586-y

CrossRef Full Text | Google Scholar

Rocha, C., Robinson, C. E., Santos, I. R., Waska, H., Michael, H. A., and Bokuniewicz, H. J. (2021). A Place for Subterranean Estuaries in the Coastal Zone. Estuarine, Coastal Shelf Sci. 250, 107167. doi:10.1016/j.ecss.2021.107167

CrossRef Full Text | Google Scholar

Rubio, F., Rolán, E., Worsaae, K., Martínez, A., and Gonzalez, B. C. (2015). Description of the First Anchialine Gastropod from a Yucatán Cenote, Teinostoma Brankovitsi N. sp. (Caenogastropoda: Tornidae), Including an Emended Generic Diagnosis. J. Molluscan Stud. 82, 169–177. doi:10.1093/mollus/eyv049

CrossRef Full Text | Google Scholar

Schönfeld, J., Alve, E., Geslin, E., Jorissen, F., Korsun, S., and Spezzaferri, S. (2012). The FOBIMO (FOraminiferal BIo-MOnitoring) Initiative-Towards a Standardised Protocol for Soft-Bottom Benthic Foraminiferal Monitoring Studies. Mar. Micropaleontol. 94-95, 1–13. doi:10.1016/j.marmicro.2012.06.001

CrossRef Full Text | Google Scholar

Scott, D. B., and Hermelin, J. O. R. (1993). A Device for Precision Splitting of Micropaleontological Samples in Liquid Suspension. J. Paleontol. 67 (01), 151–154. doi:10.1017/s0022336000021302

CrossRef Full Text | Google Scholar

Sen Gupta, B. K., Eugene Turner, R., and Rabalais, N. N. (1996). Seasonal Oxygen Depletion in Continental-Shelf Waters of Louisiana: Historical Record of Benthic Foraminifers. Geol 24 (3), 227–230. doi:10.1130/0091-7613(1996)024<0227:sodics>2.3.co;2

CrossRef Full Text | Google Scholar

Seymour, J., Humphreys, W., and Mitchell, J. (2007). Stratification of the Microbial Community Inhabiting an Anchialine Sinkhole. Aquat. Microb. Ecol. 50, 11–24. doi:10.3354/ame01153

CrossRef Full Text | Google Scholar

Sket, B. (1996). The Ecology of Anchihaline Caves. Trends Ecol. Evol. 11, 221–225. doi:10.1016/0169-5347(96)20031-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, M. R. (2017). Ternary: An R Package for Creating Ternary Plots. Zenodo 10. doi:10.5281/zenodo.1068996

CrossRef Full Text | Google Scholar

Socki, R. A., Perry, E. J. C., and Romanek, C. S. (2002). Stable Isotope Systematics of Two Cenotes from the Northern Yucatan Peninsula, Mexico. Limnol. Oceanogr. 47, 1808–1818. doi:10.4319/lo.2002.47.6.1808

CrossRef Full Text | Google Scholar

Spaw, R. H. (1978). Late Pleistocene Carbonate Bank Deposition: Cozumel Island, Quintana Roo, Mexico. Gulf Coast Assoc. Geol. Societies Trans. 28, 601–619.

Google Scholar

Stock, J. H., Iliffe, T. M., and Williams, D. (1986). The Concept of Anchialine Reconsidered. Stygologia 2, 90–92.

Google Scholar

Stuiver, M., Reimer, P. J., and Reimer, R. W. (2021). CALIB 8.2. Retrieved from http://calib.org (Accessed January 24, 2021).

Google Scholar

Tamalavage, A. E., van Hengstum, P. J., Louchouarn, P., Molodtsov, S., Kaiser, K., Donnelly, J. P., et al. (2018). Organic Matter Sources and Lateral Sedimentation in a Bahamian Karst Basin (Sinkhole) over the Late Holocene: Influence of Local Vegetation and Climate. Palaeogeogr. Palaeoclimatol. Palaeoecol. 506, 70–83. doi:10.1016/j.palaeo.2018.06.014

CrossRef Full Text | Google Scholar

Toomey, M. R., Curry, W. B., Donnelly, J. P., and van Hengstum, P. J. (2013). Reconstructing 7000 Years of North Atlantic Hurricane Variability Using Deep-Sea Sediment Cores from the Western Great Bahama Bank. Paleoceanography 28, 31–41. doi:10.1002/palo.20012

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Cresswell, J. N., Milne, G. A., and Iliffe, T. M. (2019). Development of Anchialine Cave Habitats and Karst Subterranean Estuaries since the Last Ice Age. Scientific Rep. 9, 1–10. doi:10.1038/s41598-019-48058-8

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Donnelly, J. P., Fall, P. L., Toomey, M. R., Albury, N. A., and Kakuk, B. (2016). The Intertropical Convergence Zone Modulates Intense Hurricane Strikes on the Western North Atlantic Margin. Scientific Rep. 6, 1–10. doi:10.1038/srep21728

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Reinhardt, E. G., Beddows, P. A., and Gabriel, J. J. (2010). Linkages between Holocene Paleoclimate and Paleohydrogeology Preserved in a Yucatan Underwater Cave. Quat. Sci. Rev. 29, 2788–2798. doi:10.1016/j.quascirev.2010.06.034

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Reinhardt, E. G., Beddows, P. A., Huang, R. J., and Gabriel, J. J. (2008). Thecamoebians (Testate Amoebae) And Foraminifera From Three Anchialine Cenotes in Mexico: Low Salinity (1.5-4.5 Psu) Faunal Transitions. J. Foraminiferal Res. 38 (4), 305–317. doi:10.2113/gsjfr.38.4.305

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Reinhardt, E. G., Beddows, P. A., Schwarcz, H. P., and Gabriel, J. J. (2009). Foraminifera and Testate Amoebae (Thecamoebians) in an Anchialine Cave: Surface Distributions from Aktun Ha (Carwash) Cave System, Mexico. Limnol. Oceanogr. 54, 391–396. doi:10.4319/lo.2009.54.1.0391

CrossRef Full Text | Google Scholar

van Hengstum, P. J., and Scott, D. B. (2011). Ecology of Foraminifera and Habitat Variability in an Underwater Cave: Distinguishing Anchialine versus Submarine Cave Environments. J. Foraminiferal Res. 41 (3), 201–229. doi:10.2113/gsjfr.41.3.201

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Scott, D. B., Gröcke, D. R., and Charette, M. A. (2011). Sea Level Controls Sedimentation and Environments in Coastal Caves and Sinkholes. Mar. Geol. 286, 35–50. doi:10.1016/j.margeo.2011.05.004

CrossRef Full Text | Google Scholar

van Hengstum, P. J., and Scott, D. B. (2012). Sea-level Rise and Coastal Circulation Controlled Holocene Groundwater Development in Bermuda and Caused a Meteoric Lens to Collapse 1600years Ago. Mar. Micropaleontol. 90-91, 29–43. doi:10.1016/j.marmicro.2012.02.007

CrossRef Full Text | Google Scholar

van Hengstum, P. J., Winkler, T. S., Tamalavage, A. E., Sullivan, R. M., Little, S. N., MacDonald, D., et al. (2020). Holocene Sedimentation in a Blue Hole Surrounded by Carbonate Tidal Flats in the Bahamas: Autogenic versus Allogenic Processes. Mar. Geol. 419, 106051. doi:10.1016/j.margeo.2019.106051

CrossRef Full Text | Google Scholar

Wanless, H. R. (1981). Fining-upwards Sedimentary Sequences Generated in Seagrass Beds. J. Sediment. Res. 51, 445–454. doi:10.1306/212F7CA2-2B24-11D7-8648000102C1865D

Google Scholar

Whitaker, D., and Christman, M. (2014). Clustsig: Significant Cluster Analysis. Retrieved from https://CRAN.R-project.org/package=clustsig (Accessed January 24, 2021). doi:10.4324/9781315855226

CrossRef Full Text | Google Scholar

Wickham, H. (2016). Ggplot2: Elegant Graphics for Data Analysis. Springer.

Winkler, T. S., van Hengstum, P. J., Horgan, M. C., Donnelly, J. P., and Reibenspies, J. H. (2016). Detrital Cave Sediments Record Late Quaternary Hydrologic and Climatic Variability in Northwestern Florida, USA. Sediment. Geol. 335, 51–65. doi:10.1016/j.sedgeo.2016.01.022

CrossRef Full Text | Google Scholar

Young, M. B., Gonneea, M. E., Fong, D. A., Moore, W. S., Herrera-Silveira, J., and Paytan, A. (2008). Characterizing Sources of Groundwater to a Tropical Coastal Lagoon in a Karstic Area Using Radium Isotopes and Water Chemistry. Mar. Chem. 109, 377–394. doi:10.1016/j.marchem.2007.07.010

CrossRef Full Text | Google Scholar