Geochemical signatures and nanomechanical properties of echinoid tests from nearshore habitats of Florida: environmental and physiological controls on echinoid biomineralization

Study systems and data collection

Live echinoid specimens were collected via SCUBA expeditions from 15 sites (Table 1) in three regions along the western and southern coasts of Florida: (1) central Florida Keys (FK), (2) Cedar Key (CK), and (3) Apalachee Bay at the northern Gulf coast of Florida (NG) (Fig. 1). All surveying and collecting activities were carried out within the scope of the collecting permits SAL-18-1294A-SR, SAL-19-2195-SR, SAL-19-2195A-SR, SAL-22-2195-SR, and SAL-22-2195A-SR issued by the Florida Fish and Wildlife Conservation Commission. A total of 43 specimens were analyzed for bulk geochemical analysis and four of those specimens were also subject to nanoindentation and fine scale microprobe analyses.

The three regions notably differ in salinity and, to a smaller extent, in temperature. FK is characterized by the highest and least variable salinity and temperature, NG is characterized by relatively lower and more variable salinity and temperature, and CK sites, influenced by Suwannee estuary, have the lowest salinity (Fig. 2) and temperature comparable to FK. The FK sites represent carbonate sand bottoms, whereas CK and NG sites represent mixed siliciclastic-carbonate sand bottoms. All sites occur in shallow-water, coastal areas with water depth ranging from 1 to 21 m (Table 1). The offshore site off Cedar Key (site CK-2) is outside the direct influence of Suwannee river and thus more comparable to other Northern Gulf sites. Salinity estimates for CK-2 were obtained from a comparably offshore, non-estuarine site off Crystal River with salinity estimates consistent with an in situ measurement of 31‰ collected at CK-2 during field work. In all analyses, site CK-2 is classified as ‘Northern Gulf’ (Table 1). The shallower offshore site (CK-5) was classified as ‘Cedar Key’ with salinity the same as CK-6. The alternative results, with CK-5 assigned to “Northern Gulf” using the same salinity estimates as for the nearby CK-2, are also reported below. In addition, the CK-5 specimens are highlighted with special symbols on the key figures included below. Reassigning regional membership and salinity for sites CK-5 and CK-2 changes slightly some of the plots and values of statistics reported below, but does not change any of the major conclusions of this study. Reassignments of CK-5 and CK-2 to different regions and salinity estimates were performed interactively during the analysis (see R script in Data S1), whereas the values archived in Data S2 are original values prior to the reassignment.

Four irregular echinoid taxa were selected based on their common presence in the studied regions and adequately large test size. All taxa included here possess an endoskeleton (test) consisting of interlocking plates made of high-Mg calcite.

Clypeaster subdepressus (family Clypeasteridae) is a relatively large echinoid common in soft-bottom sandy habitats of the Gulf of Mexico. The species is a mobile, semi-infaunal to epifaunal detritivore. During collecting, the specimens were typically observed only partly buried under sediment surface with the aboral side protruding above sediment surface (see also Telford, Mooi & Harold, 1987). C. subdepressus ingests coarse sand-sized particles and mostly feeds on dead algae, seagrass, and coral fragments (Telford, Mooi & Harold, 1987).

Encope spp. include two species of sand dollars: E. michelini and E. aberrans (family Mellitidae). Juveniles and subadults of these two species are difficult to distinguish morphologically and the molecular phylogeny indicates that they represent closely related species (Kroh & Mooi, 2024) that split ∼6 million years ago (Coppard & Lessios, 2017). These two species are similar in size, mode of life, feeding, and habitat preferences. They are common in soft-bottom sandy habitats of the Gulf of Mexico and often occur sympatrically. Similarly to C. subdepressus, they are mobile, semi-infaunal to epifaunal detritivores. The fossil record of the two species dates back to the middle-to-late Pleistocene (Coppard & Lessios, 2017).

Mellita tenuis is a sand dollar (family Mellitidae) common in the eastern Gulf of Mexico, where it often occurs in dense congregations (e.g., Cleveland & Pomory, 2022). The species is a mobile epifaunal to shallow infaunal detritivore that ingests sand-sized grains. The species is known from the Pleistocene of Florida (Coppard, Zigler & Lessios, 2013).

Leodia sexiesperforata is a sand dollar (family Mellitidae) that is widespread in the Gulf of Mexico. This species is a mobile, infaunal detritivore that occurs in shallow-water soft-bottom sandy habitats, often in dense congregations (Grun & Kowalewski, 2022). The species is known from the late Pleistocene fossil record of Florida (Mooi & Peterson, 2000; Portell & Oyen, 2002) and Jamaica (Mitchell, James & Brown, 2006).

Environmental estimates (salinity and temperature)

The temperature estimates used here were downloaded from NOAA (2024) using seasonal datasets with the highest available resolution (0.1°). To derive the most relevant temperature estimates, the NOAA sectors nearest to our sampling sites were selected based on the nearest geographic distance calculated using R package ‘geosphere’ (Hijmans, 2022). The distances varied from ∼2 to ∼40 km. Because all sites represent shallow water settings with fully mixed waters, water temperatures were effectively depth invariant within sites and surface water temperatures can serve as proxies regardless of the site’s depth. This assumption was confirmed by the analysis of the NOAA (2024) database: the temperature range in 0-to-25 m depth range averaged 0.05 °C and only 5% of NOAA sectors exceeded 1 °C difference within the top 25 m of the water column. Annual temperature was estimated as the arithmetic mean of the four seasonal means, while range of seasonal means was used as a measure of annual variability in temperature (Fig. 2; Data S2).

Salinity data were downloaded from Gulf of Mexico Ocean Observation (GCOOS) (2019), an online-accessible aggregator of salinity data for the Gulf of Mexico. The GCOOS sites most proximal to our sampling sites were selected based on the nearest geographic distance calculated using R package ‘geosphere’ (Hijmans, 2022). Mean and median salinity values were calculated based on all measurements available for the most proximal GCOOS station (Data S2). The distance varied from ∼0.4 to ∼21 km. The variation in salinity was measured as the interquartile range of all available measurements.

Bulk geochemical analyses

All analyzed specimens (n = 43) were live-collected individuals kept on ice and subsequently fixed in 95% scientific-grade isopropanol. Approximately 1/5 of each specimen (posterior sector cut through ambulacral lunules and running through the center of the petals, excluding lantern) (Fig. 3) were shortly bleached with sodium hypochlorite and split out into smaller fragments. These fragments were mechanically cleaned for the presence of soft tissue remnants, bleached again, and then washed with the aid of Milli-Q Water in the ultrasonic bath. These samples were finally dried at 50 °C, powdered and homogenized using an agate mortar and pestle. The resulting 43 aliquots (∼1 g) were analyzed separately for each specimen using the AQ250-EXT package (ultratrace aqua regia/ICP-AES and MS; Inductively Coupled Plasma Emission Spectrometry-Mass Spectrometry) at a commercial geochemical laboratory (Bureau Veritas Minerals). The entire set of analyzed elements (n = 53), along with their detection limits, can be found in Data S3. However, in this paper, we focus on discussing a subset of elements (Mg, strontium (Sr), barium (Ba), sulfur (S), lithium (Li), lead (Pb), zinc (Zn), nickel (Ni), cobalt (Co), manganese (Mn), boron (B), cadmium (Cd), sodium (Na)) that have often been used as environmental proxies (e.g., Iglikowska et al., 2020; Ulrich et al., 2021) and/or yielded estimates above the detection limits. All elements were converted to molar concentrations, and, as is commonly practiced, normalized to Ca concentration and expressed as element/Ca (mmol/mol). For the nanoindentation analysis described below, Mg/Ca ratios were reported using mol/mol ratios. Data quality was monitored by blind insertion of sample duplicates, internal reference materials, and the certified reference materials (including STD OREAS262 and STD BVGEO01). Analytical accuracy was within 2.5% of the value of the certified reference material for Mg, Sr, S, Na, P and within 5% of the standard value for Mn, Pb, Ba, Li, Cd, Cu and Zn.

Because the specimens were subject to destructive analyses, the remaining fragmentary material was not archived in a formal museum collection. However, it is available upon request from the first author.

Nanoindentation and fine scale microprobe analyses

Two specimens of Leodia sexiesperforata and two specimens of Encope michelini were selected for nanoindentation analyses. One fragment of each test, representing the peripheral side posterior sector, were embedded in epoxy resin, cut perpendicularly to the test margin through the interambulacral plates, and ultra-polished. A nanoindenter (NanoTest Vantage; Micro Materials, Wrexham, UK) with a Berkovitch tip was used to determine nanomechanical properties (i.e., nanohardness [H] = the material’s resistance to permanent or plastic deformation at the nano-micro level, in GPa) of the stereom following Oliver & Pharr (1992). The maximum load of 2 mN, resulting in the maximum penetration depth of ∼200 nm, was applied. Three microstructurally distinct regions from each specimen were analyzed separately, namely (1) the outermost solid, imperforate to perforate stereom layer (tubercles), (2) inner galleried stereom, and (3) inner massive and coarse labyrinthic stereom. For each stereom type, three microregions were analyzed resulting in nine separate analyses for each specimen (three analysis per region, three regions per specimen) for a total of 36 nanoscale analyses (Data S4, S5). Following the nanoindentation analysis, spot elemental measurements (Mg/Ca) in the same microregions were determined (Data S4, S5) with the aid of wavelength-dispersive spectroscopy (WDS) performed on a CAMECA SX100 electron microprobe at the Micro-Area Analysis Laboratory, Polish Geological Institute-National Research Institute in Warsaw, Poland. The following conditions were used: accelerating voltage: 15 kV; beam current: 5 nA for calcium and 20 nA for other elements; a beam diameter: ∼5 µm; standard: ‘NIST’ (serial number: 12570). Detection limit for Mg under this method was ∼0.01 wt %.

Analytical methods

The relatively small sample sizes, which reflect prohibitive costs and the time-consuming nature of both nano-scale analyses and in situ collecting of often sparse echinoid populations via SCUBA, limits the statistical power and application of multi-factor models, which otherwise would have been highly appropriate here. Consequently, only simple descriptors and tests are feasible, including simple descriptors of central tendency, two-sample tests based on medians, bivariate correlation and semi-partial correlation analyses (partialling out interactions between salinity, temperature, and depth), and exploratory ordination methods. Correlation analyses (Pearson r) were limited to assessment of individual trace element ratios against extrinsic, independently derived environmental variables (salinity, temperature, depth). Semi-partial correlations (Pearson r_SP) were estimated separately for each trace element ratio against the three extrinsic environmental variables partialled out against each other. Semi-partial correlations were performed using the R package “pcor” (Kim, 2015).

Multivariate analyses included exploratory principal component analysis (PCA) analyses based on z-standardized data (i.e., correlation PCA) to remove the effect of huge differences in variance across variables. PCA analysis was supplemented by correlation plots, loading plots, and bootstrap analysis of eigenvalues. The criterion 0.7*L, where L is a sum of proportional eigenvalues divided by number of variables (i.e., mean proportional eigenvalue), was used to decide which principal components should be considered in the interpretation (Jolliffe, 2002). It should be noted here that PCA is a parametric method that is not optimal for compositional data with variables affected by constant-sum-constraint (an issue further amplified here by Ca normalization). However, for high-dimensional data (as is the case here: 53 variables total) those effects tend to be less severe (Kucera & Malmgren, 1998). We selected PCA analysis mainly for comparative purposes as this method was previously employed for analyzing trace element ratios in echinoderm skeletons (Iglikowska et al., 2020). However, to evaluate the robustness of the results, the analysis was supplemented with nonmetric multi-dimensional scaling (NMDS). The NMDS analysis was performed on z-standardized variables using Euclidean-distance dissimilarity matrix and k = 2 dimensions (the stress <0.2 was used as our criterion for acceptable stress). NMDS was performed using the function metaMDS in the package ‘vegan’ (Oksanen et al., 2024). As shown below, the NMDS and PCA ordinations produced highly consistent ordination patterns.

In the case of nanoscale analyses, two-sample comparisons based on single factors (echinoid skeletal zone and echinoid species) were applied. Whereas single-factor analyses ignore confounding effects of other factors (e.g., taxon), the balanced sampling design mitigates this issue. For example, in the case of nano-scale analyses, both taxa include the same set of measurements across the same set of three microstructural zones, and outer and inner zones have the same number of measurements for each taxon. In addition, in the case of the nano-scale analysis, the two taxa came from different localities/environments and thus effect ‘taxon’ may reflect environmental parameters or vital effects (a conundrum explicitly considered in the interpretations below). These issues impose limitations on interpretations of effects and significance tests. Pearson correlation coefficient was used to assess the strength of association between Mg/Ca ratios and nanohardness.

To ensure conservative statistical interpretations, a Bonferroni correction was applied when appropriate with significance alpha [α] set to Bonferroni-corrected value of α_B = 0.05/k, where k is the total number of tests of a given type performed in the study. All figures (except Fig. 3) and analyses were performed in R version 4.3.1 (R Core Team, 2023). R script is provided in Data S1 and custom functions for PCA used in the script are provided in Data S6, S7. Raw bulk geochemical data (not used in the analyses) are provided in Data S3. SEM images of echinoids showing sampling points for nanoscale analyses are included in Data S5.

Trace element analysis

Correlations between the analyzed element ratios, salinity and temperature varied across elements (Table 2; Fig. 4). However, salinity and temperature were strongly correlated with each other across the sites (r = 0.84, p < 0.0001). In addition, water depth was weakly associated with lower temperature (r =  − 0.31, p = 0.04), but not salinity (r =  − 0.05, p = 0.74). After partialling out interactions between salinity, temperature, and water depth (Table 2), no notable or significant semi-partial correlations are observed for temperature (Table 2). However, after accounting for the temperature and water depth effects, S and Li still displayed positive and significant or marginally significant semi-partial correlations with salinity (Fig. 4; Table 2) with a notable effect size (r_SP > 0.5). However, these correlations decrease and cease to be highly significant (Table 3) when salinity for site CK-5 is reassigned (see ‘Material and Methods’). The subsequent analyses focus primarily on salinity (Fig. 4).

Comparisons of taxa within regions suggested that sympatric species from Florida Keys, where salinity was highest, differed in geochemical signatures for some elements. This was particularly notable for Mg/Ca and Li/Ca ratios in Florida Keys, with all L. sexiesperforata specimens having much higher values when compared to Encope spp. and C. subdepressus specimens (Figs. 4A, 4E). Similar, if less pronounced offsets were observed for Ba/Ca (Fig. 4C), Pb/Ca (Fig. 4F), Zn/Ca (Fig. 4G), B/Ca (Fig. 4J), and Cu/Ca (Fig. 4M). In contrast, the offsets in those geochemical signatures were notably smaller or absent in the case of specimens collected from the lower salinity setting of the Northern Gulf.

Within-species comparisons across different regions suggest consistent offsets in multiple element ratios. L. sexiesperforata from Florida Keys displayed higher trace element ratios (except for Mn/Ca ratio) when compared to conspecific specimens from the lower salinity settings in Northern Gulf, and in many cases those offsets were significant (Table 4). Similarly, M. tenuis displayed consistently higher Mg/Ca (Fig. 4A), Sr/Ca (Fig. 4B), S/Ca (Fig. 4B), Li/Ca (Fig. 4E), and Pb/Ca (Fig. 4K) ratios in Northern Gulf compared to the lower salinity sites at Cedar Key (although a formal analysis could not be carried out due to small sample size and uncertainty related to the salinity estimates for the site CK-5). In some cases, high intraspecific variation was observed within the same region/salinity regime, especially for Ba/Ca and Zn/Ca ratios in L. sexiesperforata (Figs. 4C, 4G).

Despite confounding effects of interspecific and intraspecific variability, when data are pooled across taxa and the specimens from the high-salinity region of Florida Keys is compared to specimens from the lower salinity regimes of Northern Gulf and Cedar Key (combined), the majority of trace element ratios are higher for Florida Key specimens (Table 5), including Mg/Ca (Fig. 4A), Sr/Ca (Fig. 4B), S/Ca (Fig. 4D), Li/Ca (Fig. 4E), Pb/Ca (Fig. 4F), Zn/Ca (Fig. 4G), B/Ca (Fig. 4J), and Cu/Ca (Fig. 4M). The Mn/Ca ratios are a notable exception showing a reverse trend (Fig. 4H; Table 5).

Principal component analyses of z-standardized elemental ratios produced an ordination consistent with correlation analyses by highlighting the distinctness of environmental settings (the three regions) and the presence of inter- and intraspecific variation within environmental settings (Fig. 5A). The first two principal components accounted for 51% of variance in the data and variables were moderately or strongly correlated either with PC1 or PC2, or both components (Fig. 5B). The long arrows that approach the outer blue circle (‘circle of correlation’; sensu Abdi & Williams, 2010) indicate the variables for which almost all variance was captured by the first two principal components. Conversely, the loading plot indicated that the variance explained by the first two components was partitioned relatively evenly across all variables (Fig. 5C). The bootstrap analysis of eigenvalues suggested that the first four principal components may be potentially informative (Fig. 5D). The correlation plot (Fig. 5B) indicated that specimens with high PC1 scores tended to have higher elemental ratios for all elements except Ba and Mn that were uncorrelated with PC1 scores. Specimens with high PC2 scores had low Ba/Ca and high Mn/Ca scores. PC2 scores were also negatively correlated with Pb/Ca and Zn/Ca ratios and showed less pronounced positive relationships with S/Ca, Sr/Ca, and Li/Ca ratios (Fig. 5B). Among the external variables, salinity displayed strong correlation with PC1 (r =  − 0.74, p < 0.0001). Temperature also displayed notable correlation with PC1 (r = 0.61, p < 0.0001), but based on semi-partial correlation analysis (Tables 2 and 3), this correlation reflected co-variation of temperature with salinity. All other correlation coefficients of environmental variables with PC1 and PC2 were <0.5 thus accounting at most for 25% of variance in the PC scores.

The ordination (Fig. 5A) indicated that the three regions were distinct in geochemical signature, with the three environmental settings largely separating along PC1 indicating that most element ratios increased with salinity. High intra-regional, intra-specific variability was observed for L. sexiesperforata from Florida Keys. Intra-regional interspecific differences in element ratios were evident in multiple cases, including differences between L. sexiesperforata and Encope spp. in Florida Keys and M. tenuis and C. subdepressus in Cedar Key (Fig. 5A). Specimens of the same species from different regions differed in geochemical signatures in the case of M. tenuis and L. sexiespeforata (Fig. 5A). In contrast, all five specimens of C. subdepressus plotted in the upper central part of the ordination indicating relatively comparable geochemical signatures regardless of the environmental setting.

Whereas PC3 and PC4 axes may be potentially informative (Fig. 5D), ordination plots of those higher components (Figs. 6B and 6C) do not indicate any notable changes to patterns observed in the first two dimensions (Figs. 5A or 6A). The NMDS ordination (Fig. 6D) is remarkably consistent with the PCA ordination indicating that the results are reproducible using multivariate ordination methods that are methodologically disparate.

A comparison of Mg/Ca and S/Ca, illustrates all above mentioned patterns for the two element ratios (Fig. 7), highlighting consistent interregional offsets within species (M. tenuis for Northern Gulf versus Cedar Key and L. sexiesperforata for Florida Keys versus Northern Gulf), high within-species within-region variability (L. sexiesperforata in Florida Keys), interspecific differences within regions (L. sexiesperforata versus other species in Florida Keys), and intraspecific invariance across regions with different salinities (Encope spp. and C. subdepressus showing similar values for Florida Keys and Northern Gulf).

Although there is a substantial variation in body size, both among and within the four echinoid taxa, the body size was not correlated notably or significantly with element ratios after accounting for interactions with environmental variables (semi-partial correlation r ≤ 0.3, p > 0.05, in all cases). Consistently with those results, specimens in the PCA ordinations (Fig. 5) did not ordinate according to size. For example, the two species with smaller body sizes (M. tenuis and L. sexiesperforata) plotted on the opposite ends of PC1 axis and were distributed widely along PC2 axis (Fig. 5). Body size was not correlated notably or significantly with PC1 (r = 0.28, p = 0.10) and showed low correlation with PC2 scores (r =  − 0.36, p = 0.03).

Micro- and nano-scale analysis

A total of 36 nano-scale analyses (Data S4), with simultaneous nano-scale measurements of Mg/Ca mmol/mol ratios (Mg) and hardness (H) were obtained for four echinoid tests (nine analyses per specimen), including two specimens of Encope michelini (one from Northern Gulf and one from Florida Keys) and two specimens of Leodia sexiesperforata (both from Florida Keys), which were also used for bulk geochemical analyses.

When data were pooled across the two analyzed taxa, notable and statistically significant differences were observed in all comparisons between the nano-scale measurements for the outer and inner zones of echinoid tests (Fig. 8). The imperforate stereom of the outer zone was characterized by significantly elevated Mg/Ca and H values when compared to the thin galleried and coarse labyrinthic stereom of the inner zone: (1) Mg/Ca: median (outer) = 0.193, median (inner) = 0.168, W = 59.5, p = 0.005; (2) H: median (outer) = 4.54, median (inner) = 3.98, W = 15.5, p = 0.00002 (α_B = 0.05/6 = 0.008; two-sample Wilcoxon rank test). The median values of Mg/Ca and H were also higher for the outer zone compared to the inner zone when data were split by taxon (Table 6).

The comparison of two taxa (Table 6) for all data indicated that median Mg was significantly higher for L. sexiesperforata than for E. michelini (median (L. sexiesperforata) = 0.192, median (E michelini) = 0.158, W = 25, p = 0.0002, α_B = 0.05/6 = 0.008; two-sample Wilcoxon rank test). In contrast, median values of H were indistinguishable statistically between the two taxa: (1) H: median (L. sexiesperforata) = 4.25, median (E. michelini) = 4.03, W = 149, p = 0.69. The consistent differences in Mg between the two taxa persisted when per-taxon median values were computed for each zone separately, with higher median values observed for L. sexiesperforata when compared to E. michelini for both the inner and outer zones (Table 6).

Mg/Ca ratios and nanohardness H were positively correlated, with moderate-to-high correlation coefficient (Fig. 9; Table 7). For pooled data and for each taxon separately, these correlations were highly significant at Bonferroni-corrected significance level (α_B = 0.05/9 = 0.0055; Table 7). In the case of E. michelini, the values obtained for specimens from Northern Gulf appeared systematically offset toward lower Mg/Ca ratios, but nanohardness was comparable for the two specimens.

Univariate analyses (Table 6, Fig. 8) and correlation analyses (Table 7, Fig. 9) summarized above consistently indicated that Mg and H were elevated for the outer zone of the echinoid test when compared to the inner zone. Concurrently, there was a consistent offset in Mg/Ca ratio between the two taxa, with Mg values elevated for L. sexiesperforata when compared to E. michelini (Figs. 8–9). However, a joint comparison of the offsets between inner and outer zones for L. sexiesperforta, and E. michelini suggested that while the two taxa differed in Mg/Ca ratios, the nano-physical properties H were very similar (Fig. 10). The nanohardness H values of the inner and outer zones were nearly identical for the two taxa despite notable interspecific differences in Mg/Ca ratios, and the magnitude of the offset in H was also comparable between the two taxa (Fig. 8). Bootstrap-based confidence intervals around medians indicated substantial uncertainty (Figs. 8–10), but this uncertainty is not sufficient in magnitude to nullify those differences.