Habitat Partitioning in the Marine Sector of Karst Subterranean Estuaries and Bermuda's Marine Caves: Benthic Foraminiferal Evidence

Introduction

Karst subterranean estuaries (KSEs) are an environment created from the two- and three-way mixing of saline groundwater, rain, and oceanic water in the subsurface on carbonate landscapes, and this hydrographic framework promotes unique physical processes, biogeochemical cycling, and biological communities. Research on subterranean estuaries has significantly increased in the last 20 years (Moore, 1999; Gonneea et al., 2014; Brankovits et al., 2017, 2018; Pain et al., 2019, 2020; Duque et al., 2020). Traditionally, the biologic-focused study (fauna and ecosystems) of caves that are flooded by a KSE have been described with the adjectives anchialine or marine cave environments (Iliffe et al., 1983, 1984; Stock et al., 1986; Bishop et al., 2015). Since the groundwater-ocean dynamics that create KSEs are linked to sea level, KSEs also vertically migrate in the subsurface in response to orbitally-forced changes in eustatic sea level (van Hengstum et al., 2019) and this is represented in the fossil record (Proctor and Smart, 1991; van Hengstum et al., 2009; Rosso et al., 2015, 2018; Guido et al., 2017). Despite advances in the last decades, there still is considerable uncertainty regarding ecosystem functioning in modern KSEs, including the impact of organic carbon on benthic community ecology (Figure 1).

Figure 1. Conceptual model of a karst subterranean estuary (KSE) and the anchialine habitat continuum created by subsurface groundwater variability (van Hengstum et al., 2019). Western North Atlantic anchialine fauna (e.g., fish, shrimps, decapods, ostracodes, foraminifera, marine Porifera, and Bryozoa) are positioned in their typical groundwater habitat based on observations from several localities, including The Bahamas (⋆), Yucatan Peninsula (•), and Bermuda (♦). Typhlatya shrimp exemplify habitat and resource partitioning in the KSE (Chávez-Solís et al., 2020).

Benthic foraminifera are useful for investigating environmental variability and benthic habitat partitioning in KSEs in both modern and ancient settings (Balduzzi and Cattaneo, 1985; Proctor and Smart, 1991; van Hengstum et al., 2008, 2009; Omori et al., 2010; Guido et al., 2017; Bergamin et al., 2018, 2020; Romano et al., 2018, 2020). Benthic foraminifera are eukaryotic protists that live among the meiofauna (0.062–0.5 mm in size) in most marine settings, and they exhibit ecological zonation from physical environmental gradients (e.g., salinity, organic matter flux, dissolved oxygen). In the deep-sea, benthic foraminifera they can account for 50% of eukaryotic biomass in sediment (Gooday et al., 1992). Foraminifera are heterotrophs that eat phytodetritus, terrestrial detritus, prokaryotes, among other trophic strategies, with only a few species living in symbiotic relationships with photosynthetic algae (Gooday, 1988, 1993; Goldstein and Corliss, 1994; Linke et al., 1995; Nomaki et al., 2005; Dupuy et al., 2010). Individual foraminifera produce a simple shell, or test, that remains within the sediment after death, so statistically significant populations can be sampled by SCUBA divers during routine ecological surveys, or their fossil remains can preserve a record of prehistoric KSEs.

Here we document how benthic foraminiferal communities correlate with particulate organic carbon source, quantity, and delivery mechanisms in Bermuda's horizontally expansive submerged caves, which represent areas of KSEs where salinity and dissolved oxygen are not limiting factors. Bermuda is a ~35 million year old volcanic seamount in the North Atlantic Ocean, that is capped by 30–80 m of limestone interspaced by paleosols (Pirsson and Vaughn, 1917; Woolard and Ewing, 1939; Plummer et al., 1976; Palmer et al., 1977; Vacher et al., 1995). The limestone units deposited during Quaternary sea-level high stands have weathered into a mature karst landscape (Mylroie et al., 1995), with abundant cave-related karst development in the oldest carbonate formations that are concentrated in the northeastern part of the island (Land et al., 1967; Plummer et al., 1976; Palmer et al., 1977; Vacher and Rowe, 1997). Bermuda is an ideal location for this study because the flooded caves are an established biodiversity hot spot of endemic anchialine fauna (Iliffe et al., 1983), and its surrounding marine environments have been the focus of long-term oceanographic study (e.g., Woolard and Ewing, 1939; Neumann, 1965; Morris et al., 1977) and foraminiferal work (Carmen, 1927; Steinker and Clem, 1984; Javaux and Scott, 2003). This study provides the broadest geographic perspective yet available for benthic foraminiferal distributions in the KSE habitat for any global island by combining foraminiferal census data from 28 new samples from the 2nd largest flooded cave (Palm Cave) into a new dataset with 100 samples from previously studied caves in Bermuda (Deep Blue, Cow Cave, Green Bay Cave).

Materials and Methods

Research Design and Study Sites

Submerged caves in Bermuda are most abundant in the oldest geologic formations (e.g., Walsingham Formation) between Castle Harbor and Harrington Sound, and between North Shore Lagoon and Harrington Sound (Figure 2). The volcanic-limestone contact has not yet been observed in any flooded cave, and cave passages are generally shallower than ~23 m below sea level. This study focuses on two large caves with physical connections into Harrington Sound (i.e., Green Bay Cave, Palm Cave), and two cave entrances that are slightly more inland (Cow Cave, Deep Blue). Harrington Sound is a shallow coastal embayment that is ~5 km², with a mean depth of 14.5 m below sea level, and a maximum depth of ~25 m below sea level in the southwest (Parsons et al., 2015). It is estimated that the seawater residence time in the mixed layer of Harrington Sound is ~30 days from tidal exchange of seawater through the narrow channel at Flatts Inlet (Morris et al., 1977), but effective tidal mixing can be as long as 4 months (Parsons et al., 2015). Harrington Sound experiences seasonal hypoxia or anoxia from water column stratification in summer (Coull, 1969; Bates, 2017), which causes benthic sedimentary and ecological zonation: a sandy shallow zone with sea grass and macroalgae (0–10 m), the Oculina coral zone (10–19 m), and the subthermocline zone (19–26 m) (Neumann, 1965). Seasonal dysoxia or anoxia have not been observed or documented in Bermuda's caves, and previous work documents that surficial foraminiferal distributions in Bermuda's caves are not typical of those observed in open ocean hypoxic to anoxic settings (Sen Gupta et al., 1996; Bernhard and Sen Gupta, 1997; Levin et al., 2009). This suggests that tidal pumping of seawater from the mixed layer in Harrington Sound into the caves is the dominant source of seawater. Tidal pumping also causes microtidal variability (<0.5 m) in the local groundwater table elevation (i.e., water table, Martin et al., 2012), which can be measured in the shallow pools that provide access to the flooded caves (Palmer et al., 1977; van Hengstum et al., 2015; Little and van Hengstum, 2019).

Figure 2. (A) Digital elevation model of Bermuda in the North Atlantic Ocean. (B) Detailed map of Harrington Sound plotted in ArcGIS. Red stars indicate study site locations.

Bermuda's caves offer a great location for investigating the marine sector of a KSE on a carbonate platform that is nearly flooded by eustatic sea-level rise. The commonly-studied flooded caves in Bermuda do not discharge a meteoric lens toward the ocean like caves on larger carbonate landscapes in the Mediterranean area (Radolović et al., 2015; Romano et al., 2020) or Yucatan Peninsula (Suárez-Morales et al., 2004; Chávez-Solís et al., 2020). This is in part because Bermuda is currently in a late stage of carbonate platform flooding relative to the position of global sea level (van Hengstum et al., 2019), and local meteoric water masses in Bermuda are spatially diminished. On the carbonate platforms of The Bahamas, for example, landmass area and mean annual rainfall alone can be used to predict meteoric lens thickness on different islands (Cant and Weech, 1986). Bermuda does have five modest meteoric lenses located in geologically-recent lithologies (Vacher, 1978), but these are not located in the primary area of cave development on the isthmus between Castle Harbor and Harrington Sound (Walsingham Formation, highest porosity: Land et al., 1967). In northeastern Bermuda where the Walsingham Formation is located, only a thin (<50 cm), seasonal and ephemeral brackish cap occurs in some cave-collapse entrances, and a limited meteoric lens can form near the “Letterbox” area of Green Bay Cave that does not directly contact the sediment-water interface, or discharge water through marine exit into Harrington Sound (see map in van Hengstum and Scott, 2011). The mineralogy of sediment in Bermuda's caves is dominated by carbonates, along with localized deposits of paleosol eroding into cave from the terrestrial surface above. Lastly, while large carbonate platforms in the Bahamas and the Yucatan frequently have large sinkholes that host extensive primary productivity and are interconnected to flooded caves (Schmitter-Soto et al., 2002; Tamalavage et al., 2018), the most comparable geomorphologic feature in Bermuda is Cliff Pool Sinkhole in the Green Bay Cave System. Even then, Cliff Pool still only supports a thin (<0.5 m) brackish water mass whose salinity depends upon seasonal rainfall and aridity.

Two of the largest cave systems in Bermuda, Palm Cave (primary entrance: 32.34°N, 64.71°W) and Green Bay Cave (primary entrance: 32.33°N, 64.74°W), were targeted for analysis of recent foraminiferal distributions (uppermost 5 cm of sediment) and compared to the impact of tidal variability on foraminiferal distributions at cave entrances to the Walsingham Cave System at Deep Blue (32.347°N, 64.711°W) and the entrance to Cow Cave (32.34°N, 64.71°W). Green Bay Cave also has the largest direct connection to the sea at Harrington Sound. According to measurements of Iliffe (2008), 1,960,000 m³ of seawater flows into Green Bay Cave at a cross sectional area of 25.5 m² during a flood tidal cycle. In contrast, Palm Cave has multiple, small subaerial entrances to the forest landscape in the Walsingham Nature Preserve, and one direct connection to Harrington Sound at Cripple Gate (cross sectional area: 0.6 m²) that permits ~1,000 m³ of seawater exchange during a flood tidal cycle (Iliffe, 2008). Both Palm Cave and Green Bay Cave have seagrass meadows (e.g., Thalassia testudinium, Syringodium filiforme) within 100 m of their marine entrances (Iliffe, 2008), in the shallow sandy benthic ecological zone in Harrington Sound (Neumann, 1965).

Benthic Foraminiferal Processing From Palm Cave

28 new surface (<5 cm depth, >30 cm³ of sediment) sediment samples were collected from various depths and locations throughout Palm Cave in August 2015 (Figure 3). At each sample station, a separate 50 mL centrifuge tube was filled with sediment by a SCUBA diver employing advanced technical cave diving procedures. Collecting this volume of sediment permitted separate sediment sub-samples from each station to be individually processed for benthic foraminifera, textural analysis, Loss-on-Ignition (LOI), and bulk organic geochemistry. The total assemblages of benthic foraminifera were analyzed at each sampling station, as some workers perceive this to be most reflective of long-term distributions (Buzas, 1968; Scott and Medioli, 1980; Debenay and Guillou, 2002). Benthic foraminifera were concentrated by wet sieving sediment sub-samples (0.63–5 cm³) over a standard 63 μm screen mesh, with the remaining coarse sediment residues sub-divided into representative aliquots using a wet splitter to enable representative census counts of ~300 individuals per sample (Patterson and Fishbein, 1989; Scott and Hermelin, 1993). While van Hengstum and Scott (2011) used a 45 μm sieve to concentrate benthic foraminifera in their study of Green Bay Cave, Little and van Hengstum (2019) documented negligible difference between resultant foraminiferal communities from Bermudian caves using 45 vs. 63 μm mesh sizes. Individual benthic foraminifera were wet picked onto micropaleontological slides and identified to the species level in most cases. Taxonomy was confirmed by comparing scanning electron microscopy of representative individuals to literature comparisons (Carman, 1933; Bermúdez, 1949; Loeblich and Tappan, 1987; van Hengstum and Scott, 2011; Little and van Hengstum, 2019). This produced a new dataset for benthic foraminifera from Palm Cave of 28 samples x 111 taxonomic observations (see Supplementary Material).

Figure 3. Detailed wall map of Palm Cave, showing explored cave passages that provide direct access to the subterranean estuary, with locations of surface sediment samples (orange circles) and cave entrances (adapted from original survey of Jason Richards). Maps and surface sample locations for Green Bay Cave, and Cow Cave and Deep Blue are provided elsewhere (van Hengstum and Scott, 2012; Little and van Hengstum, 2019).

Sediment Texture and Organic Geochemistry

Sediment texture (i.e., mean particle size) in Palm Cave was examined in separate sediment sub-samples using standard laser diffraction particle size analysis and loss on ignition (LOI) procedures. Organic matter was recorded as percent-lost relative to the original total sample weight after LOI at 550°C for 4.5 h (Heiri et al., 2001). Mean particle size for sediment sub-samples was measured with a Malvern Mastersizer 2000S laser particle size analyzer (Sperazza et al., 2004), with individual particles disaggregated in a sodium hexametaphosphate dispersant prior to analysis. The stable carbon isotopic value (δ¹³C_org) and C:N ratio of bulk sediment was measured on separate sediment sub-samples from the new sediment samples collected from Palm Cave (n = 28). In general, δ¹³C_org values of bulk sedimentary organic matter reflect the relative contributions of marine, terrestrial and pelagic organic carbon sources in coastal environments (Lamb et al., 2006; Khan et al., 2019). Terrestrial plants that preferentially use the C3 photosynthetic pathway produce organic tissues that are ¹³C-depleted and nitrogen-poor, with δ¹³C_org values generally ranging between −32 and −21‰ (Lamb et al., 2006). Sea grasses and macroalgae are common in the shallow waters from 0 to 10 m of Harrington Sound (sandy zone with seagrass and macroalgae of Neumann, 1965) where Palm Cave and Green Bay Cave have entrances (Wefer and Killingley, 1986; Fourqurean et al., 2015). Seagrasses are also important contributors of carbon isotopically enriched organic matter to the sediment in nearshore environments relative to marine plankton (Fry et al., 1977; Anderson et al., 1992; Kennedy et al., 2010). Based on measurements on plant tissues around Bermuda coastal waters, seagrasses (e.g., Syringodium and Thalassia: mean of −7.4‰, n = 455) and macro algae like Halimeda (mean of −14.3‰, n = 60) also have relatively more ¹³C-enriched δ¹³C_org values (see Table 1 in Burgett et al., 2018). The δ¹³C_org values of particulate organic carbon (POC) in the Sargasso Sea are similar to elsewhere in the tropical North Atlantic and range from ~-21 to −23‰ (Jeffrey et al., 1983; Altabet, 1990), and these values primarily reflect the ranges of marine planktonic algae (France, 1995). For samples from Palm Cave, carbonates were first digested with 1.0 M HCl for 12 h, and the remaining sediment residue was rinsed, desiccated, powdered, and weighed into silver capsules. Stable carbon isotope ratios (δ¹³C_org) were measured against international standards at Baylor University Stable Isotope Laboratory on a Thermo- Electron Delta V Advantage Isotope Ratio Mass Spectrometer. Final isotopic ratios are reported relative to the standard delta (δ) notation relative to Vienna Pee Dee Belemnite (VPDB) for carbon (expressed as parts per mil: ‰) with an uncertainty of ± 0.1‰, with the C:N ratio and total organic carbon (TOC, weight %) calculated from the analysis of the residual sediment fraction after acidification.

Data Analysis

A final dataset was compiled for further multivariate statistical analysis that contained 128 total samples: (a) the 28 new samples from Palm Cave (this study), 74 samples from Green Bay Cave (van Hengstum and Scott, 2011), 18 samples from Cow Cave (Little and van Hengstum, 2019), and 7 samples from Deep Blue (Little and van Hengstum, 2019). From Little and van Hengstum (2019), we use the >45 μm database, which has higher census counts and thus likely better represents the sampled areas. All taxonomic observations were grouped to the generic level to limit the impact of taxonomic bias imparted by different microscopists at the species level. The exception was the brackish water indicator Trochammina inflata, which was left at specific level to differentiate from co-generic marine taxa (e.g., T. globularis). Also, the final dataset amalgamated the taxonomically similar genera Pattelina, Heteropatellina, and Patellinoides (all members of the order Spirillinida). These genera were differentiated by Little and van Hengstum (2019, Deep Blue and Cow Cave), but lumped together by van Hengstum and Scott (2011, Green Bay Cave). The final combined dataset was 128 samples by 89 taxonomic units at the generic level (see Supplementary Material).

Numerical analyses were completed in R Studio software (version 1.1.383; R Studio, 2016). The final combined dataset of benthic foraminifera was used to calculate rarefied species richness (S) and the Shannon-Wiener Diversity Index (H), using the “rarefy” and “diversity” functions in the R package “vegan” (Gotelli and Colwell, 2011; Oksanen et al., 2013), and raw counts were converted to proportional abundances using the “decostand” function in the “vegan” package. Unconstrained Q-mode cluster analysis was used to help divide the dataset into ecologically meaningful benthic foraminiferal communities by grouping statistically similar populations (Legendre and Legendre, 1998). Prior to cluster analysis, the proportional abundances were square root transformed to minimize the impact of dominant species and to better compare community structure (Legendre and Legendre, 1998). Samples were then subjected to hierarchical cluster analysis using an unweighted paired group averaging algorithm and the Bray-Curtis dissimilarity index using the “simprof” function in the “clustsig” package (Clarke et al., 2008; Whitaker and Christman, 2014), and the “vegdist” function in the “vegan” package (Oksanen et al., 2013). The proportion of dominant taxa in the dendrogram (5% proportional abundance in at least 1 sample) were illustrated between the different samples by generating a vertical stratigraphic plot using the “stratplot” function in the “rioja” package (Juggins, 2015). Finally, physical environmental parameters were graphically summarized in box-and-whisker plots: δ¹³C_org, C:N, bulk organic matter-LOI, TOC, mean particle size, and mean groundwater depth, and compared to the benthic foraminiferal communities grouped by unconstrained Q-mode cluster analysis.

Results: Benthic Foraminiferal Communities

Similar to other Bermudian caves investigated (Green Bay Cave, Deep Blue, and Cow Cave), Palm Cave has abundant benthic foraminifera with a mean concentration of 7256 individuals per cm³ of sediment (range: 458–14,603). No new endemic species of foraminifera are described, but the identified benthic foraminiferal communities differ from those found in other Bermudian coastal environments (e.g., mangroves, carbonate lagoons, or reefs: see Javaux and Scott, 2003). Seven communities of benthic foraminifera can be recognized at a dissimilarity index of 0.44 on the dendrogram produced from unconstrained Q-mode cluster analysis (Figure 4), but this division increases to nine groupings if Community 5 is subdivided at a dissimilarity of index of ~0.36 (discussed below). Two samples were plotted at the bottom of the dendrogram and were deemed most dissimilar to the other samples: ISC-72 and CC5-ITA. From the original dataset, however, the very high relative abundance of Spirophthalmidium emaciatum (70%) in ISC-72 suggests this sample is most closely related to Community 6. Similarly, CC5-ITA had a very high abundance of Trochammina, Entzia, and Miliammina relative to other samples in Community 1. With currently available data, we consider that the samples CC5-ITA and ISC-72 are just endmember examples in Community 6 and 1, respectively.

Figure 4. Benthic foraminiferal communities in Palm Cave, Green Bay, Cow Cave, and Deep Blue as defined by groupings in a dendrogram produced from Q-mode cluster analysis. A 7-cluster (dissimilarity 0.44) or a 9-cluster division (dissimilarity 0.36) of the dendrogram correlates well with the measured environmental variables. Only taxa with >5% in at least one sample are shown to emphasize significant fauna.

Community 1 characterizes the intertidal zone of Cow Cave, and this lower diversity (H = 2.2) community is dominated by taxa known to tolerate periodic aerial exposure related to tides, high quantities of bulk organic matter in the sediment, and lower salinities. Dominant taxa include Trochammina inflata (mean: 18.3%), Miliammina fusca (mean: 13.1%), and Entzia macrescens (mean: 8.7%; Table 1). Tidal variability of the water table causes periodic exposure of these shallow benthic areas (Figure 5), and this community likely also experiences repeated brackish conditions from rain and surface runoff (Little and van Hengstum, 2019). Community 1 is characterized by the highest bulk organic matter content (mean: 25.5%), which is predominantly derived from terrestrial sources based on the mean δ¹³C_org values and C:N ratio (mean: −25.9‰, mean C:N ratio: 18.3; Figure 6).

Table 1. Selected summary statistics for environmental characteristics, and arithmetic mean of the percent abundance of dominant taxa in each benthic foraminiferal community.

Figure 5. Box-and-whisker plot of summary statistics (mean is red dot, median is central line) for benthic environmental characteristics associated with each benthic foraminiferal community. Despite the trend or gradient observed in δ¹³C_org values, C:N ratios, and bulk organic matter content, the brown (blue) colored boxes identify benthic foraminiferal communities that prefer terrestrial (marine) organic matter sources delivered to the benthos.

Figure 6. Crossplot of δ¹³C_org values and C:N ratio of bulk sediment from the Palm Cave (this work), Green Bay Cave (van Hengstum and Scott, 2011), Cow Cave, and Deep Blue (Little and van Hengstum, 2019). Samples are grouped according to their benthic foraminiferal communities. Colored ellipses correspond to 95% confidence intervals around the mean values for selected groups. The values of δ¹³C_org and the C:N ratio of other organic carbon pools are compiled from the literature (Lamb et al., 2006; Kennedy et al., 2010; Fourqurean et al., 2019).

Community 2 includes subtidal areas of Deep Blue in the saline groundwater mass that are exposed to normal sunlight. The benthos has lower bulk organic matter (mean: 18.8%) than Community 1, but the organic carbon remains primarily sourced from the adjacent terrestrial surface (mean: −26.1‰, C:N ratio mean: 11.1; Table 1). This community has higher species richness (S = 28) and diversity (H = 2.7) than Community 1, which is most likely because these benthic area are constantly submerged. This community is dominated by the infaunal genus Bolivina (mean: 16.6%), T. inflata (mean: 12.2%), and H. anderseni (mean: 7.5%, Table 1). Rosalina is also dominant (mean: 10.5%), which is primarily Rosalina globularis, and small-sized Triloculina oblonga (mean: 12.2, Supplementary Material).

Communities 3 and 4 occur in a transect of increasing depth into the saline groundwater at Cliff Pool Sinkhole within the Green Bay Cave (AC-04 to AC-11 of van Hengstum and Scott, 2011), and include the shallowest samples from Cliff Pool Sinkhole (ML-1, ML-2, ML3). The shallow sinkhole itself receives direct sunlight, and shallow sampling sites have taxa that must be tolerant of brackish conditions. Communities 3 and 4 start in open water, and form a transect down into the aphotic zone of Green Bay Cave. Community 3 is dominated by Bolivina spp. (mean: 20.2%), H. anderseni (mean: 21.4%), and Rosalina spp. (mean: 11%). While Ammonia (mean: 11.6%) has been grouped to the generic level, van Hengstum and Scott (2011) observed that Community 3 was characterized by Ammonia beccarii var. tepida. In contrast, Community 4 has higher species richness (S = 30) and diversity (H = 2.6), and Bolivina spp. (mean: 26%) and Rosalina spp. (mean: 11%) remain dominant (Figure 4; Table 1). However, Community 4 also includes other taxa like phytodetritus feeding Melonis barleeanum (mean: 5.8%), Siphogenerina columellaris (mean: 4.9%), Siphonina (mean: 4.2%), and Nonionella (mean: 4.7%, Figure 4). Community 3 contains a mean of 14.1% of bulk organic matter, and a mean δ¹³C_org value of −24.2‰ and a mean C:N ratio of 11 (Table 1), whereas Community 5 has higher bulk organic matter (mean 17.5%) content.

Most of the samples collected from throughout Palm Cave, Cow Cave, and Green Bay Cave collectively form Community 5 at a Bray-Curtis dissimilarity index of 0.44. Community 5 is located in benthic areas that: (i) are only impacted by saline groundwater, (ii) are not subaerially exposed anytime during the tidal cycle, and (iii) are in the aphotic zone. While all the areas characterized by Community 5 samples have similarities, important differences emerge with the physical environmental characteristics that warrant dividing further at a Bray-Curtis dissimilarity index of 0.36 (Figure 4). For example, Communities 5A, 5B, and 5C show a trend toward more enriched δ¹³C_org values and decreasing bulk organic matter content (Figure 5). Community 5A is located in the subtidal areas of the roofed Cow Cave, and is located in shallow subtidal pools that are proximal to the adjacent terrestrial landscape. Except for the samples collected specifically in Cliff Pool Sinkhole (ML-1, ML-2, and ML-3), most of Communities 5B and 5C are located deeper in the saline groundwater. Community 5A is also a transition between benthic foraminiferal communities characterized by more depleted δ¹³C_org values (Communities 1–4) vs. more enriched δ¹³C_org values (Communities 5B, 5C, 6, and 7), and similarly with C:N ratios (Figure 5). Despite a general linear trend in the δ¹³C_org values of the sediment (Figure 5), there is a distinct change in the benthic foraminiferal communities from the dominance of Bolivina and Rosalina in Community 5A to Spirophthalmidium, Spirillina, and Pattelina-like genera (Patellina, Heteropatellina, Patellinoides) dominating Communities 5B and 5C (Figure 4; Table 1).

Community 6 is located in areas most distal to any cave exit (either terrestrial or marine). Relative to most other areas, Community 6 is characterized by enriched δ¹³C_org values (mean: −18.5‰) and lower C:N ratios (mean: 8.6). The sediment where Community 6 is found has the lowest total organic carbon content (mean: 5.1 ± 2.2%) and bulk organic matter observed (mean: 6 ± 2.3%). While all cave areas are dominated by carbonate and organic matter, Community 6 is dominated by fine-grained carbonate particles (fine silt: 5.9 ± 1.3 μm, Table 1). Community 6 is less diverse than Community 5 (S = 22, H = 2.1), and is dominated by Spirophthalmidium emaciatum (mean 20.4%), Spirillina vivipara (mean: 19.6%), Patellina (mean: 9.6%), and Bolivina (mean: 4.3%).

Community 7 occurs in a transect from open water in the shallow (0–10 m) sandy ecological zone in Harrington Sound (E-47, E-48, E-49, and E-50), through the marine cave entrance into Green Bay Cave and stops at the Rat Trap area (see Figure 2 in van Hengstum and Scott, 2011). This community is characterized by lower bulk organic matter on average (mean: 10%), but high variability (range: 3.2–20%) depending upon the spatial variability in carbonate sand content (mean particle size: 311 μm, Table 1). The δ¹³C_org values (mean: −18.6‰) of the sediment are more enriched relative to Communities 1–4, with a lower C:N ratio (mean: 8.1). This community of benthic foraminifera has a higher diversity (H = 2.5) and species richness (S = 29), and is dominated by common carbonate lagoon genera Quinqueloculina (mean: 26%), but also Bolivina (mean: 12%), and Triloculina (mean: 11.4%).While Ammonia (mean: 9.7%) has been grouped to the generic level, van Hengstum and Scott (2011) observed that Community 7 was characterized by Ammonia beccarii var. parkinsoniana.

Discussion: Drivers of Foraminiferal Distributions

Our previous work (i.e., van Hengstum and Scott, 2011) divided the benthic foraminiferal communities in Green Bay Cave into those that were termed anchialine (terrestrially-dominated) vs. submarine (marine-dominated). This was based on the ideas of “terrestrial vs. marine influence” that were originally presented in Stock et al. (1986, more specifically), see the footnotes on page 91). This partly reflects the challenges faced by early workers with evaluating environmental and biological differences between: marine flooded caves like those found in northeastern Bermuda that have negligible meteoric lenses, vs. (ii) those with extensive meteoric lens development such as observed in caves the Yucatan Peninsula. However, a “terrestrial vs. marine” distinction is too simple for describing both benthic and pelagic habitat partitioning in subterranean estuaries, and the dependance of physical processes, biological communities, and biogeochemical cycles on specific groundwater masses (Brankovits et al., 2017, 2018; Chávez-Solís et al., 2020) that must vertically migrate through geologic time (Riedl and Ozretić, 1969; van Hengstum et al., 2019). For example, recent work with Typhlatya shrimp surveyed over a large area in the Yucatan Peninsula clearly documents pelagic habitat partitioning with respect to groundwater masses and carbon sources (see Figure 2 in Chávez-Solís et al., 2020). After applying broader environmental concepts of a KSE (Figure 1), the sampled marine caves of Bermuda include only the marine sector of KSEs that are not oxygen or salinity limited. As discussed below, modern benthic foraminiferal communities in the marine sector of KSEs in Bermuda's can be correlated with exposure caused by tidal variability of the water column, and the quantity, source, and transport mechanisms of particulate organic matter that is delivered to the benthos.

Marine Benthic Habitats Dominated by Terrestrial-Sourced Organic Carbon

In general, benthic foraminiferal Community 1 (mean: −25.9‰), 2 (mean: −26.1‰), 3 (mean: −24.2‰), and 4 (mean: −23.2‰) all occupy sediment that has more depleted δ¹³C_org values and lower C:N ratios (Figures 5, 6). Based on a comparison with literature values, this likely indicates that terrestrial organic matter is the dominant source of organic carbon delivered to these areas. In opposite, Communities 5, 6, and 7 occupied sediment with more enriched δ¹³C_org values (mean: −21.6 to −17.6‰, Table 1), which suggests that marine carbon sources are more dominant. The proximity between sample stations and cave openings (i.e., exits) to the forest landscape is clearly impacting organic matter delivery to the benthos. Detrital carbon is a key resource for benthic foraminifera in marine environments (Gooday, 1993; Jorissen et al., 1995), so it is unsurprising that foraminifera are also responding to different carbon sources in flooded marine caves. Cave exits provide a point source for terrestrial sedimentation into an otherwise marine environment, and diverse organic carbon pools increase the opportunity for resource partitioning and foraminiferal diversity. While terrestrial detritus eroding into the caves from the adjacent landscapes is likely a dominant carbon source for Communities 1–4 (δ¹³C_org values: −26.1 to −23.2‰), we cannot rule out the contribution of organic matter from marine plankton in the water column at Cliff Pool Sinkhole (Green Bay Cave) as also impacting the proximal Communities 3 and 4. Since marine POC from marine plankton has a δ¹³C_org value of −21 to −23‰ (Jeffrey et al., 1983; Altabet, 1990; France, 1995), the δ¹³C_org value of benthic sediment alone cannot differentiate the relative contributions of carbon from terrestrial detritus vs. marine plankton. Deposition of POC from marine plankton can also explain the modestly lower and stabilized C:N values in Communities 3 and 4 relative to Communities 1 and 2 (Figure 5). This is because vascular terrestrial detritus typically has much higher C:N ratios (Lamb et al., 2006), which is likely to collect in cave pools closer to the forest, whereas the deeper water column at Cliff Pool Sinkhole can support higher primary productivity. The uncertainty in these results is an opportunity to use more specific organic biomarkers to resolve carbon source contributions to the benthos. Overall, the more carbon isotopically depleted samples suggest an increased proportion of terrestrial detritus in Communities 1–4 relative to Communities 5–7, as all areas likely receive inputs of POC from marine plankton in response to active seawater-groundwater circulation between Harrington Sound and the caves (discussed below).

Similar to the result of Little and van Hengstum (2019), benthic habitats dominated by terrestrial organic carbon deposition can be further separated by tidal exposure: Community 1 is exposed to tidal variability of the water column, whereas Communities 2, 3, and 4 are all subtidal (Figure 5). The micro-tidal changes to the water table elevation at Cow Cave are sufficient to expose Community 1 during low tides (Figure 5, Little and van Hengstum, 2019). Community 1 is dominated by agglutinated, brackish-tolerant species Trochammina inflata, Miliammina fusca, Entzia macrescens, and Helenina anderseni (Table 1). These taxa are common in marshes, which also experience tidal exposure (Hayward et al., 1999; Javaux and Scott, 2003; Camacho et al., 2015). Community 1 likely also experiences intermittent brackish conditions from runoff during pluvial events given the close proximity of Community 1 to the terrestrial landscape (Figure 7). Taxa like E. macrescens and Miliammina fusca are useful sea-level indicators due to their significant correlation with tidal elevation (Scott and Medioli, 1978; Guilbault et al., 1995; Gehrels and van de Plassche, 1999). In the studied Bermuda caves, Miliammina fusca is the primary taxon in intertidal areas dominated by terrestrial organic carbon, with E. macrescens occurring in both intertidal (Community 1) and subtidal (Community 2) areas (Figure 4).

Figure 7. Conceptual model of idealized benthic foraminiferal communities in the marine sector of karst subterranean estuaries in Bermuda, along with idealized changes in stable carbon isotopic value of bulk organic matter (δ¹³C_org), organic matter quantity, mean particle size of substrate, tidal exposure in subaerial pools, and the magnitude of nutrient and particle delivery from ocean-groundwater exchange through marine cave openings. Representative values of δ¹³C_org values for source organic matter summarized from the literature (Goericke and Fry, 1994; Kennedy et al., 2010; Fourqurean et al., 2019), and ecological zones in Harrington Sound from Neumann (1965).

Subtidal benthic areas dominated by the deposition of terrestrially-sourced organic carbon could be further divided by light exposure, depth into the subsurface, and particle size (Communities 2, 3, and 4, Figure 7). None of the recovered benthic foraminifera in the caves employ photosymbionts (e.g., Amphistegina), therefore, foraminifera are all feeding on materials delivered to, or otherwise available in, the substrate (e.g., algae, bacteria, detritus, or dissolved organic matter). While these communities are all dominated by Rosalina and Bolivina (Table 1), differences in the secondary taxa are relevant to the community structure. It is important to note that Cliff Pool Sinkhole has a water column open to sunlight, and eyewitnesses report that the sinkhole can appear seasonally greenish from primary productivity. This potentially explains faunal differences between Cliff Pool Sinkhole (increased abundance of Trichohyalus aguayoi and Helenina in Community 3) vs. Deep Blue (increased abundance of Trochammina and Entzia in Community 2, Figure 4), with the potential for increased resource partitioning and foraminiferal diversity deeper into the saline groundwater mass at Cliff Pool Sinkhole (Communities 3 and 4). Other secondary taxa that are also present include Melonis barleeanum (Caralp, 1989) and Globocassidulina subglobosa, both of which are known to respond favorably to seasonal delivery of phytodetritus in the open ocean. Their increased percent abundance in Communities 3 and 4 proximal to Cliff Pool Sinkhole is likely related this additional food resource produced in the dimly-lit water column, and is supported by modestly lower C:N ratios (Figure 5). This likely indicates that additional environmental variables are impacting the benthic foraminiferal communities that are not quantified in this study (i.e., dissolved nutrients at the sediment-groundwater interface, pore water chemistry). Based on experimental cultures, the epifaunal Rosalina globularis prefers warm and saltier waters (optimal: 27°C, >25 psu) for reproduction (Saraswat et al., 2011, 2015), which are found in the relatively shallow (<20 m depth) marine habitats in Bermuda's caves. The most abundant genus is the infaunal Bolivina, which is likely capitalizing on any available food resource in the benthos. While the impact of particle size on benthic foraminifera in other environments remains a debated issue (Debenay et al., 2001; Diz et al., 2004; Armynot du Châtelet et al., 2009), our results suggest that substrate texture in KSEs may also be influencing benthic habitat variability and foraminiferal community structure. Community 3 was associated with coarser substrates (mean particle size: >271.2 μm) compared with Community 4 (mean particle size: 149.6 μm).

These results now clearly inform the environmental conditions of Pleistocene-aged KSEs that were colonized by benthic foraminifera, which have been subsequently preserved in fossil cave deposits. For example, Proctor and Smart (1991) documented fossil benthic foraminifera in elevated caves in England, which are most likely dated Marine Isotope Stage 11 (424,00–374,000 years ago). Detailed foraminiferal counts were not provided, but fine textured sediment contained echinoid spines, Ammonia beccarii, Cassidulina, and Bolivina. On the Yucatan Peninsula in Mexico, van Hengstum et al. (2009) described Bolivina-dominated foraminiferal communities in Carwash Cave (older than Last Glacial Maximum 20,000 years ago), and there are Helenina- and Bolivina-dominated Pleistocene-aged foraminiferal deposits in Maya Blue Cave and Vaca-Ha (van Hengstum and E. Reinhardt, unpublished observations). Based on modern benthic foraminiferal communities in Bermuda, these examples of fossil benthic foraminifera were all living in sediment that was dominated by terrestrial organic matter inputs, with additional taxonomic variability in the Pleistocene foraminiferal community most likely related to the salinity of the contemporaneous groundwater mass flooding that specific cave.

Marine Benthic Habitats Dominated by Marine-Sourced Organic Carbon

Communities 5, 6, and 7 occur in benthic habitats that primarily receive marine sources of organic carbon, which is most likely from marine plankton based on the δ¹³C_org values from −22.3 to −15.1‰ and C:N ratios of 8–10 (Figures 6, 7). Marine particulate organic carbon generally has values from −18 to −23‰ in lower latitudes (Jeffrey et al., 1983; Altabet, 1990; Goericke and Fry, 1994). An important component of particulate organic matter is marine plankton, which has a much lower C:N ratio (<10) relative to terrestrial organic carbon (Lamb et al., 2006; Figure 6). Seagrasses and macroalgae are abundant in the shallow sandy ecological zone of Harrington Sound (upper 10 m; Neumann, 1965), and consistent with global trends (Kennedy et al., 2010), seagrasses in Harrington Sound have enriched δ¹³C_org values and C:N ratios exceeding 18 (Fourqurean et al., 2019). Two locally common seagrass species, for example, Thalassia testudium has a mean δ¹³C_org value of −7.5‰ (n = 351) and a mean C:N ratio of 21 (n = 364), and Syringodium filiforme has a mean δ¹³C_org value of −6.3‰ (n = 369) and a mean C:N ratio of 18 (n = 362) (Fourqurean et al., 2019; Figure 6). Sediment below seagrass in Harrington Sound is estimated to have a C:N ratio of 11, and δ¹³C_org value of −13 to −16‰, based on a global mean Δ¹³C_{seagrass−sediment} value of 6.3‰ (Kennedy et al., 2010). While some contribution of organic matter from seagrasses and macroalgae cannot be ignored, the low C:N ratios associated with Communities 5, 6, and 7 (~8–10, Figure 7; Table 1) indicates that a more nitrogen-rich organic matter source like marine plankton is likely being delivered to these communities.

Community 7 is located in the main tunnel in Green Bay Cave that is closest to Harrington Sound, and extends from open water into an area colloquially referred to as the “Rat Trap” (see map in van Hengstum and Scott, 2011). This tunnel is characterized by attenuating light availability (disphotic zone), and past the Rat Trap is the aphotic zone. This tunnel is dominated by Quinqueloculina, Ammonia, Elphidium, and Peneroplis (Figure 6), which are common taxa in the shallow carbonate lagoons throughout the tropical North Atlantic Ocean (Javaux and Scott, 2003). While this study grouped all taxa at the genus level, van Hengstum and Scott (2011) noted that Community 7 is associated with Ammonia beccarii var. parkinsoniana. This tunnel coincides with areas that receive classic taxonomic investigation in marine caves, where strong gradients in light and tidal currents exist. The coarser substrate is likely promoted by winnowing-away of finer sediment particles in response to the large volume of seawater circulating through this tunnel during each tidal cycle (~1,960,000 m³: Iliffe, 2008). Coarser substrates also occur in other Mediterranean marine caves with higher current velocities from tidal action (Radolović et al., 2015). It is likely that tidal and storm currents also have transported some benthic foraminifera tests from Harrington Sound into the cave tunnel, however, current velocities notably decrease through the rest of Green Bay Cave past the narrowing of cave passage at the Rat Trap. Interestingly, PCS-15 from Palm Cave was not similar to Community 7, despite its proximity to the opening with Harrington Sound (i.e., Cripple Gate). Perhaps this is related to the much smaller connection of Palm Cave to Harrington Sound at “Cripple Gate,” and the much lower volume of water exchanged with Harrington Sound during tidal cycles impacting benthic habitat conditions.

Community 6 is most distal from any cave entrance (or exit, Figure 7), is associated with areas where endemic anchialine cave animals (e.g., Procaris chacei) are more likely to be observed, and the fine-grained substrate (mean particle size: 5.3 μm) has the lowest bulk organic matter content (mean: 6%) and TOC (mean: 5%). Community 6 has the highest percent abundance of Spirophthalmidium emaciatum, members of the order Sprillinida (e.g., Patellina, Spirillina, and Mychostomina), and Rotaliella arctica. Observation of Spirophthalmidium emaciatum outside of marine caves is limited, and includes New Guinea (Haig, 1988), the Mediterranean (Kuhnt et al., 2007), and the deep Pacific (Brady, 1884). In general, the thin-walled and fragile test of S. emaciatum is unlikely to endure taphonomic reworking. This taxon is very closely related (or a ecophenotypic variant) of S. acutimargo, which constitutes up to 3% of communities in deep-marine carbonate oozes (>2000 m water depth) from the North Atlantic Ocean near the Azores (Hermelin and Scott, 1985). While common in Community 6 (mean water depth: 15.4 m), Rotaliella arctica is also observed at 4450 m water depth on the Bermuda Rise (Pawlowski, 1991), and at 1980 m water depth in the Arctic (Scott and Vilks, 1991). Interestingly, dominant taxa in Community 6 (Rotaliella, Spirillinids, and Spirophthalmidium) are not common in carbonate sediments from shallower locations (300–1000 m) off the Little Bahama Bank in The Bahamas (Martin, 1988), or on the shallow (10–40 m) flooded carbonate Serranilla Bank (Triffleman et al., 1991) elsewhere in the Tropical North Atlantic Ocean. In caves, a deep-sea affinity is well-known in higher taxonomic groups like crustaceans (Iliffe et al., 1983). In the Mediterranean, S. emaciatum has been recently described in benthic foraminiferal communities from a marine cave with carbonate sedimentation on the Murica coastal region in Spain (Cueva Gemelas: Bergamin et al., 2020), but not in marine caves on eastern Sardinia with siliciclastic sedimentation and evidence for the discharge of meteoric water (Romano et al., 2018, 2020). On the island of Cozumel (Mexico) where Aerolito Cave is located (Mejía-Ortiz et al., 2007a,b; CalderóN-GutiéRrez et al., 2018), caves passages flooded with oxygenated saline groundwater and benthic carbonate sediment also host populations of S. emaciatum (van Hengstum, personal observation). Murray (2006) states that Spirillina and Patellina prefer hard substrates, so it remains possible that the cave walls/ceiling may be providing an ideal natural hard substrate for these taxa. However, our research design is unable to assess this possibility (sampling benthic sediment, not cave walls). As suggested by Fichez (1991), perhaps complex organic matter is a critical resource for Community 6, which can only be evaluated by resolving biogeochemical cycling of nutrients in these areas (Fichez, 1990, 1991; Airoldi and Cinelli, 1996; Jian et al., 2020). Finally, Community 6 also occurs very close to the shoreline in Green Bay Cave (samples ISC-61 or ISC-73 See fig 2 in van Hengstum and Scott, 2011), yet does not receive enough nutrients to promote higher diversity and species richness like Community 5 from direct seawater transfer through the limestone pores. From a conservation perspective, these results indicate that any introduced physical cave openings from anthropogenic activity will change nutrient supply and communities of benthic meiofauna in marine sectors of KSEs.

Community 5 (including: 5A, 5B, and 5C) is colonizing benthic habitats experiencing more regular tidal groundwater-seawater exchange and likely higher sedimentation rates (Community 5B and 5C) or some terrestrial input (Community 5A). This is likely providing increased nutrient availability, and in turn is promoting higher species richness and diversity than benthic habitats colonized by Community 6 (Table 1). Community 5 is not primarily influenced by: (i) strong tidal currents and attenuating light availability like Community 7, (ii) exposure from water table variability like Community 1, or (iii) increased organic matter input from the adjacent terrestrial landscape like Communities 2–4. Dominant taxa in Community 5 include Globocassidulina subglobosa that commonly occurs in carbonate sediments >1000 deep of the Little Bahama Bank (Martin, 1988), and Melonis barleeanum who is successful in areas with regular delivery of marine organic matter (Caralp, 1989). In Green Bay Cave, Community 5 is proximal to the large entrance (exit) into Harrington Sound. Tintinnid loricas (pelagic ciliates) are also deposited in the sediment with Community 5 in Green Bay Cave, and they were likely transported into the cave and settled out of the suspended sedimentary load. This observation supported the hypothesis that Community 5 in Green Bay Cave receives increased deposition of suspended sediments as the velocity of groundwater currents caused by tidal flow attenuate in the subsurface (van Hengstum and Scott, 2011), similar to what has been observed in marine caves on the Mediterranean coast of France (Fichez, 1991). An important taxonomic caveat is that small miliolids identified as Triloculina in Bermuda caves are not well-developed adult specimens with clear phenotypic characteristics like those common in subtidal marine carbonate lagoons (Javaux and Scott, 2003). Thus, it cannot be ruled-out that these are juvenile Quinqueloculina (Schnitker, 1967). Overall, these smaller cave miliolids are perhaps responding to a nutritional deficiency that is not constrained in this study.

Further evaluation of the physical environmental characteristics (Table 1; Figures 5–7) supports a tripartite division of the dendrogram (Figure 4, Bray-Curtis Similarity Index: 0.36). As previously discussed, the subtidal benthic foraminifera in Cow Cave (Community 5A, Cow Cave) is transitional between the foraminiferal communities that favor primary inputs from marine vs. terrestrial organic matter (Figure 5). While these samples are from shallow subtidal settings, benthic areas in Cow Cave receive protection from direct terrestrial sedimentation by the rock roof of the cave (i.e., geomorphology). Faunistically, samples from Palm Cave (Community 5C) closely compare with benthic areas in Green Bay Cave (Community 5B, Figure 4). However, finer substrates (Community 5B) had a greater abundance of Reophax spp. (mean: 3.1%) and Bolivina spp. (mean: 11.4%), while coarser substrates (Community 5C) had a greater abundance of Spirillina vivipara (mean: 13%) and Patellina-like taxa (mean: 7.8%). However, a careful evaluation of the cross plot between δ¹³C_org and C:N (Figure 6) indicates that organic matter in Palm Cave (Community 5C) is slightly more carbon isotopically depleted with higher C:N values than Green Bay Cave (Community 5B). At least two hypotheses emerge to explain these differences between Community 5C and 5B that are related to cave morphology: either (i) proportionally less carbon isotopically enriched organic matter (e.g., seagrass, algae) is transported into Palm Cave through the smaller opening at Cripple Gate vs. the large marine opening of Green Bay Cave (Figure 3), or (ii) more ¹³C-depleted terrestrial organic matter is eroding into Palm Cave through its more numerous openings (exits) to the terrestrial landscape above (Figure 2). This suggests that differences in tidal exchange of seawater and particulate organic carbon transport through Cripple Gate (Palm Cave) and Green Bay Entrance (Green Bay Cave) may be influencing benthic habitat partitioning between caves, but other geochemical proxies are required to resolve these uncertainties. The collective role of cave geometry and volumetric exchange of seawater during tidal cycles between saline groundwater and Harrington Sound appear important for generating differences in the benthic foraminiferal community structures, and by corollary, benthic habitats in the oxygenated marine sectors of KSEs.

Conclusions

• Benthic foraminifera in the marine sector of karst subterranean estuaries (KSEs) in Bermuda can be grouped into 9 different communities that are correlated to organic matter source and quantity, tidal exposure, texture of the carbonate sediment substrate (i.e., mean particle size size), proximity to cave exits, and the relationship between cave morphology and tidally-forced circulation between coastal seawater and saline groundwater.

• The source and quantity of particulate organic carbon delivered to the benthos is significantly correlated to benthic foraminiferal community structure. A primary division emerges in foraminiferal communities correlated with organic carbon inputs primarily from the terrestrial surface with some marine plankton (mean δ¹³C_org values <-23.2‰, C:N ratios >11), vs. organic carbon inputs dominantly from marine plankton (mean δ¹³C_org values >-21.6‰, C:N ratios <10) that are transported into the caves from groundwater-seawater circulation from the adjacent carbonate lagoon (i.e., Harrington Sound).

• Benthic areas with more depleted δ¹³C_org values of sediment are dominated by Rosalina, Bolivina, and Triloculina. In contrast, benthic areas with more enriched δ¹³C_org values of sediment are dominated by Patellina, Spirillina, Rotaliella, Miliolinella, and Spirophthalmidium.

• Benthic areas with the least amount of organic carbon in the sediment (mean: TOC 5.1%, mean bulk organic content: 6%) are the least diverse, and dominated by Spirophthalmidium emaciatum (mean: 31%, maximum: 70%).

• Differences in the size and morphology of direct connections between Harrington Sound and Palm Cave vs. Green Bay Cave are impacting benthic foraminiferal community structure.

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

PvH designed the research. JC collected the samples from Palm Cave and completed laboratory activities. PvH and JC jointly interpreted that data, generated the artwork, and co-wrote the final manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by student research grants from the American Academy of Underwater Sciences (JC, 2016), the Geologic Society of America (JC), the Bermuda Aquarium, Museum and Zoo (JC), and resources for publishing from Texas A&M University at Galveston.

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

We acknowledge the support of S. R. Smith, G. Nolan, B. Williams, with sample collections approved through the support of a Special Permit (2015, Number 150301) issued to PH from the Department of Environmental Protection in the Bermuda Government. This research is only possible through the long-term support from the Tucker family, the Bermuda National Trust, and members of the Bermuda Zoological Society. Additional technical support was provided by S. Little, A. Tamalavage, R. Sullivan, T. Winkler, T. Iliffe, and T. van Hengstum. A constructive evaluation by two reviewers and D. Brankovits improved the final manuscript.

Supplementary Material

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

References