Ground-truth surveys of seagrass habitat
To obtain georeferenced field data on benthic cover levels from habitats of the Bahama Banks, we employed two similar, in-water survey and image approaches: (1) swimmer-based photo-transects; and (2) tow board photo transects (Supplementary Fig. 6), resulting in a total of 2542 surveys.
For (1), free-divers swam over the bottom of the seafloor at a fixed height with a digital camera (Canon 5D mIV, GoPro Hero) set to capture images manually. Photographs were captured using automatic settings in a 1.0 m × 1.0 m footprint, 1.5 m above the seafloor following [39]. A center console vessel was used to run the transects at distances of 5–7 km, whereby the free-diver would capture successive photos at a horizontal distance of between 400–800 m, and the location was logged using either a handheld GPS (Garmin GPS 73) or a boat-mounted GPS with a depth sounder (Garmin EchoMap DV). Transect locations were chosen based on a priori local expert knowledge of varying benthic cover in the region. Surveyed areas included: southern New Providence (24.948862°, −77.387834°), southeast of New Providence (24.980265°, −77.229168°), south of Rose Island (25.066268°, −77.160063°), the middle Great Bahama Bank (24.735355°, −77.212998°), and the northern Exumas (24.729973°, −76.889488°). For (2), snorkeling observers were pulled from a research vessel on tow boards affixed with underwater action cameras (GoPro Hero 3+) traveling at ~1 m/s. The start and end of a tow were delineated with either a handheld GPS (Garmin eTrex 30) or a boat mounted GPS with depth-finder (Garmin EchoMap DV), and tows proceeded in a straight line recorded by the GPS. Cameras recorded images at 0.5 Hz throughout the tow, starting in conjunction with creating a waypoint. Samples (i.e., paired image and geolocated point) were sub-selected from the tow once movement began, at the midpoint of a tow, and immediately before movement stopped. Images were manually quality controlled such that if a selected image contained obstructions or was out of focus, the nearest clear image was selected to replace it. If no images within 10 s were clear (i.e., 10 m maximum spatial error), the sample was discarded. If the GPS track contained gaps or segments larger than 10 m, only images/point pairs at the start and end waypoints were sampled.
Surveys focused on historical fishing grounds for queen conch (Lobatus gigas) between 2015 and 2018 following the sampling design and methods of ref. 32. A stratified random design was used to allocate 6000 m2 of observation effort into each cell of a 1’ by 1’ grid placed over each fishing ground. This effort was split into multiple tows between 200 and 1000 m in length, thus images were separated by at least 100 m.
Fishing grounds extended from the edge of a deepwater sound to between 7 and 10 km up the bank and were limited to the depths used by freediving fishers. Surveyed fishing grounds included: the Exumas (24.382207°, −76.631058°), the southwestern Berry Islands (25.455529°, −78.014214°), south of Bimini (25.375592°, −79.187609°), the Grassy Cays (23.666864°, −77.383547°), the Joulter Cays (25.321297°, −78.109251°) and the southeast tip of the Tongue of the Ocean (23.376417°, −76.621943°). For details on image processing, see section on remote sensing below.
Sediment coring
To gather the sediment cores analyzed for organic carbon content on the Bahama Banks, we collected samples from various benthic habitats that included varying densities of seagrass habitat (Thalassia testidinum and Syringodium filiforme). We percussed, via SCUBA, an acrylic cylinder tube perpendicular to the seafloor into marine sediment until rejection at various penetration depths up to 30 cm. The sample was then extracted vertically from the marine sediment and capped at the bottom to avoid loss of material. This sample was then transported vertically through the water column to a research vessel where it was removed from the coring device and immediately capped on top with an air-tight cap. Compression rates were negligible (~5 cm) across the first 5 cores, and as such were not subsequently measured. The samples were then labeled, photographed, geotagged, and the first 30 centimeters of each core was extruded. To complete the extrusion process, we placed each sample on top of a capped piston device in the same orientation as collection (deepest portion of collected sediment still on the bottom). The bottom cap was removed to thread the acrylic cylinder tube onto the piston device and then was lowered to various measured lengths to collect corresponding depth sections of the sediment core. These sections were sliced (every 1–5 centimeters), labeled, and placed into whirl pack bags to collect the wet weight of each sample. All samples were then frozen and stored for future laboratory analyses. All samples were dried in a laboratory oven at 55 °C for 48 h until constant dry weights were reached. The samples were then weighed to collect their corresponding dry weights. The dry bulk density (DBD) was calculated by diving the sample dry weight (g) by the sample volume (cm3). The samples were then further ground with a mortar and pestle until a homogeneous fine grain size was achieved. Sediment samples collected from the Exuma Cays (142 samples from 16 cores) were analyzed for Corg content. Sediment samples were weighed accurately into silver capsules and acidified with 4% HCl until no effervescence was detected in two consecutive cycles. The samples were then dried in a 60 °C oven overnight, encapsulated into tin capsules and analyzed using an Organic Elemental Analyzer Flash 2000 (Thermo Fisher Scientific, Massachusetts, USA). We then conducted a standard loss on ignition (LOI) methodology at our laboratory facility (Braintree, Massachusetts, USA) for all the samples. Each sample was subsequently sub sampled with 5–15 grams of representative material and placed into a ceramic crucible to collect its mass. The crucibles were then loaded into a separate muffle laboratory oven and heated at 550 °C for 6 h. Upon completion of this muffle, the crucibles were then immediately weighed to collect the LOI of organic material from each sample, defined as the weight lost in the muffle (g) divided by the subsample dry weight (g). A fitted regression between the Corg and LOI from the Exuma Cays cores was generated (Supplementary Fig. 7), and then used to predict the sediment Corg contents from LOI measurements in the Grand Bahama cores. Sediment Corg stocks were quantified by multiplying Corg and DBD data by soil depth increment (1–5 cm) of the sampled soil cores. The cores from the Exuma Cays (15 cm) and Grand Bahama (30 cm) were collected with different depths, we therefore fitted a regression between Corg stock in 15 cm-depth and Corg stock in 30 cm-depth for the Grand Bahama cores (Supplementary Fig. 8) and used this regression to extrapolate Corg stock of the Exuma Cays cores into 30 cm-depth. Moreover, to allow direct comparison among other studies27, the Corg stock per unit area was standardized to 1 m-thick deposits by multiplying 100/30.
Tiger shark tagging
The research and protocols conducted in this study complies with relevant ethical regulations as approved by the Carleton University Animal Care Committee. The shark data used in this paper were collected as part of a multi-year, long-term research program evaluating the interannual behavior and physiology of large sharks throughout the coastal waters of The Commonwealth of The Bahamas23. All sharks were captured using standardized circle-hook drumlines33 on the Great and Little Bahama Banks throughout the country, focusing efforts in three primary locations: off New Providence Island, the Exuma Cays, and off West End, Grand Bahama, from 2011–2019. All sharks were secured alongside center console research vessels and local dive boats, where their sex, morphometric measurements, and blood samples were taken. A mark-recapture identification tag was applied to the shark at the base of the dorsal fin. Some of the sharks sampled in the present study were also tagged with a coded acoustic transmitter which was surgically implanted ventrally into the peritoneal cavity and then sutured, as part of a concurrent study on shark habitat use and residency within the region23.
Pop-off archival satellite tags were affixed to eight tiger sharks (seven female, one male; 298 ± 28 cm total length; mean ± SD) in The Bahamas from 2011–2019, permitting measurements of swimming depth and water temperature recorded at either 4-min (Sea-Tag MODS, Desert Star Systems LCC, USA) or 10-s intervals (miniPAT tags, Wildlife Computers, USA). Pop-off satellite tags were inserted into the dorsal musculature of the sharks using stainless steel anchors and tethers. All pop-off satellite tags were either recovered manually, permitting access to the full time-series, or popped-off and transmitted their data to an Earth-orbiting Argos satellite, resulting in a subset of the full time-series (transmission frequencies: 2.5 min [miniPAT], 10 min [PSATGEO], daily average [Sea-Tag MOD]). Tiger shark positions were estimated from the satellite data using tag-specific proprietary state space algorithms from Wildlife Computers (GPE3; based on ref. 34) and Desert Star Systems35. With miniPAT tags, positions were further filtered to remove the least reliable positions (<0.1 observation score). Tracking durations with reliable positioning estimates were variable (mean = 144 days; 44 to 376 day range). Descriptive statistics of depths experienced by sharks (n = 176,206) were generated and depth use patterns were plotted for the periods where reliable positioning data were available (Supplementary Figs. 4 and 5). Shark satellite positions were filtered to include just the positions located on the Bahamas Banks. Kernal utilization distributions (KUD) were calculated using the adehabitatHR package36, from which 95 and 50% KUD polygons were extracted. This was conducted with satellite tracking data from 5 sharks for which there were sufficient datapoints for this analysis.
Camera tag biologger packages were affixed to a subset of 7 mature tiger sharks on the Great and Little Bahama Banks from 2016–2020, using two methods: (1) capture and release using hook and line, following the same methods as above; and (2) in-water placement on free-swimming tiger sharks. In (1), we built custom camera-tag packages using a positively buoyant material (Diab Syntactic © non-compressible foam). All cameras used were forward-facing, and uniidirectional, with the exception of one unit which was a 360-degree camera24. Two asset recovery tags were secured to the center of the biologger payload using clear silicone: a satellite tag (SPOT-386A, Wildlife Computers, Redmond, WA, USA) and a VHF radio tag (F1840B, Advanced Telemetry Systems, Isanti, MN, USA). Two stainless steel nuts were added to the package to provide forward-facing ballast to reduce the buoyancy from the camera housing, thus allowing the tag to float on the surface in a manner which maintained the vertical orientation of the satellite and radio tag antennae in air. The entire biologger package was attached to the left side of the shark’s dorsal fin by drilling two small holes and threading two connected, biodegradable cable ties through and around the package. The heads of the cable ties were then joined together via the eyes of a dissolvable galvanic timed-release swivel (A2 model, Neptune Marine Products, Port Townsend, WA, USA), which would eventually corrode in seawater after an estimated period of ~24 h (swivel was pre-dissolved to permit a short-term deployment), thereby allowing the positive buoyancy of the package to cause it to naturally release and come off the animal. Once attached, the camera was activated for recording and the shark was released. The entire time to collect all animal data, apply tags and attach the biologger was 12 min. In (2), small action cameras (GoPro Session) were inserted into custom float packages attached to stainless-steel clamps and were placed firmly on the dorsal fins of free-swimming tiger sharks. Dissolvable galvanic swivels were used as above, programmed to corrode after ~6–12 h. A single VHF transmitter tag as above was included in the package to aid in asset recovery.
Collectively, camera tag deployments ranged from 55.43–117.06 min of recorded footage, with a mean of 80.81 ± 11.65 min (Additional Supplementary Information). All camera tags were ultimately recovered by traveling to the most recent ping from the satellite tag while actively using a VHF radio to locate the pings from the VHF transmitter. Swimming speed estimates were obtained from the camera tag-packages by visual inspection of the footage and conversion from tailbeat frequency (TBF). Footage was played back at 1.5x speed, with the sharks’ head in the lower center of the screen. To estimate TBF, an observer manually tallied the number of times the shark moved its head from the right-hand side of the screen to the left-hand side and back to the right-hand side, counting this as one tailbeat. The total number of tailbeats per minute was recorded. Speed (S; m s−1) was estimated as S = SL*TL*TBF, where SL is average stride length, 0.36 body lengths/tail beat, which was digitized and interpolated from a previous study on tiger sharks37, and TL is total length in m. This was also used to calculate the total distance traveled by each shark. Classification of time spent over seagrass habitat was determined by counting the cumulative duration (minutes) when each tiger shark was seen swimming over any type of seagrass habitat, regardless of species or density. Total percentage of time over seagrass was then calculated for each shark by dividing time by the total track length (minutes), and the geolocation where each animal was tagged and released was used for spatial reference.
Empirical remote sensing
Remote sensing was used to map benthic seafloor habitat extent across the large spatial scales as seen in the study area38,39,40. Seven Landsat 8 (OLI) tiles (path/row, see Supplementary Table 3) covering study sites around The Bahamas were used in this study. Landsat 8 was chosen as it has among the best potential to facilitate large-scale seagrass density mapping due to its strong temporal resolution and algorithms for water-based corrections. Landsat has equivalent spectral characteristics to many other sensors (e.g., Sentinel-2, Hyperion, SPOT, ASTER, CBERS). We recognize Landsat has slightly lower spatial resolution than Sentinel-2. Images from 2019 and 2020 which had minimal to no cloud over near shore areas, which fell within the depth range of seagrass species, were chosen from web server of the United States Geological Survey (USGS; https://earthexplorer.usgs.gov/). More tiles would have been obtained; however, quality cloud-free images were not available for certain areas, such as the southern Bahamas. It was assumed that there were no significant changes between image acquisition years (less than two; Supplementary Table 3) and file data collection periods. Earth Resources Observation and Science ortho-rectified and terrain corrected (L1T) all Landsat imagery (https://www.usgs.gov/core-science-systems/nli/landsat/landsat-levels-processing). The visible (red, green and blue) and near-infrared bands, with 30 m spatial resolution of OLI data, were involved in cloud detection and land mask, while only water penetrating visible bands were used for seagrass mapping. Prior to seagrass mapping, the Fmask cloud detection algorithm41 was used for masking clouds, shadows, and land cover from each Landsat tile. Next, all raw OLI images (visible bands) were converted into top-of-water reflectance for radiometric correction following standard procedure suggested by the USGS (https://www.usgs.gov/core-science-systems/nli/landsat/using-usgs-landsat-level-1-data-product). For this study, the simple Dark Object Substract 1 (DOS1) atmospheric correction method plugin42,43 was implemented in the QGIS (v. 3.18). A water column correction44 was applied to radiometric and atmospheric corrected imagery. The depth invariant bottom index retrieved from the bi-plot of the reflectance of visible band-pairs15, for blue-red, red-green and green-blue, were layer stacked. Finally, all tiles were joined in a mosaic and cropped a subset of 170,388 km2 prior to seagrass classification. Excluding cloud, shadows, and land, about 88,000 km2 of ocean area were used for seagrass cover mapping.
A total of 2542 georeferenced, in situ field photos were collected uniformly along the bank of the Bahamas. The dominant benthic vegetation species belonging to the field photo were identified by expert knowledge. Seagrass cover percent for each photo was precisely estimated using image thresholding technique available in ImageJ (v. 1.53e)45, a java-based image processing software (https://imagej.nih.gov/ij/). All color (RGB) photos were converted into 8-bit grayscale prior to thresholding. Sand cover can easily be determined through automatic thresholding technique, from where fraction seagrass and non-seagrass cover can easily be determined through substruction of sand cover from total area of each photo. All field data thus were divided into four seagrass cover-types (I represents <25%, II represents <50%, III represents <75%, and IV represents <100% seagrass cover), non-seagrass and submerged sand substrates (Supplementary Table 4). Half of field data were used for training OLI data and half for accuracy assessment. Whereas the 1 m × 1 m ground-truthing images were much smaller than the 30 m × 30 m pixel size of the remote sensing product, we found that adjacent pixels tended to be homogeneous as spatial variability in seagrass configuration is small along the Bahamas Banks where environmental gradients shaping seagrass ecosystems are smooth (Supplementary Fig. 6).
We then used machine learning neural network (NN) algorithms46,47 for seagrass classification, which was carried out in ENVI (v. 5.3). For the NN algorithm, the key parameters optimum for classification of seagrass cover classes were: training threshold contribution (0.8), training rate (0.1), momentum (0.9), and number of hidden layers (1 for non-linear). A confusion matrix was calculated using ground truthed ROIs47 to assess the accuracy of NN parameters. The accuracy measures48 were expressed in terms of overall accuracy, producer’s and user’s accuracies, and kappa coefficient. The accuracy of classified Landsat 8 images produced by using NN analysis showed evidence for misclassification between individual seagrass density classes, thus yielding a relatively acceptable (70.2%) overall accuracy; however, “seagrass” was correctly mapped (Supplementary Table 5). Consequently, the Kappa value was also found to be low (0.6, i.e., <0.8), which was expected for mapping a large region like The Bahamas, at the spatial resolution of Landsat (30 m). Lower user’s (40%) and producer’s (56%) accuracies were achieved for the discrimination of sparse “Sg-I” class (<25% cover). Higher accuracies were recorded for medium density “Sg-II and III” classes, compared to the highest density “Sg-IV” class. The NN classifier can be considered effective since a lower number of predicted seagrass class pixels were assigned as either non-seagrass or sand pixels (Supplementary Table 6), resulting in an acceptable seagrass cover class map without over- or under-estimating seagrass cover areas for The Bahamas.
Comparison of remote sensing products
The current mapping efforts were integrated with several previous seagrass estimates18,19,20 to generate an integrated, composite range of estimates of seagrass coverage on The Bahamas Banks. This approach did two key things: (1) compared overlaps from each remote sensing product, to generate estimates of agreement on seagrass habitat; and (2) generated composite maps and calculated the spatial extent predicted from each category. This resulted in a range of estimates, from high to low probability. Seagrass shapefiles from ref. 20 were downloaded from the Allen Coral Atlas website (https://allencoralatlas.org/atlas/#5.39/24.3807/-76.0918). This product reported good accuracy for Caribbean seagrass (up to 67% match between observed and predicted), although it used ground-truthed sample images from The Dominican Republic and US Virgin Islands for their Bahamas predictions, which is not authentically representative. From empirical products18,19, seagrass estimates were derived from images (Fig. 1). Images were georeferenced using raster georeferencer in QGIS (Version 3.12.3). The remainder of analyses were conducted with R49 via RStudio50. Color bands were used to extract the seagrass estimates from each image ([19]—greyscale values 80–120 were assigned to seagrass; [19]—green values ≤ 170 dense seagrass [>70%], red bands ≤ 145 sparse seagrass [<70%]). Bands were determined by visual inspections for consistency with existing images.
To generate an overall weighted estimate of seagrass distribution, estimates from current mapping18,19 were integrated into a standardized raster grid. The grid was generated across the entirety of the banks with a resolution of 0.01° lat/lon, producing cell sizes of 1.23 km2. The resulting products had similar estimated seagrass extents to those reported in empirical research ([18] = 65,436 km2, estimated = 63,841 km2; [19] = 37,000 km2, estimated = 33,952 km2). Estimates from the three sources were used to generate an ensemble seagrass distribution estimate via a weighting scheme, where dense seagrass from current estimates (>50% seagrass) and dense seagrass from [18] (>70%) were given two votes, and low density estimates from these sources, as well as seagrass (density unspecified) from [19] were assigned one vote (Supplementary Table 2). Estimates from each source spanned varied regions of the banks; therefore, for each raster cell, the percent votes from the total available were calculated (Supplementary Table 2). Various cutoff percentage points were used to estimate the potential seagrass coverage on the banks based on the weight of evidence.
Post hoc analysis of ensemble model accuracy
Classifying linear habitat tows: Image streams of the benthos were collected by action cameras mounted beneath tow boards operated by skin diving observers in search of queen conch in The Bahamas. Images were taken at a frequency of 2 Hz. Boards were dragged at between 0.5 and 1.5 m/s depending on sea surface conditions and at a relatively constant speed chosen by the observer at the start of each tow. The start and end locations of each tow were marked via GPS with a resolution of ~2 m. The continuous nature of the image stream which included shifts in visibility and depth and caused a constantly fluctuating view of the benthos made categorical assignments qualitative. The dominant habitat type and approximate coverages were estimated from an initial view of the benthos, and then major habitat transitions were noted in subsequent images. Here we define a transition as a notable shift in habitat type (i.e., seagrass to sand) or coverage (i.e., 25 to 75% seagrass coverage) that persisted for more than 10 images (~5 m of bottom). At each transition, the sequence of the image was noted along with the new habitat type(s) and coverage(s). We used the same seagrass categories as those used in our remote sensing satellite imagery model: SG-I (>0 and <25%), SG-II (<50%), SG-III (<75%) and SG-IV (≥75%). We further divided the NSg category into C-I (<25% hard substrate with living coral), C-II (>25% hard substrate with living coral), MA-I (<25% macroalgae), MA-II (>25% macroalgae), G-1 (<25% gorgonian plain), G-II (>25% gorgonian plain) and included a hardbottom category for pavement. Only living benthic coverage was classified (i.e., random seagrass blades that were not rooted were not counted toward coverage). The same observer (A. Kough) classified all the habitat, had been a participant in collecting tow data at all sites and was familiar with how the benthos appeared both on camera and in person.
Comparing linear ground-truth data to integrated, ensemble model: Tows proceeded linearly at an approximately constant velocity thus position could be inferred from image order in between two geolocated points. For example, on a tow with 1200 images, image 300 corresponded with a location 25% between the start point and the end point. Points for each transition were calculated via image order and lines containing the habitat classes were created as shapefiles in ArcGIS. The rasterized ensemble estimate, which integrated four remote sensing estimates (including our own empirical estimate), was transformed into a polygon layer and then a georeferenced intersect was calculated between the Ensemble polygons and the habitat tow lines in ArcGIS. Benthic classes from the tow data were reassigned into seagrass (SG) or not seagrass (NSG) to correspond with the Ensemble’s predictions.
We chose a random subset of 197 tows covering 69 km from available data in the Berry Islands, Grassy Cays, Exuma Cays, Bimini and Jolter Cays to classify for ground-truthing. This resulted in an analysis of 126,969 additional benthic images, which fell into 97 cells within the ensemble model, with a mean of 2.5 tows/cell. Ultimately, we obtained 909 segments of the benthos along the tows and split by the Ensemble cells. Each cell was assigned the class of the Ensemble prediction, as our goal was to verify the voted classification. If part of a ground-truthing tow fell within a cell and contained any of the habitat type that matched the Ensemble prediction, it was considered a plausible verification of the classification. In addition, we calculated the ratio of the amount of towed distance that matched the Ensemble’s prediction against the total towed distance within each cell as a measure of agreement. Results suggested that the Ensemble’s predictions were verified by ground-truthing in most cells (~80%, Supplementary Fig. 3). For example, 83 of the Ensemble cells containing tow data for ground-truthing were classified as SG. Tow data verified that 65 of these cells contained areas dominated by seagrass. Further, in 50 of these cells most habitat encountered in the tows was seagrass.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
