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 (δ13Corg 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 δ13Corg 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.

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 (δ 13 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 δ 13 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 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 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.

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 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), groundwaterdependent habitats (Bauer-Gottwein et al., 2011), and other native ecosystems, like coral reefs and coastal vegetation (e.g., mangroves, seagrass) (Muckelbauer, 1990).

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 Losson-Ignition procedure (Sieve-first LOI) . 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 3 (1.25 cm 3 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 3 (expressed throughout as D >63µm mg cm −3 ). 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%. 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 , Carwash , Yax Chen (Collins et al., 2015), Maya Blue (Pohlman et al., 1997), and Cenote Bang (Brankovits et al., 2017).
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 (δ 13 C org ) and the C:N ratio. In coastal settings, δ 13 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 3 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 δ 13 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 δ 13 C values for a two-point calibration, are within 0.2‰ (1σ) δ 13 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 δ 13 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 δ 13 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. 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).
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 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: 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: where δ sediment is the measured δ 13 C org value of the sediment sample, δ terrestrial is the average δ 13 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 δ 13 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 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 (δ 13 C org ), and C:N ratios. Radiocarbon and microfossil sample positions are marked on relevant cores.
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 5 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 3 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 3 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 3 . 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 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.
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 δ 13 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), 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).

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 redissipation 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 −1 . 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 −1 . 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 −3 ) 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 −3 ± 40, whereas the rest of the 22 cm core has a mean D >63µm is 28 mg cm −3 ± 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 −3 ± 79 (n 163 cm) versus 67 mg cm −3 ± 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 δ 13 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 δ 13 C org values become more enriched (i.e., positive) ( Figure 3). As a result, the OM-limited sedimentary sections near the top are relatively 13 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 δ 13 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 δ 13 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 δ 13 C org values indicates that 13 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 δ 13 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).

Foraminiferal Communities
Benthic foraminifera were abundant in sediment from sampled areas in Aerolito (Figure 2). Two groupings can be identified in Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 the dendrogram produced from unconstrained Q-mode cluster analysis (Figure 6), which generally separates core-top samples (i.e., less bulk OM, more depleted δ 13 C org values) from the bottom of the cores (i.e., higher bulk OM, more enriched δ 13 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. 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) FIGURE 4 | Isotope mixing model results. Downcore trends in the percent terrestrial organic carbon (OC) as calculated from the δ 13 C org values using a simple twosource mixing model. The terrestrial endmember value is δ 13 C org −28.2‰, as the mangrove vegetation reported by Gonneea et al. (2004), while the marine endmember is δ 13 C org −20.1‰, as the marine POM along the coast of the Cozumel Channel reported by Brankovits et al. (2017).

FIGURE 5 | (A)
Cross-plot between stable carbon isotopic values (δ 13 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 δ 13 C org values.  Table S3). In contrast, the OMenriched 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 −3 ) than the OM-enriched Community (mean: 915 individuals cm −3 ), 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%).

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 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).
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 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 hurricaneinduced 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 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  and Green Bay Cave  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., 2019). We speculate that a similar mechanism is also occurring in 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 δ 13 C org values of sediment OM.
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 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 13 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 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 13 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 δ 13 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 δ 13 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 δ 13 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 δ 13 C org values) could also explain the decrease in the fraction of terrestrially derived OM in the sediments. However, the OMenriched sediment in Aerolito has depleted δ 13 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.
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 Collins et al., 2015). Higher OM content downstream is associated with high C:N and depleted δ 13 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-tosea 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 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.
Frontiers in Environmental Science | www.frontiersin.org May 2021 | Volume 9 | Article 670914 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;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 rainfallinduced 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, nonmiliolid 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 δ 13 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 δ 13 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 δ 13 C org values ( Figure 7B), which suggests that marine-derived POM delivery is beneficial to miliolid taxa in these environments. Sediment δ 13 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 Cresswell et al., 2017;Little and van Hengstum, 2019;Cresswell and van Hengstum, 2021), where habitats with more increased terrestrial-based OM (depleted δ 13 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 δ 13 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 suboxic 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 sealevel 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 δ 13 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 13 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 δ 13 C org values), the benthic foraminiferal community in Aerolito was composed of taxa preferentially found in sediment with 13 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.