Subsurface Microbial Community Composition in Anchialine Environments Is Influenced by Original Organic Carbon Source at Time of Deposition

Introduction

Coastal sediment contains an estimated 10³ to 10⁹ prokaryotes/gram, and these microorganisms play pivotal roles in the global biogeochemical cycling of carbon (Jorgensen, 2012; Orsi, 2018), sulfur (Holmkvist et al., 2011), and nitrogen (Kuypers et al., 2003). Microbial community diversity in sediment has been characterized by sequencing of the 16S ribosomal RNA (rRNA) gene in coastal marine environments (Koretsky et al., 2005; Lenk et al., 2012; Ul-Hasan et al., 2019), freshwater lakes (Nold et al., 2010; Kojima et al., 2012), and the deep-sea (Reed et al., 2002; Engelen et al., 2008; Inagaki et al., 2015; Labonté et al., 2017). Metagenomics and transcriptomics analyses showed that coastal sediment microbes are physiologically versatile and functionally diverse (Dyksma et al., 2016; Chen et al., 2017; Yu C. et al., 2018). However, the factors influencing microbial community composition and distribution in the sediment below anoxic water columns remain understudied, and are mostly unknown in anchialine environments.

The subsurface mixing of rain and marine water on carbonate platforms creates karst subterranean estuaries (van Hengstum et al., 2019; Brankovits and Pohlman, 2020; Cresswell and van Hengstum, 2021), which are analogous to other coastal estuaries by having an upper meteoric water mass of varying salinity buoyed on a denser saline groundwater mass below. The steep environmental gradients created by these variable water mass conditions create subsurface habitats and environments frequently prefaced with the adjective ‘anchialine’ (Stock et al., 1986), especially in the Atlantic realm. Blue holes and sinkholes not only provide access points for human hydrographic observation of karst subterranean estuaries (Smart et al., 1988), they also enable terrigenous chemicals, nutrients and organic matter to enter the otherwise subsurface anchialine environments (Collins et al., 2015b), and exposure to sunlight also promotes organic matter production by photosynthesis (Tamalavage et al., 2018). Sinkholes and blue holes are abundant across the Bahamian archipelago, and they are located in both subaerial (‘sinkhole’) and shallow subtidal settings (‘blue hole’). Structurally, Bahamian carbonate platforms began forming in the Jurassic from the long-term accumulation of shallow water carbonates (Mullins and Lynts, 1977), with sinkholes and blue holes subsequently created by subsurface dissolution and collapse of limestone bedrock (Mylroie et al., 1995). Owing to the high primary porosity of the antecedent carbonate and the present elevation of relative sea level, sinkholes on the terrestrial landscape are now aquatic environments as their depth extends below the water table and into the local groundwater. With sufficient local rainfall and aquifer recharge, groundwater on low-lying Bahamian islands is stratified with an upper meteoric water mass (typically oxic, lower salinity) divided from a lower saline groundwater water mass (anoxic, seawater), both of which destabilize and mix with increasing proximity to the shoreline. Both terrestrial organic matter coming from runoff and organic matter produced in situ (primary productivity) contribute to the particulate organic matter pool that is ultimately deposited within a sinkhole basin (van Hengstum et al., 2010; Collins et al., 2015a; Tamalavage et al., 2018; van Hengstum et al., 2018).

While the drivers of subsurface microbial communities remain poorly understood, several hypotheses have been proposed, which included varying sediment composition, enzymatic activity, or organic matter bioavailability. For example, evidence has been provided that subsurface microbial community structure in deep-sea anoxic sediment from the Juan de Fuca ridge may be caused by variations in sediment chemistry and texture (e.g., clays, carbonate-content, and organic carbon), as opposed to porewater geochemistry (Labonté et al., 2017). Subsurface microbial community composition has also been linked to enzymatic activity (Böer et al., 2009) and labile organic matter (Franco et al., 2007). Variations within microbial community structure in subtidal sands significantly correlated to enzymatic activities, specifically chitobiase, α-glucosidase, and phosphatase activity, as well as microphytobenthic algae productivity (Böer et al., 2009). The seasonal variation of chlorophyll a, a proxy for the availability of labile organic matter, slightly shifted the subsurface microbial community structure within sediment from the North Sea (Franco et al., 2007). In riparian settings, organic matter provenance between different terrestrial organic matter (OM) pools (e.g., soil, vegetation) was shown to influence the subsurface microbial community (Fagervold et al., 2014).

In anchialine environments, microbial communities in the water columns vary greatly between sinkholes (Sahl et al., 2011; Kajan et al., 2022), and they tend to follow oxygen concentration and density stratification (Gonzalez et al., 2011; Davis and Garey, 2018; He et al., 2019). Other than site and the type of watershed, salinity is one of the strongest drivers of the microbial community composition (Hoffman et al., 2018). Inferred microbial functions showed high potential for sulfur oxidation and nitrate reduction in sinkholes, with nitrite and sulfite reduction associated with the photic zone (He et al., 2019). Moreover, there are striking differences between the microbial communities found in the water column and those found in microbial mats and surface sediment (Sahl et al., 2011). Surface sediment microbial communities also depend on shoreline proximity (i.e., coastal or inland) (Suárez-Moo et al., 2022). However, little is known regarding sediment microbial communities in anchialine environments, or the factors shaping their composition and distribution. In this study, we aimed to explore the subsurface prokaryotic community in Blackwood Sinkhole on Abaco Island in the northern Bahamian Archipelago. Based on work introduced above, we hypothesize that changes in the structure and microbial composition will be mainly related to organic matter provenance at time of deposition. Aside from this work, terrestrial landscape changes and productivity changes in the water column at Blackwood Sinkhole are well documented in previous studies (van Hengstum et al., 2016; Tamalavage et al., 2018; Tamalavage et al., 2020; Fall et al., 2021). To test these ideas, we (a) described the composition and structure of the microbial communities using high-throughput sequencing the V4-V5 hypervariable regions of the 16S rRNA gene, and (b) compared to results to known changes in sediment and organic matter supply in Blackwood Sinkhole over the last ~1900 years that occurred from known evidence of local landscape and sinkhole productivity changes (Tamalavage et al., 2018; Tamalavage et al., 2020).

Study Site

Blackwood Sinkhole is located ~220 m from the shoreline on the northeastern coast of the Great Abaco Island on the Little Bahama Bank (26.79, 77.42) (Figure 1). Blackwood Sinkhole is 32 m in diameter, and this groundwater-fed basin has a stratified water column that has remained stable during the last 3000 years (van Hengstum et al., 2016; Tamalavage et al., 2018; Fall et al., 2021). The sinkhole has an upper ~10 m of meteoric water (1.4 psu), which transitions to a mixing zone below (10-15 m below sea level), and this rests above anoxic saline groundwater (~15-40 m below sea level, 39.9 psu) (van Hengstum et al., 2016). The groundwater stratification and geomorphology of Blackwood Sinkhole has promoted excellent sediment preservation from a lack of vertical sediment mixing from either invertebrate bioturbation or waves. Basal sedimentary deposits in Blackwood Sinkhole are carbonate (karst) gravel, which abruptly transition to laminated sapropel (i.e., high organic matter) interbedded with carbonate horizons that have accumulated over the last 3000 years (van Hengstum et al., 2016). The laminated sapropel has a total organic carbon (TOC) content ~10% (Tamalavage et al., 2018). However, sedimentation rate and site of primary deposition has not been laterally or temporally uniform through time. This is because (a) the bottom geometry of the sinkhole basin is complex since it was inherited from the collapse of an original limestone cave, and (b) the sediment sources and delivery mechanisms during the last 3000 years are linked to changing primary productivity, differential erosion of the vertical sinkhole walls, and organic matter transport and deposition from adjacent wetland and surrounding forest habitats (Tamalavage et al., 2018). Using detailed radiocarbon dating on other core samples, constant sedimentation (~0.4 mm yr⁻¹) was initiated 3000 Calibrated years Before Present (cal yrs BP) around the sinkhole periphery, with evidence of coarse-grained particle deposition during intense hurricane strikes on Abaco Island (van Hengstum et al., 2016). However, sedimentation in the center of the sinkhole was delayed until ~1900 cal yrs BP, with higher sedimentation rates (1.2 mm yr⁻¹) than the periphery and far less erosional inputs from the vertical sinkhole walls (Tamalavage et al., 2018).

Figure 1 Map of Blackwood Sinkhole within the context of the Bahamas region. (A) The western North Atlantic region noting the location of the Blackwood Sinkhole on the Little Bahama Bank. (B) Islands on the Little Bahama Bank, and the position of Blackwood Sinkhole on Great Abaco Island. (C) Aerial photograph of Blackwood Sinkhole and the relation (220 m) to the Atlantic Ocean. Red dots indicate coring locations, including BLWD-C7, used in this study. Indicating the location of the core sampled in this study, BLWD-C7. (D) Radiograph and picture of core BLWD-C7. Figure adapted from (Tamalavage et al., 2018).

Based on a detailed analysis on the changes in organic matter (OM) provenance through time using stable carbon isotopes (δ¹³C_org) and the C:N ratio, Tamalavage et al. (2018) documented how the primary source of OM accumulating in Blackwood Sinkhole changed through time. While OM was always typically derived from primary productivity and erosion of terrigenous materials, the dominant source was not always constant. For example, OM was primarily derived from the terrestrial landscape from ~3000 to ~1500 cal yrs BP. At ~1500 years ago, however, OM deposition was dominated by constituents of authigenic OM produced in situ. Finally, during the last 1000 years, organic matter input into Blackwood Sinkhole was dominated by inputs from an adjacent wetland that colonized the epikarst surface at around the same time. These wetlands colonized the favorable habitat that was created when concomitant regional sea-level and groundwater-level rise flooded topographic lows adjacent to Blackwood Sinkhole (Tamalavage et al., 2018; Fall et al., 2021).

Material and Methods

Sample Collection and Handling

A 7.6 cm diameter push core (BLWD-C7; 26.79°N, 77.42°W, depth: 34 ± 1 m below sea level) was collected on 29 July 2018 with a polyvinyl chloride pipe, using advanced technical scuba diving procedures that followed guidelines established by the American Academy of Underwater Sciences. The periphery of the sinkhole bottom was targeted for core sampling in an attempt to re-collect the last 3000 years of sedimentary deposition (van Hengstum et al., 2016). In the lab, the 89-cm long core was sectioned lengthwise into a working and archive half (subsequently stored at 4°C), photographed, radiographed, and the lithology was qualitatively described (Schnurrenberger et al., 2003). The working half was split into 24 stratigraphic layers, based on sediment color and texture, and stored at –80°C for further analyses (Supplementary Table 1).

To validate the age of the core and compare to previous results, five terrestrial plant macrofossils (e.g., leaves, twigs) were sent for radiocarbon dating of at the National Ocean Sciences Accelerator Mass Spectrometry facility at Woods Hole Oceanographic Institution (Woods Hole, MA, USA, 20.1−22.2 cm, 32.9−35 cm, 52.1−54.2 cm, 64.5−67.7 cm, and 79.5−86.9 cm). Conventional radiocarbon ages were calibrated into Calendar years Before Present (cal yrs BP, present is considered 1950 Common Era) with IntCAL13 (Reimer et al., 2013). A final downcore Bayesian age model for BLWD-C7 was computed using the R program Bacon v2.2 (Blaauw and Christen, 2011) to provide probability estimates at each core depth (Supplementary Figure 1) using three of the five samples (32.9−35 cm, 52.1−54.2 cm, 64.5−67.7 cm), and the surface designated as −68 Cal yrs BP to mark the date of collection in the year 2018. Two samples dates were rejected from the final radiocarbon age model (Figure 2B; Supplementary Figure 1). The shallowest radiocarbon result from the base of the prominent sapropel horizon at 20.1 to 22.2 cm in BLWD-C7 was not used in the age model because it provided a calibrated age result that was too old relative to previous higher-resolution dating results (Supplementary Table 2). Our previous, more detailed, stratigraphic dating results indicate a much younger age for this sapropel horizon (291-223 Cal yrs BP, Tamalavage et al., 2018), so the dated terrestrial macrofossil must have been an older plant fragment that resided on the terrestrial surface prior to its erosion into Blackwood Sinkhole. Similarly, the deepest radiocarbon calibration result from 79.5−86.9 cm in BLWD-C7 was rejected because it caused significant change in sedimentation rate that was inconsistent with previous results (Supplementary Figure 1), and as such is also suspected to be an old, reworked terrestrial plant fragment.

Figure 2 (A) Radiograph and picture of core BLWD-C7. The white star indicates the sapropel layer that is synchronous across all cores from Blackwood Sinkhole. (B) Downcore age model for BLWD-C7 showing core to be up to ~1900 years old. Radiocarbon results were calibrated into calendar years before present (present being 1950) with IntCAL13 and visualized using Bacon v2.2. The calibrated ¹⁴C dates are shown in the transparent blue regions, the grey coloring indicates more likely calendar ages, and the red line shows the best model based on the weighted mean age for each depth. (C) C:N ratio based on the percentage of (D) total nitrogen (%) and (E) total organic carbon (%), in accordance to depth in centimeters.

Measurement of Carbon and Nitrogen Content

Total carbon and nitrogen content was measured on 24 subsamples downcore from BLWD-C7 (Table 1) that were also targeted for 16S rRNA gene analysis. First, the subsamples were freeze dried overnight, homogenized, and 2 to 6 mg of ground sample were placed into tin capsules and measured on a Costech instruments ECS 4010 CHNSO Analyzer (Costech Analytical Technologies) to measure total carbon (TC) and total nitrogen (TN). Data calibration was determined relative to acetanilide and standard reference material for organics in marine sediment according to the National Institute of Standards and Technology (NIST). To measure organic carbon (with mass correction applied), the samples were acidified using 8 mL of 1M HCl for 24 h or until effervescence ceased, then desiccated at 60°C, and re-homogenized. The ground, acidified samples were weighed (0.5 to 1.8 mg) into silver capsules then processed on the CHNSO analyzer (Tamalavage et al., 2018). The potential loss of carbon from the direct acidification process was corrected by multiplying the percent of sample remaining (post-acidification weight subtracted from pre-acidification weight/pre acidification weight) (Tamalavage et al., 2018). The atomic C:N ratio was determined using the organic carbon (OC) acidified values divided by the unacidified Total Nitrogen (TN) values and multiplied by the molecular weight ratio (14.01/12.01) (Tamalavage et al., 2018). Precision on replicates measurements of C and N was within ± 2% weight percent.

Table 1 Downcore organic matter analysis [OC% (organic carbon corrected to the initial volume of sediment), TN% (total nitrogen) and C:N atomic ratios] as well as sequence-based calculations [Number of sequences (N_seq), sequence coverage, observed richness (S_OBS), and inverse Simpson index] for core BLWD C7.

16S rRNA Gene Sequencing and Analysis

Total DNA was extracted from separate sediment subsamples (1–2g) from each of the 24 horizons using the DNeasy PowerSoil kit (Qiagen, USA) (Supplementary Table 1) according to the manufacturer’s protocol, except more material was used than the recommended 0.25 g, and then stored at –20°C until PCR amplification. As recommended for sediment and low biomass samples (Sheik et al., 2018), a negative control, or extraction blank, using 0.5 mL of 0.2 µm filtered sterile MilliQ water was also extracted to identify possible contamination from the ambient laboratory and kit reagents. The Earth Microbiome Project PCR primers 515F-Y (5′-GTGYCAGCMGCCGCGGTAA-3′) and 926R (5′-CCGYCAATTYMTTTRAGTTT-3′) were used to target the V4-V5 region of the 16S rRNA gene (Parada et al., 2016). Thermal cycling was performed under the following conditions: initial preheating for 3 min at 94°C; 35 cycles of denaturation at 94°C for 1 min, annealing at 50°C for 1 min, and extension at 72°C for 1 min 45 s; final extension at 72°C for 10 min. PCR reactions were completed in triplicates and the products were pooled and cleaned using the MinElute PCR Purification Kit (Qiagen, USA). Amplicons (100-200 ng) were sent to the Texas A&M AgriLife Bioinformatics and Genomics facility for library preparation with the Amplicon library preparation kit (Illumina) and sequencing with Illumina MiSeq with 250 bp paired-ends.

Analyses of the 16S rRNA gene amplicons were completed using the software mothur as presented in the MiSeq Standard Operating Procedure (SOP) tutorial (Kozich et al., 2013), which included reducing sequencing and PCR errors, processing the improved sequences, running an alignment using the reference SILVA alignment (v132), removing poorly aligned sequences and undesirables, and pre-clustering the sequences into amplicon sequence variants (ASVs). The remove.lineage command using all the lineages observed in the negative control at taxonomic level 6 (genus) was used to remove any potential contaminants from final analyses. For the non-metric multidimensional scaling (NMDS) analysis, the taxa (class level) abundances were square root transformed and the environmental variables (geochemical data) were z-transformed, which allowed the variables to be compared using a similar scale. The NMDS analysis was performed using the envfit function (permutations = 9999) available in the vegan R package (Oksanen et al., 2019) and plotted with ggplot2 (Wickham, 2016). The ANOVA tests were performed using the aov command to fit an analysis of variance model in R (RCore Team, 2017). Finally, stratigraphically-constrained (i.e., age constrained according to radiocarbon results) hierarchical clustering analysis was performed using the package ‘rioja’ in R, which can provide a semi-quantitative tool to assess intervals of time when the microbial community structure was similar (Juggins, 2017). The sequences are available in GenBank under the BioProject #PRJNA639820. Accession numbers for each sample are available in Supplementary Table 3.

Results

Sedimentation

Basal sediment in BLWD-C7 is carbonate gravel, above which is laminated algal sapropel interbedded with layers of calcium carbonate (Figure 2A). This general lithology is similar to previous results (van Hengstum et al., 2016; Tamalavage et al., 2018). Most notably, Blackwood Sinkhole has a prominent layer of non-laminated algal sapropel that is located near the top of the stratigraphic architecture (see photographs in Figure 2 of van Hengstum et al., 2016), which was identified in BLWD-C7 from 23–25 cm. Tamalavage et al. (2018) provided a more precise median age range for deposition of this non-laminated algal sapropel layer in an alternate core sample of 291-223 cal yrs BP (BLWD-C3). The 2σ uncertainty age envelope for this layer with the radiocarbon age results for BLWD-C7 (this work) is similar to previous age results. However, the new radiocarbon results indicate that the specific location where BLWD-C7 was collected only started accumulating sediment ~1900 years ago instead of ~3,000 yrs BP as observed in BLWD-C2 from Blackwood Sinkhole (Figure 2B). This result further underscores that bathymetric complexities created by the original cave collapse impacted subsequent sediment focusing and deposition (van Hengstum et al., 2016).

In core BLWD-C7 (this study), C and N contents show variation in accord with the primary mineralogical and organic matter variability in the stratigraphic layers (Table 1 and Figure 2). In general, sedimentary layers with high TOC content (%) also have a corresponding high C:N ratio (Table 1 and Figure 2). A high C:N ratio (>20) usually indicates that the source of carbon fluxed to the sediment was likely dominated by terrestrial organic carbon input, whereas a lower C:N ratio (<15) indicates a higher proportion of carbon from primary productivity [e.g., phytoplankton, (Stenuite et al., 2007)] (Lamb et al., 2006; Fuentes et al., 2016). The C:N ratio fluctuates in Blackwood Sinkhole, with the lowest value at horizon 20.1−22.2 cm (C:N ratio of 14.4) and the highest value at 72.5–76.9 cm (C:N ratio of 39.5) (Figure 2). The sapropel layer seen at approximately 23–25 cm contains 15% TOC and has a C:N ratio of 19, indicating that the organic matter has both terrestrial and primary productivity sources (Tamalavage et al., 2018).

Microbial Community Composition Analysis Using 16S rRNA Gene Sequencing

We sequenced on average 749,985 reads for each of the 24 samples, for a total of 20,669,626 raw reads. After trimming and deduplication, between 277,000-551,107 reads per sample, or a total of 13,907,092 reads remained (Supplementary Table 3). We ran a negative control, which returned 3,764 sequences, to identify possible extraction kit contamination (Supplementary Figure 2). Contaminants made up ~5% of the entire dataset, and they consisted of classes and species known to be common contaminants of DNA extraction and sequencing kits reagents (e.g., Ralstonia, unclassified Xanthobacteraceae, Afipia, Streptococcus, Sphingomonas, unclassified Enterobacteriaceae, Pseudomonas, Bradyrhizobium, Micrococcus, Cutibacterium, Paenibacillus, Mesorhizobium, Pelagibacterales Clade I, Marivivens, Lactococcus, unclassified Burkholderiaceae, unclassified Sphingomonadaceae, and unclassified Rhizobiaceae) (Salter et al., 2014; Glassing et al., 2016). Using a conservative approach, ASVs belonging to the same genus as the contaminants were removed. After removal of poorly aligned and chimeric sequences, and decontamination, we were left with a total of 4,055,715 ASVs for taxonomic assignment (Supplementary Table 3).

Main taxonomic classes observed in all samples are common in anoxic sediment, such as wetlands and marine sediment, and belong to Dehalococcoidia, Gammaproteobacteria, Bathyarchaeota, and Campylobacteria (Figure 3 and Supplementary Figure 3) (Wang et al., 2012; Biderre-Petit et al., 2016; Franco et al., 2017). In addition to Bathyarchaeota, a number of other archaeal phyla (e.g., including Euryarchaeota, Nanoarchaeota, and Crenarcheaota) were also present (Figure 3). When looking at a deeper taxonomic level, negligible changes were observed within the class Dehalococcoidia, except a slight increase in the class Sh765B-AG-111 was observed down core (Supplementary Figure 3).

Figure 3 Class level prokaryotic relative abundance for each layer of the sediment core. Groups 1, 2, and 3 are divided by the dashed lines. Taxa highlighted in color are the ones considered in the discussion.

We also compared the samples using a stratigraphically-constrained cluster analysis, which provides a semi-quantitative time-wise (i.e., core depth) division of changes in the microbial community structure (Figure 4). Interestingly, the output dendrogram can be assessed in two different ways. First, three groups emerge on the dendrogram at Bray-Curtis Dissimilarity Index of 2.3: (a) group 1 (0-32.8 cm) spanning 545 cal yrs BP until present, (b) group 2 (32.9-42.3 cm) spanning 545-930 cal yrs B, and (c) group 3 (43 to 86.9 cm) spanning 930 to 1908 cal yrs BP. For example, Gammaproteobacteria was mainly present in Group 1 and Group 3. Gammaproteobacteria were mainly dominated by the order Thiomicrospirales in Group 1 and displayed greater taxonomic variations in Group 3, especially in the deeper part of the core (Supplementary Figure 3). Campylobacteria was mainly present in Group 3, where it was dominated by the families Sulfurovaceae and Arcobacteraceae. In Group 1, there were more unclassified Campylobacterales members and Thiovulaceae, which were absent from Group 3. Other taxa present showed no major affiliation with any group. Secondly, the dendrogram can be perceived slightly broader, and at a Bray-Curtis Dissimilarity Index of 2.75 two groups emerge: (a) group 1 (0-42.3 cm) spanning the last 1000 years, and (b) group 2 (43-86.9 cm) contains samples than are older than 1000 years.

To determine if the beta diversity, or the distribution of the microbial community, was related to carbon content, the relative abundance of the different taxonomic classes was analyzed in relation to total organic carbon content (%), total nitrogen (%), and C:N ratio in an MSDS (Figure 5). The samples clustered based on their age group. While total organic carbon (p = 0.737) and total nitrogen (p = 0.258) are not statistically correlated to the relative abundance of the microbial taxa, the C:N ratio is weakly statistically correlated to the relative abundance for the microbial taxa (p = 0.096) (Figure 5). Calculated diversity indices include both taxon richness and evenness, both of which varied throughout the core. The shallowest samples, BLWD-C7 2–3 cm, 6–13.5 cm, and the three samples belonging to Group 2, 35.2–39.4 cm, 39.6–42.3 cm, and 43–45 cm, displayed the highest diversity indices (Inverse Simpson), with 126.4, 92.6, 495.6, 326.8, and 152.5, respectively (Table 1). There was no correlation between the inverse Simpson diversity index and C:N ratio values (ANOVA p = 0.897).

Figure 4 Stratigraphically-constrained hierarchical cluster analysis of microbial community downcore based on the three previously determined age groups based on source of organic matter. Blue bars, linked by a grey area line, represent the number of ASVs for each class.

Figure 5 Non-metric multidimensional scaling (NMDS) ordination of each sample based on their microbial composition. Groups previously identified are labeled as orange circles (Group 1), blue triangles (Group 2), and teal squares (Group 3). Each line is correlated to the ordination and represents the direction and strength of the carbon influence (C:N ratio, OC, and TN). C:N ratio is significantly correlated (p = 0.096).

Discussion

The most interesting result is the timewise division of the subsurface prokaryotic community structure that generally correlates to previous evidence of changes in OM source and deposition in Blackwood Sinkhole. According to the 16S rRNA gene data, the most prevalent classes within the entire sediment core are Dehalococcoidia, Bathyarchaeota, Campylobacteria, and Gammaproteobacteria. These classes exhibit meaningful timewise distributional changes in the ~1900 year old succession of BLWD-C7, suggesting that the type of organic matter at time of burial provided resources for the enrichment of these microbial communities.

At the broadest level, the dendrogram can be divided at Bray-Curtis Dissimilarity Index of 2.75 (Figure 4), which creates two groups: (a) group 1 (0-42.3 cm) spanning the last 1000 years, and (b) group 2 (43-86.9 cm) contains samples than are older than 1000 years. In this greater division of the dendrogram, the older Group 2 has a higher proportion of Campylobacteria, whereas the younger Group 1 is characterized by increased proportion of Gammaproetobacteria. This bipartite division of the dendrogram compares well with the distinct intervals of OM sourcing in Blackwood Sinkhole determined from a multicore analysis of δ¹³C, δ¹⁵N, and C:N (Tamalavage et al., 2018). The oldest sediment deposited in Blackwood Sinkhole from 2990-1520 cal yrs BP had OM that was primarily derived from terrestrial sources. From 1502-1009 cal yrs BP, OM derived from freshwater primary productivity produced in situ in the upper part of the stratified water column. At 1300 years ago, however, a wetland began colonizing the topographic lows adjacent to Blackwood Sinkhole (van Hengstum et al., 2016). Pollen grains preserved in the stratigraphy of Blackwood Sinkhole indicate that new vegetation on the landscape which increased around ~1000 years ago included increased abundance of mangrove species (e.g., Lagunclularia, Avicennia) and emergent macrophytes (e.g., Typha, (Fall et al., 2021)). As a result of wetland expansion, by ~1000 years ago, organic matter produced within the wetland eroded and transported into Blackwood Sinkhole, and became the dominant source of OM into the system (Tamalavage et al., 2018). It would appear that increased delivery of wetland-derived organic matter at ~1000 cal yrs BP generally favored expansion of subsurface populations of Grammaproteobacteria, at the expense of Camphylobacteria. In addition, it seems that the expanded primary productivity in the upper water column from 1500 to 1000 years ago, which would provide more bioavailable organic carbon, was not as meaningful for the eventual subsurface prokaryotic community structure.

In balance, however, a bipartite division of the dendrogram does not completely explain the changes observed in the subsurface microbial community structure. If the dendrogram is alternatively divided at a Bray-Curtis Dissimilarity Index of 2.3, three different communities of prokaryotes emerge (as depicted in Figure 4): Group 1 (0-32.8 cm) spanning 545 cal yrs BP until present, (b) group 2 (32.9-42.3 cm) spanning 545-930 cal yrs B, and (c) group 3 (43 to 86.9 cm) spanning 930 to 1908 cal yrs BP. The numerical separation is most likely caused by the salient decrease in Camphylobacteria in sediment that was deposited from 545 to 930 cal yrs BP.

The relative abundance of Dehaloccoidia and Bathyarchaeota showed little variation throughout the core. However, studies completed in anoxic sediment elsewhere have often discovered an abundance of the class Dehalococcoidia (phylum Chloroflexi) (Cetecioğlu et al., 2009; Biderre-Petit et al., 2016). Specifically, members of the class Dehalococcoidia are known to only inhabit environments with strict anoxic conditions, such as river sediments, aquifers, and sludge (Löffler et al., 2013). Based on metabolic potential found in metagenome assembled genomes (MAGs), some strains are known to rely only on halogenated compounds for anaerobic respiration, while others are anaerobic acetate fermenters (Hug et al., 2013). Members of the phylum Chloroflexi are commonly found in anoxic environments, including within the water column (Cetecioğlu et al., 2009). Analysis of Chloroflexi MAGs from aquifer sediments (Hug et al., 2013) and a bioreactor (Kindaichi et al., 2012) also revealed that some strains may also have the ability to degrade cellulose, which is the main component of plant cell walls. Chloroflexi members may be thriving on the terrestrial input of organic matter in Blackwood Sinkhole. Bathyarchaeota are widespread in anoxic sediment in both marine and freshwater (Zhou et al., 2018). Their wide distribution may be caused by the fact that these microbes are metabolic generalists able to anaerobically utilize detritus (proteins and polymeric), aromatic compounds, methane and methylated compounds, and a wide range of organic matter (Zhou et al., 2018). Bathyarchaeota are able to break down carbohydrates of terrestrial origin (Lazar et al., 2016) and lignin (Yu T. et al., 2018), which are important sources of organic matter in Blackwood Sinkhole. Metagenomic analyses of the microbial communities inhabiting Blackwood Sinkhole may shed more light on the potential impact of the members of the Dehalococcoidia and Bathyarchaeota in this anoxic system influenced by both aquatic and terrestrial OM.

Campylobacteria was generally more abundant in Group 3, coeval with deposition of OM primarily from the terrestrial surface and primary productivity. Campylobacteria are nitrogen fixers and sulfate reducers commonly found in sulfur-rich environments such as hydrothermal vents, or salt marshes (Drexler et al., 2020), since Campylobacteria has been observed with the root system of Spartina (McClung and Patriquin, 1980; Miller et al., 2007). While we do not have companion measurements of porewater sulfate concentrations, sulfate has been found in higher concentrations in anoxic subsurface environments, and increasing with depth in another stratified anchialine sinkhole in a riverbed (Bottrell et al., 1991; Davis and Garey, 2018), suggesting that the higher prevalence of Campylobacteria in Group 3 may be influenced by sulfate concentrations. The higher relative abundance of Thiovulaceae (specifically members the family Sulfurimonas, a widespread sulfur-oxidizing species found in sediments, hydrothermal vent fields, aquifers and subsurface environments) in the upper part of the core in Group 1 also supports the idea that the sulfur cycle is a determining factor in shaping the microbial communities (Supplementary Figure 3).

Gammaproteobacteria was abundant in Group 1 and towards the end of Group 3, which were both characterized by OM originating from wetlands and mangroves, or terrestrial input, respectively, which are sources of carbon that are less labile, therefore not as easily degraded by microorganisms than carbon of photodetritus origin (Guillemette et al., 2013). Gammaproteobacteria, known to be facilitating carbon fixation (Dyksma et al., 2016), have also been found in high abundance in coastal sediment from sites in Western Europe, Germany, and Australia (Lenk et al., 2011; Dyksma et al., 2016). In the layers where carbon isn’t as labile, therefore less available for catabolism, Gammaproteobacteria could be playing a key role in the system by fixing CO₂ to produce usable organic carbon compounds, but more studies are required to confirm.

Alphaproteobacteria were found only in the top horizon (2–3cm) and the bottom 20 cm (72.5–86.9 cm) of the core. These bacteria are known for a variety of metabolic abilities within anoxic sediment, including nitrogen fixation, ammonia oxidation, and methylotrophy (Acosta-González and Marqués, 2016; Williams et al., 2017). Specifically, the order Rhodobacteraceae within the Alphaproteobacteria class are associated with marine algae and seagrasses in coastal sediments (Raulf et al., 2015; Park et al., 2016). Alphaproteobacteria and other microbes, especially nitrate reducers, are commonly found in the water column. They are known to decrease with subsurface depth (Lloyd et al., 2020), which can explain their dominance only in the top most sediment. Our system is marine groundwater fed, as such, an influx of fresh nutrients at the base of the sediment could fuel bacteria such as Alphaproteobacteria. This groundwater input could explain why Alphaproteobacteria were seen at the bottom of the core as well as at its surface, but a detailed nutrient analysis would provide more information regarding Alphaproteobacteria distribution in the subsurface.

Conclusions

Overall, the entire core displayed variability in carbon source (terrestrial/aquatic derived), which promotes a diverse subsurface microbial community. The microbial communities were dominated by members of the Chloroflexi, Bathyarchaeota, Proteobacteria, and Campylobacteria phyla, all of which are commonly found in coastal anoxic sediment. The microbial communities varied based on landscape change and variability in in situ productivity, suggesting that the source of carbon at time of deposition is a driver of the microbial community’s presence within this anoxic environment at Blackwood sinkhole. Further exploration of microbial communities’ relationships to carbon in other sinkhole environments will expand our understanding of the role of microorganisms living in anoxic sediment within the carbon cycle.

Data Availability Statement

The data presented in the study are deposited in the NCBI repository, BioProject accession number PRJNA639820, available at: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA639820.

Author Contributions

JL and PvH designed the study. CR generated the data, performed the analyses and wrote the first draft of the manuscript. AT helped with total C and N measurements, and manuscript preparation. PvH collected the sample and participated in the manuscript preparation. JL edited the manuscript, with contributions from all. All authors contributed to the article and approved the submitted version.

Funding

Funding for this research was provided by the Center for Dark Energy Biosphere Investigations (Research Grant to JL and PvH), Texas Sea Grant (Graduate student Grant in Aid to CR), and the National Sciences Foundation (EAR-1703087 to PvH).

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

This is C-DEBI contribution # 593.

Supplementary Material

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

References