Dataset: 16S rRNA sequences from subsurface microbial communities (Blackwood Sinkhole)

Methods are detailed in Risley et al., submitted. Briefly, a core was recovered from Blackwood Sinkhole, Bahamas, and divided into 24 stratigraphic layers. For each layer, DNA was extracted and the 16S rRNA gene was amplified and sequenced to characterize the microbial communities in relation to the C:N ratio.

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) (van Hengstum et al., 2016). Blackwood Sinkhole is 32 m in diameter, and this groundwater-fed basin has a stratified water column. The sinkhole hydrography is characterized by an upper ~10 m of meteoric water (1.4 psu), which transitions to a mixing zone below (10-15 meters below sea level), and this rests above anoxic saline groundwater (~15-40 meters 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 mixing from either invertebrate bioturbation, or physical mixing from wave or wind action. Basal sedimentary deposits in Blackwood Sinkhole are carbonate (karst) gravels, which abruptly transition to laminated sapropel (i.e., high organic matter) interbedded with carbonate horizons that were deposited over the last 3,000 years (van Hengstum et al., 2016). The laminated sapropel has a high total organic carbon (TOC) content ~10% (Tamalavage et al., 2018). However, sedimentation rate and location of sedimentation in the sinkhole has not been laterally or temporally uniform over time. This is because (a) the bottom geometry of the sinkhole basin is complex as inherited from the collapse of an original limestone cave, and (b) the sources and delivery mechanisms of sediment during the last 3000 years are linked to primary productivity, erosion of the vertical sinkhole walls, and deposition of organic matter derived from adjacent wetland habitats (Tamalavage et al., 2018). Using detailed radiocarbon dating on other core samples, constant sedimentation (~0.4 mm yr⁻¹) was initiated 3,000 Calibrated years Before Present (Cal yrs BP) on a core collected from the sinkhole periphery, with evidence of coarse-grained particle deposition during intense hurricane strikes on Abaco Island (van Hengstum et al., 2016). In contrast, sedimentation in the center of the sinkhole was delayed until ~1,600 Cal yrs BP, and sedimentation rates in the center were higher than the periphery (1.2 mm yr⁻¹).

Based on a detailed analysis on the changes in organic matter provenance through time using a 3-endmember mixing model using stable carbon isotopes (δ¹³Corg) and the C:N ratio, there were three periods of time in Blackwood Sinkhole when organic matter was dominated by a different source (Tamalavage et al., 2018). The oldest part of the record (Group 3, 1520-2990 Cal yrs BP) was characterized by primarily terrestrial organic matter deposition, and the middle part of the record was dominated by inputs from primary productivity (Group 2, 1009-1502 Cal yrs BP). Finally, the last millennium was dominated by organic matter inputs from an adjacent wetland (Group 1, 0-1008 Cal yrs BP). These wetlands were emplaced on the epikarst surface in response to concomitant regional sea-level and groundwater-level rise inundating topographic lows on the landscape adjacent to Blackwood Sinkhole, which created favorable wetland habitat (Tamalavage et al., 2018).

Sample collection and handling
A 7.6 cm push core (BLWD-C7; 26.79°N, 77.42°W) was collected on July 29, 2018 with a polyvinyl chloride pipe, using advanced technical scuba diving procedures following 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 3,000 years of sedimentary deposition (van Hengstum et al., 2016). In the lab, the core was sectioned lengthwise into a working and archive half (subsequently stored at 4℃), 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℃ for further analyses.

To validate the age of the core and compare to previous results, five samples (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) of terrestrial plant macrofossils (e.g., leaves, twigs) were selected for radiocarbon dating of at the National Ocean Sciences Accelerator Mass Spectrometry facility at Woods Hole Oceanographic Institution (Woods Hole, MA, USA). 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 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 (2018 year of collection). The shallowest radiocarbon result from the base of the prominent sapropel horizon at 20.1 to 22.2 cm in BLWD-C7 provided a calibrated age result that exceeded 500 Cal yrs BP, yet previous results indicate a much younger age for this horizon (291-223 Cal yrs BP, Tamalavage et al., 2018). Similarly, the deepest radiocarbon calibration result from 79.5−86.9 cm in BLWD-C7 caused significant change in sedimentation rate, and as such is suspected to be old and reworked terrestrial plant remains. These dates were rejected from the final radiocarbon age model.

Measurement of carbon and nitrogen content
Total carbon and nitrogen content measurements were performed on subsamples from the 24 sediment horizons from BLWD-C7. 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 CHN analyzer. 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 Total Nitrogen (TN) values unacidified 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.

16S rRNA gene sequencing and analysis
Total DNA was extracted from separate sediment subsamples (1-2 g) from each of the 24 horizons using the DNeasy PowerSoil kit (Qiagen, USA) and stored at -20˚C until PCR amplification. A negative control using 0.5 mL of 0.2 µm sterile MilliQ water was also extracted to identify possible contamination from the ambient lab and kit reagents. The PCR primers 515F (5'-GTGYCAGCMGCCGCGGTAA-3') and 806R (5'-CCGYCAATTYMTTTRAGTTT-3') were used to target the V4 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 was 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. Operational taxonomic unit (OTU)-based analyses using sequences clustered with the split method argument were also performed (i.e., rarefaction curves and heatmap of shared OTUs). Using a Bray-Curtis dissimilarity matrix, the beta-diversity between samples was examined and ordinated by non-metric multidimensional scaling (NMDS) in R (RCore Team, 2013), with overlaying the carbon content parameters (TC, TN, and C:N ratio) applying the 'envfit' function from the 'vegan' package (Dixon, 2003). Analysis of stratigraphically-constrained (i.e., age constrained according to radiocarbon results) hierarchical clustering using the package 'rioja' in R (Juggins, 2017). The sequences are available in GenBank under the BioProject #PRJNA639820.