Ontogenetic Patterns of Elemental Tracers in the Vertebrae Cartilage of Coastal and Oceanic Sharks

Introduction

Marine predators, including coastal and oceanic sharks, play critical roles in shaping ecosystem structure and function, primarily through direct (predation) and indirect (risk/behavior) effects on prey populations (Heithaus et al., 2008). Most shark species demonstrate longer lifespans, slower reproductive development, and lower fecundity than most exploited teleost fishes (Hoenig and Gruber, 1990). These traits enhance their susceptibility to population declines due to stressors such as habitat loss, climate change, and fishing pressure (Worm et al., 2013). The decline and loss of sharks in coastal waters can have profound effects on local ecological dynamics (Myers et al., 2007; Ferretti et al., 2010), so it is essential to understand the unique life history characteristics of these species to build effective management and conservation strategies. Central to this effort will be increasing our knowledge of species- and population-specific movement and habitat use patterns, migration routes, degrees of site fidelity, and population connectivity (Speed et al., 2011; Chapman et al., 2015).

Most estuarine and coastal environments are highly productive and diverse, and many shark species use these regions for reproduction, feeding, and as juvenile nursery grounds (Knip et al., 2010; Bethea et al., 2015). The reliance on nearshore habitats by a variety of shark species means that interspecific co-occurrence is likely, along with potential ecological interactions (Matich et al., 2017a; Heupel et al., 2019). The implications of multiple-predator co-occurrence can be widespread, including altered predation pressure for shared prey populations (Sih et al., 1998). However, some degree of differentiation in life history among sympatric predator species is expected, since resource partitioning often occurs to diminish potential competitive interactions (Papastamatiou et al., 2006). Therefore, examining the life history patterns of multiple shark species simultaneously can provide a far more informative view of ecosystem function than single-species studies alone.

Since large predators are often highly mobile and long-lived, tracking their movements and habitat use over extended periods of time and space can be prohibitively challenging (e.g., Block et al., 2011). One method for examining shark life history that has been developed in recent years is the analysis of natural biogeochemical tracers in vertebrae. Elements that are present in the environment are incorporated into the vertebral hydroxyapatite matrix during the biomineralization process, as the vertebrae grow concentrically over time. Although studies validating elemental uptake pathways into elasmobranch vertebrae are currently limited, concentrations of certain elements are thought to be associated with environmental variables such as salinity (Sr, Ba; Tillett et al., 2011) and temperature (Ba, Mg; Smith et al., 2013), dietary intake (Zn, Mn; Mathews and Fisher, 2009), maternal loading (Zn; Raoult et al., 2018), and unresolved physiological controls (McMillan et al., 2017). Consequently, chemical analysis of vertebrae sampled continuously along their radial growth axis may proxy how environmental conditions changed throughout an individual’s life (Scharer et al., 2012; Smith et al., 2013). This method has been used on sharks to retrospectively determine natal origin (Lewis et al., 2016; Smith et al., 2016; TinHan et al., 2020), movements across salinity gradients that correspond with ontogeny and reproduction (Tillett et al., 2011), and lifetime interactions with oceanographic features (Mohan et al., 2018).

This study focused on an assemblage of coastal and oceanic shark species in the northwestern Gulf of Mexico (nwGOM), from Galveston, TX, United States to the Mississippi River Delta, LA, United States (Figure 1). The biological diversity and habitat heterogeneity created by estuarine-oceanic gradients in the nwGOM makes it a particularly interesting area to examine the life histories of sharks. The selected species inhabit a gradient of habitats from freshwater-influenced estuaries to the open ocean, and experience varying degrees of potential co-occurrence and overlap. In approximate order of most estuarine to most oceanic, the species of interest in this study are Bull Sharks (Carcharhinus leucas), Bonnethead Sharks (Sphyrna tiburo), Blacktip Sharks (Carcharhinus limbatus), Spinner Sharks (Carcharhinus brevipinna), and Shortfin Mako Sharks (Isurus oxyrinchus). Using elemental tracers in vertebrae, the objectives of this study are to identify patterns in the elemental signatures of vertebral cartilage of five shark species across ontogeny, and to use the observed patterns to infer aspects of each species’ life history, such as movements and habitat use patterns.

Figure 1. Study area in the northwestern Gulf of Mexico (area within dashed square). The region where all sharks were collected is within the solid black square. Black contour lines begin at the continental shelf break (200 m depth).

Materials and Methods

Study Area and Species

Sharks used for this study were collected from the nwGOM, specifically near Galveston, Texas and the Mississippi River Delta, Louisiana (Figure 1). This region hosts a wide range of environmental conditions due to coastal and oceanic currents and freshwater input from estuaries. Bathymetry varies across the nwGOM; the continental shelf (to 200 m depth) extends far offshore (up to ∼100 nautical miles) in the western part of this region, but the shelf break is closer to shore (∼10 nautical miles) near the Mississippi River Delta (Bryant et al., 1990). The largest ocean circulation system in the GOM, the Loop Current and associated eddies, can extend into the study area near the Mississippi River Delta (Weisberg and Liu, 2017), but seasonal wind patterns and large riverine inputs generally drive nearshore currents throughout the region (Smith and Jacobs, 2005). Individuals of each shark species collected from Galveston and Louisiana were pooled for analyses, despite some differences in habitats and environmental conditions between those sampling regions. Our study aimed to describe movement and habitat use patterns across a larger spatial scale (the nwGOM), as opposed to comparing smaller regional patterns.

The species included in this study exemplify the diversity of sharks in the nwGOM. C. leucas are euryhaline, and juveniles utilize low salinity habitats as nursery grounds (Froeschke et al., 2010a; Heupel et al., 2010; Matich and Heithaus, 2015). Late juveniles and adults generally inhabit coastal and offshore waters, but adult females return to estuaries to give birth (Tillett et al., 2012). S. tiburo inhabit estuarine and coastal environments, generally exhibiting residency within a home estuary but using a variety of habitats within coastal and estuarine zones throughout their lifetime (Heupel et al., 2006). In the nwGOM, S. tiburo are found primarily in moderate to high-salinity waters (20–40 ppt) and near tidal inlets (Froeschke et al., 2010b; Plumlee et al., 2018). C. limbatus and C. brevipinna both inhabit bays and coastal estuaries as juveniles, where they prefer moderate to high salinities, warm temperatures, and moderate to deep depths (Ward-Paige et al., 2015; Plumlee et al., 2018). Both species are known to travel long distances within coastal waters (Kohler et al., 1998), and C. limbatus exhibit evidence of reproductive philopatry (Hueter et al., 2005; Keeney et al., 2005). I. oxyrinchus are pelagic and highly migratory, making large-scale movements across ocean basins. I. oxyrinchus mating and pupping ground locations in the GOM are not well understood, but there is some evidence to suggest females give birth offshore (Casey and Kohler, 1992; Gibson et al., 2021).

Vertebrae Collection

Individuals from each species were opportunistically collected in the study region between 2014 and 2017. Depending on the number and quality of samples available, 9–12 individuals per species were chosen for analysis. The selected samples, representing the largest individuals in the collection, were evenly distributed between males and females, and were collected during similar time periods (Table 1). At least five thoracic vertebrae were removed from each individual, and all vertebrae were frozen at −20°C for storage and were thawed prior to processing.

Table 1. Sample sizes (Count) of individuals analyzed with LA-ICP-MS by species and sex, with years of collection, average fork length (FL) ± 1 SD in cm, and range of ages in years (sexes pooled).

Once thawed, vertebral columns were submerged for 20–30 s in boiling water to aid in removing excess muscle and connective tissue. Individual vertebral centra were then separated, cleaned, and dried for at least 24 h. A single cleaned vertebral centrum per individual was then cut along its longitudinal axis using an IsoMet low-speed diamond blade saw (Buehler, Illinois Tool Works Inc.), removing a 2 mm cross-section from the center. The two sides of a sectioned centrum are essentially identical, since the vertebrae grow concentrically outward over time. Cross sections from each sample were thus cut in half to isolate one side of the centrum and were mounted on a glass petrographic slide using Crystalbond thermoplastic cement adhesive (Figure 2).

Figure 2. Example of a sectioned C. leucas vertebral centrum, shown with transmitted light (left) for viewing growth bands (birth band and first growth band shown with red lines), and reflected light (right) for viewing LA-ICP-MS transect path.

Trace Element Sampling

Vertebrae elemental concentrations were measured from sectioned vertebral centra in the direction of radial (outward) growth using an Elemental Scientific NWR193UC (193 nm wavelength, <4 ns pulse width) laser system coupled to an Agilent 7500ce inductively coupled mass spectrometer (LA-ICP-MS) at The University of Texas at Austin. The laser system is equipped with a large format two-volume laser cell with fast washout (<1 s), which accommodated all vertebrae samples and standards in a single loading. Laser ablation parameters were optimized for sensitivity and signal stability from test ablations on representative unknowns: 60% laser power, 10 Hz repetition rate, 25 × 100 μm aperture, 15 μm/s scan rate, He flow of 850 mL/min, and Ar flow of 800 mL/min. Prior to analysis, samples and standards were pre-ablated to remove potential surface contamination. Laser analyses of unknowns were bracketed hourly by standard measurements (USGS MAPS-4, MACS-3 and NIST 612, typically measured in triplicate for 60 s). Baselines were determined from 60-s gas blank intervals measured while the laser was off and all masses were scanned by the quadrupole. USGS MAPS-3 (synthetic bone) was used as the primary reference standard and accuracy and precision were proxied from 39 replicates of NIST 612 and USGS MACS-3 (synthetic aragonite) analyzed as unknowns. Analyte recoveries for NIST 612 and USGS MACS-3 were typically within 5% of GeoREM preferred values. Oxide production rates, as monitored by ThO/Th on NIST 612, averaged 0.34% over the analysis periods. Laser energy densities over the analytical sessions averaged 3.70 ± 0.04 J/cm² for line traverses. The quadrupole time-resolved method measured 11 masses using integration times of 10 ms (^24–25Mg, ⁴³Ca, ⁵⁵Mn, ⁸⁸Sr), 20 ms (⁷Li, ^66–68Zn, ¹³⁸Ba), 25 ms (¹³⁷Ba). The sampling period of 0.24 s corresponds to 89% quadrupole measurement time, with data reporting every 3.63 μm at the scanning rate of 15 μm/s. Time-resolved intensities were converted to concentration (ppm) equivalents using Iolite software (Hellstrom et al., 2008), with ⁴³Ca as the internal standard and a Ca index value of 35 weight% (Mohan et al., 2018). Specifically, the counts-per-second (CPS) of each element were ratioed to the ⁴³Ca CPS at each time point along the transects, which was compared to the known element:Ca ratio in the reference standard. Many studies of otolith and vertebrae chemistry have converted ppm to element:Ca molar ratios, but this was not conducted for these element concentrations since the ppm units were already standardized to ⁴³Ca. Assuming a constant index value of Ca (35 weight%) may influence the reliability of the calculated element concentrations, given there is the potential for variability in Ca within and among vertebrae. However, an examination of Ca CPS along a transect in this study revealed limited Ca variability, which was unlikely to drastically affect element concentration estimates. Future studies of the variation in Ca weight% within and among shark vertebrae would enhance the precision of LA-ICP-MS results.

Aging and Life Stage Determination

Following LA-ICP-MS analysis, digital images of vertebral centra were obtained using transmitted light on a dissecting microscope mounted with a camera to visualize the opaque bands in the corpus calcareum that correspond to growth (growth bands, Figure 2). The age of each shark was determined by counting the visible growth bands, with four independent readers conducting blind counts and subsequently resolving any discrepancies. The birth band was identified as the first growth band accompanied by a change in the growth axis angle. Each individual growth band after the birth band was also identified and marked during the aging process. For all species but I. oxyrinchus, annual deposition of growth bands (one growth band per year) was assumed (McMillan et al., 2017). For I. oxyrinchus, we assumed a deposition rate of two growth bands per year for the first 5 years of life, then annual deposition for the remaining years (Wells et al., 2013; Kinney et al., 2016). The distance (μm) from the start of the LA-ICP-MS laser transect path to the birth band was measured for each transect using ImageJ v1.53a (Schneider et al., 2012), and any prenatal trace element data were excluded from further analyses. The birth band was thus used as the starting point of each transect. Distances (μm) from the birth band to each growth band were then measured, which provided an age (in years) at each measurement along a given transect. Finally, the distance (μm) from the birth band to the visible edge of the corpus calcareum was measured to determine the end of the transect, and all data beyond that point were excluded from analyses.

To examine ontogenetic patterns of elemental signatures, life stages were determined for each shark species based on literature-derived values of age at 50% sexual maturity (Table 2). Age at maturity is sex-specific for C. limbatus, C. brevipinna, and I. oxyrinchus, but not for C. leucas or S. tiburo, which was reflected in life stage determination. We isolated four distinct life stages: young-of-the-year (YOY), early juvenile (EJ), late juvenile (LJ), and adult (AD). The juvenile stage, after surviving as a YOY but before reaching maturity, is prolonged for many of the study species, and can include periods of ontogenetic shifts in habitat use and feeding. To examine this period in greater detail, we partitioned it into EJ (first half of juvenile stage) and LJ (second half) for all species but S. tiburo. Since S. tiburo mature quickly (age 2), the juvenile stage was not partitioned and did not include an EJ stage. Life stage designations were then added to the trace element transect data to align with the distance at each age based on growth band measurements. Since not all individuals of a given species were the same age, the number of individuals within each life stage differed (i.e., some sharks were captured prior to reaching maturity). Sample sizes for each life stage per species are listed in Table 2.

Table 2. Life stage classifications for each shark species by sex: young-of-the-year (YOY), early juvenile (EJ), late juvenile (LJ), and adult (AD).

Statistical Analyses

Six elements were used for statistical analysis: Li, Mg, Mn, Zn, Sr, and Ba. Because we obtained comparable profiles for elements measured with different isotopes, for plotting and statistical analyses we selected the isotope with highest natural abundance (²⁴Mg, ⁶⁶Zn, ¹³⁸Ba). Element concentration values were smoothed to remove the effect of outliers by conducting a 15-point rolling median followed by a 15-point rolling mean along the time-series of each transect, a filter length corresponding to approximately 50 μm of vertebrae material (1–2 months of growth). Smoothed trace element concentration means and 95% confidence intervals were plotted against distance (increasing age) from the birth band to visually and qualitatively compare elemental patterns within and among species.

Prior to statistical analyses of life stages, one mean value was calculated per element per life stage for each individual shark. While this reduces the amount of information provided by the transects, it is necessary to avoid issues with autocorrelation and lack of independence among data points. Since the original data are in the form of a time series (i.e., highly autocorrelated and thus each observation is not independent from the others), using multiple data points within each individual’s life stage would be inappropriate for most statistical tests. A multi-element principal component analysis (PCA) was conducted for each species, with life stage as the grouping variable. This provides a visualization of how vertebrae elemental composition differs among life stages within each species by reducing the multidimensional dataset to two dimensions; namely, by plotting the principal components that explain the most variability in the data (Hotelling, 1933). One PCA was conducted per species using the “prcomp” function in the “stats” package in R version 3.5.1 (R Core Team, 2019), with data points being zero-centered and scaled to unit variance prior to analysis. Although PCAs are valuable for examining overall elemental composition of vertebrae across ontogeny, an Analysis of Variance (ANOVA) framework was also employed to determine differences in individual elements across life stages. A series of ANOVA tests were conducted for each species to determine which elements differed across life stages (n = 6 tests per species, one per element). All models were checked for the assumptions of homoscedasticity and normality of residuals using Levene’s test and Shapiro–Wilks test, respectively. Most models fit these assumptions.

Results

Plots of the mean concentration along the LA-ICP-MS transects revealed species-specific elemental patterns through life (Figure 3). We use distance from the birth band to approximate time from birth to death of each individual. Li concentrations generally decreased over time for all species except C. leucas, which increased over time (from ∼1 to 1.5 ppm, Figure 3A). S. tiburo consistently had the highest Li concentrations of all species (∼2.5 to 2 ppm), followed by C. brevipinna (∼2 to 1.5 ppm), then C. limbatus and I. oxyrinchus (both ∼1.5 to 1 ppm, Figure 3A). Mg concentrations were very similar among species, and generally exhibited a declining trend over time (from ∼4000 to ∼3500 ppm) but with high variability (Figure 3B). Mn concentrations were similar for all species except S. tiburo, which exhibited much higher concentrations that increased throughout life (from ∼75 to ∼200 ppm, Figure 3C). For C. leucas, C. limbatus, and C. brevipinna, Mn concentrations peaked early in life (∼50 ppm), then declined and stabilized for the remainder of life (∼25 ppm, Figure 3C). I. oxyrinchus had the lowest Mn concentrations of all species, which were stable over time (∼10 ppm, Figure 3C). Zn concentrations were also similar among species, and were generally high early in life (∼30 ppm) followed by a gradual decline, ultimately stabilizing late in life (∼20 ppm, Figure 3D). A notable exception was I. oxyrinchus, which had consistently higher Zn concentrations than the other species (from ∼45 to ∼30 ppm) and exhibited large oscillations in Zn late in life (Figure 3D). Sr concentrations were generally stable over time for all species except C. leucas, which exhibited increasing Sr concentration over time (from ∼1750 to ∼2300 ppm, Figure 3E). S. tiburo had the highest Sr concentrations (∼2250 ppm) until the end of life, when C. leucas surpassed them (∼2300 ppm, Figure 3E). The lowest Sr values were seen in both C. limbatus and C. brevipinna (∼1600 ppm), with I. oxyrinchus in the middle (∼1900 ppm, Figure 3E). Finally, Ba concentrations were low and stable for C. limbatus, C. brevipinna, and I. oxyrinchus over time (all ∼10 ppm, Figure 3F). For C. leucas, Ba remained stable and higher than all other species throughout early life (∼30 ppm), then gradually decreased over time (to ∼10 ppm, Figure 3F). The opposite was true for S. tiburo, whose Ba concentrations increased over time and became much higher than all other species late in life (from ∼25 to ∼150 ppm, Figure 3F).

Figure 3. LA-ICP-MS transects for each individual element, with species represented by color. Solid lines represent the mean value of the element from the birth band (distance = 0) to capture for each species individually. Shaded areas represent 95% confidence intervals around the mean.

Differences emerged among the study species in terms of the overall elemental composition of vertebrae among life stages based on PCA (Figure 4). C. leucas life stages separated into two distinct groups: YOY/EJ and LJ/AD (Figure 4A). Principal component 1 (PC1) explained the majority (65.2%) of the variation among data points, while PC2 explained much less variation (14.8%). The elements with the most influence along PC1 were Sr and Li (positive loadings), while Mn, Ba, Mg, and Zn had negative loadings for that axis. S. tiburo had almost no separation among life stages along either PC axis (Figure 4B). PC1 explained 47.5% of the variation among data points, while PC2 explained 23.8%. Loadings were negative along PC1 for all elements. C. limbatus exhibited distinct separation of the YOY life stage, while LJ and AD overlapped nearly completely (Figure 4C). The EJ ellipse fell between YOY and the LJ/AD grouping. PC1 explained 40.1% of the variation among data points, while PC2 explained 21.8%. The elements with the most influence along PC1 were Ba (positive loading), and Zn and Li (negative loading), which appeared to explain much of the variation among life stages in contract to values of Mg and Sr, whose vectors aligned with variation within each life stage. C. brevipinna exhibited a similar pattern to C. limbatus, with YOY separating from the other life stages, LJ and AD overlapping, and EJ in the middle (Figure 4D). However, the separation between these groupings was less defined, with EJ overlapping more with the LJ/AD group. PC1 explained 34.7% of the variation among data points, and PC2 explained 28.5%. The most influential elements along PC1 were Ba (positive loading), and Mn and Li (negative loadings). Along PC2, Mg and Zn exhibited positive loadings while Sr was negative. For I. oxyrinchus, the YOY and EJ ellipses overlapped, and the LJ overlapped slightly with all four life stages (Figure 4E). AD was separate from YOY and EJ. PC1 explained 46.6% of the variation among data points, while PC2 explained 19.4%. All elements besides Ba influenced PC1 (negative loadings), while Ba exhibited a negative loading along PC2.

Figure 4. Visualized principal component analysis for each species, with life stage (YOY, young-of-the-year; EJ, early juvenile; LJ, late juvenile; AD, adult) as the grouping factor (represented by color). Principal components 1 and 2 are represented as the x and y axis, respectively, and the percent of total variation explained by each PC is listed. Loadings for each element are represented by black arrows, and are labeled per element.

Analysis of individual element concentrations revealed species-specific ontogenetic patterns across life stages (Figure 5). Statistically significant results are reported here, while a full summary of results can be found in Table 3. For C. leucas, five elements differed significantly among life stages: Mg (p = 0.035), Mn (p < 0.001), Zn (p = 0.022), Sr (p < 0.001), and Ba (p < 0.001, Table 3 and Figures 5B–F). For S. tiburo, no elements differed significantly among life stages (Table 3 and Figures 5A–F). The mean concentrations of four elements differed significantly among life stages for C. limbatus: Li (p < 0.001), Mn (p < 0.001), Zn (p < 0.001), and Ba (p < 0.001, Table 3 and Figures 5A,C,D,F). The mean concentrations of three elements differed significantly among life stages for C. brevipinna: Li (p < 0.001), Mn (p < 0.001), and Ba (p < 0.001, Table 3 and Figures 5A,C,F). For I. oxyrinchus, four elements differed significantly among life stages: Li (p < 0.001), Mg (p = 0.025), Mn (p < 0.001), and Zn (p = 0.011, Table 3 and Figures 5A–D).

Figure 5. Boxplots visualizing differences among life stages for each element and species, with life stage represented by color. Boxes show the 25th percentile, median, and 75th percentile, with whiskers representing 1.5 times the interquartile range. Outliers (data points that fall outside the range of the whiskers) are represented by black circles. Asterisks represent species with statistically significant differences among life stages (ANOVA). The asterisks (*) indicate statistical significance.

Table 3. Results from ANOVA tests comparing mean elemental values among life stages for each species and element independently.

Discussion

Patterns of element concentrations in the vertebral cartilage of five shark species may be linked to their movement and habitat use patterns through life. The documented ontogenetic patterns of vertebral element composition can be partly explained by known life history traits. However, some unexpected species-specific elemental concentration patterns may be linked to physiological and biological controls, as opposed to purely reflecting ambient environmental conditions. While much is left to learn regarding the pathways and rates of element uptake in elasmobranch vertebrae, the results presented here highlight the applicability of this technique as a method for understanding the life histories of estuarine, coastal, and oceanic sharks.

The most euryhaline of the study species, C. leucas exhibit a complex life history that was reflected in ontogenetic patterns of vertebral element concentrations. A clear difference in overall element composition was observed between early life (YOY and EJ) and later life (LJ and AD), and all but one of the elements (Li) differed significantly across life stages. Since C. leucas are known to use freshwater regions in estuaries as nursery grounds early in life (Heupel et al., 2010), many of the observed patterns of element composition in their vertebrae can be explained by changes in their exposure to freshwater as they age. Sr and Ba generally increase and decrease, respectively, with increasing salinity (Dorval et al., 2005), and are considered fairly reliable proxies for salinity history in elasmobranch vertebrae (McMillan et al., 2017). Our results suggest that C. leucas inhabit low-salinity estuarine regions throughout the first 4–5 years of life, followed by increasing use of marine habitats as they transition to adulthood. Similar elemental patterns have been observed in the vertebrae of C. leucas in Australian coastal waters, including evidence suggesting periodic returns to freshwater habitats for pupping in adult females (Tillett et al., 2011; Werry et al., 2011). The consistent ontogenetic shift in Sr and Ba concentrations in C. leucas vertebrae within multiple ocean basins lends support to the use of this method as a tracer of salinity history in euryhaline elasmobranchs. Additionally, since C. leucas exhibit movements between vastly different water types with different trace element composition (fresh, brackish, and seawater), they are an ideal species for vertebral chemistry due to the relative clarity of trends across ontogeny.

In contrast to the clear ontogenetic patterns of element composition in C. leucas vertebrae, S. tiburo exhibited very little change in elemental signatures throughout their lifetime. Two elements, Mn and Ba, increased near the end of the S. tiburo life history transects, but this trend was inconsistent among individuals and not supported by statistical comparisons among life stages. S. tiburo occupy estuarine and coastal habitats and are thought to be resident within a home estuarine/coastal region (Heupel et al., 2006). While there is evidence that adult S. tiburo migrate seasonally along the United States. East Coast (Driggers et al., 2014), this is unresolved in the nwGOM. Furthermore, S. tiburo distribution does not appear to differ between seasons (spring and fall) along the Texas coast (Froeschke et al., 2010b), suggesting long-distance migration may not occur in this region. The lack of clear ontogenetic changes in vertebrae element composition may indicate that S. tiburo in the nwGOM mainly reside in specific estuarine and coastal regions, and their habitat use does not change dramatically across ontogeny. Comparing ontogenetic changes in vertebrae trace elements of S. tiburo from multiple locations throughout the GOM and western Atlantic may help to characterize regional differences in migration patterns. S. tiburo also exhibit natal philopatry, meaning adults return to regions near their own birthplace, and their population is highly structured in the western Atlantic and GOM based on genetic markers (Escatel-Luna et al., 2015; Portnoy et al., 2015; Gonzalez et al., 2019). Collectively, the apparent limited movement and dispersal patterns of S. tiburo in the nwGOM could mean that they are prone to localized depletion via anthropogenic impacts (Walker, 1998; Hueter et al., 2005). Vertebral chemistry may be a useful tool in identifying species and populations that are prone to similar risks.

Carcharhinus brevipinna and C. limbatus exhibited similar life history concentration patterns for most elements examined, suggesting their habitat use patterns likely overlap throughout ontogeny. Both species use coastal bays and estuaries as nursery grounds, and they often co-occur in the same habitats (Castro, 1993; Hueter and Manire, 1994; Parsons and Hoffmayer, 2007). In coastal environments, juveniles and subadults of both species exhibit preferences for moderate to high salinities (20–30 ppt), warm temperatures (>25–30°C), and close proximity to tidal inlets (Ward-Paige et al., 2015; Plumlee et al., 2018). Stable post-YOY elemental concentration patterns for both species suggest preferences for specific environmental conditions that do not change drastically with ontogeny. However, the notable observed shift in elemental concentrations following the YOY life stage for both species may indicate movement away from nursery areas after their first year of life. Adult C. limbatus and C. brevipinna are known to be highly mobile, traveling long distances along coastlines in the GOM and western Atlantic (Castro, 1993; Kohler et al., 1998; Kajiura and Tellman, 2016), which may explain the changes in trace element profiles following their YOY stage. Considering the clear similarities in vertebral chemistry between C. limbatus and C. brevipinna in this study, the likelihood of these two species sharing the same habitats throughout ontogeny is high, which may influence their ecological roles in their shared habitats (Matich et al., 2017b). Further study of the overlap patterns between them, including habitat use and diet, would enhance our understanding of how these similar species coexist.

Similar to C. limbatus and C. brevipinna, I. oxyrinchus life history transects were generally stable over time, but ontogenetic shifts in element composition and concentrations emerged when comparing life stages. The two elements most closely linked with salinity, Sr and Ba, did not differ across life stages, suggesting that I. oxyrinchus do not inhabit estuarine or low-salinity regions like river plumes. This result is consistent with currently known habitat use patterns of this species: They generally inhabit coastal and pelagic zones, which have low environmental variability compared to highly dynamic estuaries. The general stability of element concentrations along I. oxyrinchus life history transects may reflect a relatively narrow range of suitable environmental conditions, such as temperature and salinity (Vaudo et al., 2016). It is also important to note that I. oxyrinchus are the only lamnid shark capable of regional endothermy examined in this study (Wolf et al., 1988). The relative stability of temperature within the bodies of I. oxyrinchus may contribute to stability in elemental concentrations if temperature exerts control over element uptake into vertebrae. Tagging studies have provided some clarity on the movement and habitat use patterns of I. oxyrinchus in the GOM, which may explain the observed differences in vertebrae elemental composition across life stages. For example, juvenile I. oxyrinchus tagged off the coast of Mexico appear to remain in the GOM year-round, largely within continental shelf habitats (Vaudo et al., 2017). Conversely, sexually mature adult I. oxyrinchus in the GOM are documented making long-distance migrations, as far as northeastern United States. Atlantic coastal waters in the summer and returning to the GOM in the winter (Gibson et al., 2021). Differences in migration patterns between juvenile and adult I. oxyrinchus in the GOM may have contributed to the differences in elemental signatures observed between life stages.

Although many of the observed differences in trace element concentrations among these five species could be explained by differences in life history and habitat use, some elements showed unexpected patterns. For example, Sr concentrations are thought to be associated with salinity history, but were only moderate for I. oxyrinchus which are the most oceanic of the study species. We would expect their Sr concentrations to be the highest of all the study species, since seawater generally holds higher Sr concentrations than fresh or estuarine waters (Dorval et al., 2005). Similarly, S. tiburo vertebrae had elevated concentrations of specific elements (Li, Mn, Sr, Ba) compared to other study species that reside in similar coastal and estuarine regions. Element uptake in these sharks may therefore be influenced by species-specific physiological and biological mechanisms in addition to reflecting ambient water concentrations and environmental history. In the case of S. tiburo, we hypothesize that the elevated concentrations of Li, Mn, Sr, Ba may be related to the difference in diet between S. tiburo and the other study species. S. tiburo consume primarily benthic prey such as crabs and shrimp (Bethea et al., 2007; Plumlee and Wells, 2016; Kroetz et al., 2017), while the other study species consume fish and other pelagic prey (Stillwell and Kohler, 1982; Bethea et al., 2004; Plumlee and Wells, 2016; TinHan, 2020). Trace element uptake in elasmobranch vertebrae is thought to be at least in part driven by dietary sources (Mathews and Fisher, 2009), so differences in diet and water column position (benthic vs. pelagic) could potentially influence differences in vertebral chemistry among species. Elucidating the species-specific mechanisms that influence the concentrations of trace elements would enhance the utility of vertebral chemistry as a method for direct comparisons of habitat use patterns among species.

Another element with an unexpected pattern was Mg, as concentrations were very similar and exhibited similar patterns of decreasing concentration over time for all species regardless of differences in known habitat use patterns among them. One experimental validation study using Round Stingrays (Urobatis halleri) concluded that Mg concentrations in vertebrae decreased with increasing temperature, and that uptake was not mediated by somatic growth rates or vertebral accretion rates which also increase with temperature (Smith et al., 2013). The maximum age of the sharks in this study was 17 years (C. brevipinna), and water temperatures in the nwGOM have increased over that timescale (by approximately 0.05°C per year from 1963 to 2015, Turner et al., 2017). It is therefore possible that the consistent decrease in Mg in the vertebrae of these coastal and oceanic sharks is reflecting temperature increases due to global climate change. The relationship between Mg uptake and environmental conditions requires further study, but these results represent the possibility of reconstructing water temperature changes over time in long-lived elasmobranchs.

Similarities also emerged among most of the species in their patterns of Mn and Zn concentrations, with both elements being elevated early in life followed by a decline. A similar pattern was observed for multiple shark species from Australia, including C. brevipinna, where Zn concentrations were very high pre-birth and declined post-birth (Raoult et al., 2018). The elevated Zn concentrations observed in the early life of all species in this study match the post-birth pattern observed by Raoult et al. (2018), and we hypothesize that this may be the result of maternal loading. Maternal loading is the process by which pregnant females transfer elements and molecules to their embryos via lipid mobilization (Addison and Brodie, 1987). This process is known to occur in sharks with regard to trophically derived contaminants in liver and muscle tissue (Lyons et al., 2013; Mull et al., 2013), and trace elements in the vertebrae of shark embryos (Coiraton and Amezuca, 2020). Both Zn and Mn are thought to be primarily incorporated into elasmobranch vertebrae via dietary sources as opposed to ambient environmental conditions (Mathews and Fisher, 2009), so it is reasonably likely that these trophically derived elements would undergo maternal loading during embryonic development. Additionally, Mn has been linked to dissolved oxygen levels in the otoliths of teleost fish (Limburg and Casini, 2019) and bivalve shells (Schöne et al., 2021). It is therefore possible that exposure to hypoxia may have influenced elemental concentrations in the vertebrae of the sharks in this study. Hypoxic conditions are widespread and highly variable in the nwGOM (Rabalais and Turner, 2019), but the influence of hypoxia exposure on shark vertebrae chemistry remains unresolved.

The use of vertebral chemistry can enhance our understanding of habitat use patterns through ontogeny in elasmobranchs, but interpretation of element concentrations would benefit from experimental validation of uptake routes and rates among taxa. We are aware of only two published experimental studies to date that have examined the effects of environmental conditions (temperature, pH, water chemistry) on element uptake in elasmobranch vertebrae, including one stingray and one oviparous shark (Smith et al., 2013; Pistevos et al., 2019). The lack of information regarding the complexities of environmental, physiological, and biological controls of trace element uptake limits our ability to make inferences about habitat use and movement. Experimental validation is likely hindered in part by the difficulties inherent to husbandry and maintenance of animals that are as slow-growing, long-lived, and large as many of the shark species in the nwGOM. S. tiburo are relatively small sharks with rapid time to maturity, which makes them an ideal candidate for experimental validation of the effects of extrinsic (environmental conditions) and intrinsic (physiological and biological) factors on trace element incorporation in shark vertebrae. An alternative to controlled tank studies is mark-recapture approaches that have been applied to sharks in the eastern Pacific (e.g., Mohan et al., 2018). Future validation studies would strengthen our ability to interpret the life history transects presented here, which would lead to a better understanding of the movements and habitat use patterns of coastal and oceanic sharks.

The results of this study exemplify the diversity of habitat use, movement, and migration patterns of sharks in the nwGOM, and provide evidence that vertebral chemistry is a promising method for reconstructing aspects of shark life history. Examining the movements of large, long-lived animals is extraordinarily challenging (Block et al., 2011), but chemical analysis of vertebrae eliminates the issues of long-term tracking (as with telemetry studies) or high recapture rates (as with conventional tagging). There is still considerable work to be done to validate the mechanisms of trace element uptake into elasmobranch vertebrae, but those studies would allow us to revisit our results and shed new light on the life histories of these sharks. Anthropogenic impacts including fishing mortality and habitat loss have resulted in declining populations of many shark species in the GOM (Baum and Myers, 2004), and the IUCN lists all of the species in this study as either near threatened (C. leucas and C. limbatus, Burgess and Branstetter, 2009; Simpfendorfer and Burgess, 2009), vulnerable (C. brevipinna, Rigby et al., 2020), or endangered (I. oxyrinchus and S. tiburo, Rigby et al., 2019; Pollom et al., 2020). Conservation and management of these species will require understanding the intricacies of their habitat requirements and movement patterns, especially for those that migrate through multiple jurisdictions (Rooker et al., 2019). Trace elements in vertebral cartilage, especially following validation studies, may prove to be instrumental in conserving populations of threatened elasmobranchs.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

Ethical review and approval was not required for the animal study because this study used pre-existing samples. No animals were handled during the course of this study.

Author Contributions

ML and RW conceived the study. BF collected and provided shark samples. ML, JM, and NM conducted LA-ICP-MS sampling. ML led the statistical analyses with assistance from TT, JM, TR, and RW and drafted the manuscript. JM, TT, and TR created the figures. All authors contributed equally to manuscript revisions.

Funding

This work was supported by a grant to ML by Texas Sea Grant Grants-in-Aid of Graduate Research (GIA-2019-2021).

Conflict of Interest

BF was employed by Fisheries Research Support, L.L.C.

The remaining 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

We thank the many volunteers at Texas A&M University at Galveston who helped to collect shark samples from Galveston docks, and Louisiana Department of Wildlife and Fisheries employees C. Marshall and P. Kent who helped to collect shark samples from Louisiana docks. Special thanks to E. Gutierrez and M. Curtis for assisting with sample collection and processing, and to R. Ware (NSF Research Experience for Undergraduate Students 2019) for helping with vertebrae aging.

References