The quality of the fossil record across higher taxa: compositional fidelity of phyla and classes in benthic marine associations

Materials

To evaluate fidelity across multiple higher taxa simultaneously, concurrent live and dead samples were collected in a series of dredges over four field seasons (June 2011, November 2011, May of 2012, and April 2013) in Onslow Bay, North Carolina (U.S.A.). Dredging was conducted at 52 localities (Fig. S1), 16 of which were sampled in multiple years, resulting in a total of 220 samples collected from a variety of habitats, depths, and distances from shore (see Tyler & Kowalewski, 2018 and Supplemental Appendix for additional sampling details). Samples included invertebrates from six phyla (Annelida, Arthropoda, Brachiopoda, Cnidaria, Echinodermata, and Mollusca). All live macro-invertebrates that could be seen without the aid of a microscope were counted and identified to the lowest taxonomic level (typically species). Encrusting species (e.g., bryozoans, barnacles, sponges) were excluded. Concurrent with live surveys, two bulk samples of dead material were collected at each locality (7.5 liters) and wet-sieved (mesh size 4.76 mm). As many organisms have multiple skeletal elements (e.g., valves, plates, spines) that may disarticulate after death, each element was corrected by dividing element counts by the number of elements estimated to be contained within a single live individual (Tyler & Kowalewski, 2018). Skeletal fragments were excluded to prevent distortion of the DA composition which can result from inflating abundance estimates of species easily identifiable from small fragments (Kowalewski et al., 2003). Abundance data are available at https://github.com/tylercl/Multi-Taxic-Fidelity (DOI: 10.5281/zenodo.7871639).

Compositional and diversity fidelity

Prior to all analyses small samples with fewer than 50 individuals were excluded (36 out of 52 sites were retained in the final analysis). As samples from DAs typically contain more individuals, sample standardization was performed using rarefaction, and samples were rarefied to the smallest of the two assemblages in any given live-dead pair. The fidelity of richness across and within phyla, classes, and species was estimated using Spearman’s and Pearson’s measures of correlations. For phyla and classes, taxon-focused estimates were derived for data pooled across all localities with live-dead pairs of points representing one regional estimate of species richness of a given higher taxon. For locality-focused estimates, species richness was estimated separately for each higher taxon and fidelity was measured as agreement between live and dead species richness within a given higher taxon across localities. Accelerated bias-corrected methods were used here to obtain improved estimates of confidence intervals obtained from bootstrap-derived distributions (e.g., Chernick & LaBudde, 2011). Thus, to assess fidelity within phyla and classes, 95% confidence intervals were estimated using bootstrapping with an accelerated bias correction (1,000 iterations). Pearson’s correlations of richness, Spearman’s rank correlations of proportional rank abundance (calculated using the PaleoFidelity R package Kowalewski, 2023), and Bray-Curtis similarity were calculated using sample standardized data. The three measures of correlation/similarity tend to co-vary with one another but Pearson’s and Bray Curtis provide estimates weighted toward abundant taxa and are more sensitive to outliers. Spearman’s rank correlation weights all taxa equally regardless of their abundance.

To determine to what extent DAs provide meaningful estimates of diversity across multiple higher taxa simultaneously, the following diversity estimates were calculated (two sites with fewer than 20 individuals were removed to retain meaningful numbers of individuals across phyla): richness (S), Shannon’s H, Simpson’s D, and Pielou’s evenness (J). To assess the reliability of DAs in capturing diversity, three assemblages were compared: (1) the multi-taxic LA and the molluscan DA; (2) the multi-taxic LA and non-mollusc DA; and (3) the mollusc DA and the non-mollusc DA. Mean diversity between assemblage pairs was compared using a Kruskal-Wallis test with a Bonferroni correction using the rstatix package (Kassambara, 2021).

To assess whether compositional fidelity differed across localities due to skeletal robustness (i.e., preservation biases), data were divided into two groups: molluscs, which have more durable skeletal components, and non-molluscs, with a larger proportion of species that are either soft-bodied or have less robust skeletal elements. This coarse grouping allows adequate retention of within-group sample sizes within localities. Similarity and rank-order agreement were then re-calculated for each subset separately using rarefied richness. Two sites with fewer than 20 individuals were excluded to retain meaningful numbers of samples and sites across phyla, leaving 50 sites.

To further assess compositional fidelity, the expected ‘perfect’ fidelity was estimated separately for each live-dead comparison. Fidelity estimates are highly sensitive to small sample sizes and unbalanced sampling. Unbalanced sampling is common when sampling the marine benthos (including this analysis) because live specimens tend to be far more rare relative to dead remains. It is therefore instructive to compare observed values of statistics of interest (e.g., Spearman’s rank correlations) to those expected if the fidelity were perfect and departures from perfect concordance (e.g., Spearman’s Rho = 1) was due to sampling alone. The estimates were obtained by randomization via reassignment of group membership by random resampling without replacement (e.g., Kowalewski & Novack-Gottshall, 2010; Manly & Alberto, 2021). For a given live-dead comparison both live and dead specimens were pooled together to create a single species distribution and specimens were then assigned randomly to live and dead samples using their original sample sizes. Because this randomization samples a single pooled distribution, the fidelity should be perfect, and any departures (e.g., Spearman’s Rho < 1) would reflect (and estimate) sampling effects (see also Kowalewski, 2023; PaleoFidelity R Package). Each randomization was based on 1,000 replicate samples.

Differences in multivariate dispersions between live and death assemblages were examined using a modified analysis of Homogeneity in Multivariate Dispersions (HMD; Tomasovych & Kidwell, 2011), which can be used to quantify the effects of premortem and postmortem variation. Analyses were conducted separately for the multi-taxic, mollusc, and non-mollusc assemblages. Two common dissimilarity measures were calculated, Jaccard dissimilarity which relies on presence/absence data, and Bray-Curtis, which utilizes species abundance data. Both measures were calculated using rarefied standardized data (six sites with fewer than 30 individuals were excluded, leaving 46 sites). Abundance was standardized by dividing the abundance of a given species by the total locality abundance. HMD was performed using a custom function in R (R Core Team, 2019) developed by Tomasovych & Kidwell (2011).

Diversity fidelity

At the coarsest, regional scale, a single live-data comparison can be carried out for higher taxa by pooling data across all localities. In terms of agreement in relative abundance of phyla, mollusc specimens dominated both the LA and DA (Figs. 1A–1B). With the exception of arthropods, rank-order abundance in the DA remained consistent with the LA. The live-dead agreement in abundance was lower at the class level, mostly reflecting discordances in bivalves (most abundant in the DA but ranked fourth in the LA) and arthropods (most abundant in the LA but ranked sixth in the LA). Nevertheless, four out of the five classes that were most abundant in the DA were also among the top five in the LA. Moreover, gastropods and polychaetes were ranked second and third in terms of specimen abundance in both the DA and LA samples (Figs. 1E–1F).

In terms of agreement in species diversity, molluscs were the most diverse phylum in both the LA and DA assemblages (Figs. 1C–1D; Table 1). The non-standardized species richness of phyla was positively correlated between the LA and DA assemblages (Spearman’s Rho = 1, p = 0.003; Pearson’s r = 0.91, p = 0.01; Fig. 2A; Table S1), with molluscs dominating both assemblages. In fact, the relationship was monotonic with a perfect live-dead agreement in rank abundance of the five phyla. The live-dead agreement in non-standardized species richness was also high when comparing classes, with the diversity of classes in the DA being a good predictor of diversity of classes in the LA (Pearson’s r = 0.79, p = 0.004; Spearman’s Rho = 0.79, p = 0.004; Fig. 2C; Table S1). As with molluscs, the bivalves and gastropods were the most diverse classes in the DA.

Phyla with high abundance in the LA also had high richness (Pearson’s r = 0.96, p = 0.003; Spearman’s Rho = 0.96, p = 0.0005), as did classes (Pearson’s r = 0.97, p < 0.001; Spearman’s Rho = 0.94, p < 0.001). This relationship was similarly strong in the DA for both phyla (Pearson’s r = 0.99, p < 0.001; Spearman’s Rho = 0.75, p = 0.05) and classes (Pearson’s r = 0.93, p < 0.001; Spearman’s Rho = 0.95, p < 0.001).

At the locality level, the correlation between the rarefied number of species in the multi-taxic LA and DA was strong and significant (Rho = 0.71, p < 0.001, 50 sites; Figs. S2A–S2B). Similarly, when assemblages were constrained to molluscs, locality level richness in mollusc LAs was strongly correlated with richness in sympatric mollusc DAs (Rho = 0.56, p < 0.001, 36 sites; Figs. S2C–S2D). The relationship was weaker and non-significant for non-molluscs (Rho = 0.02, p = 0.21, 17 sites; Figs. S2E–S2F). However, this could reflect the small number of samples in the non-mollusc assemblage (17 sites). The mollusc DAs were less informative as a proxy for site-level richness in the multi-taxic LA (r = 0.57, p = 0.02, 36 sites; Figs. S2G–S2H). When locality-level estimates of sample-standardized alpha diversity were obtained separately for each phylum (Fig. 2B), the correlation estimates were highest for molluscs and arthropods, but slightly lower than the multi-taxic correlation estimates reported above (Fig. 2B). The live-dead correlation was low for all other phyla (Fig. 2B). Similar results were obtained for locality-level sample standardized estimates of alpha diversity for classes, with highest positive correlations obtained for bivalves and echinoids (Fig. 2C). Non-sample standardized locality-level estimates of alpha diversity within phyla and classes (Tables 1 and 2) resulted in similarly strong correlations between live and dead molluscs and bivalves, but also produced elevated correlation values for arthropods, crustaceans, and echinoids. Results were comparable regardless of which correlation method was used (Spearman versus Pearson).

Compositional Fidelity

At the regional scale with data pooled across all localities, compositional agreement for multi-taxi data is low (Fig. 3; Fig. S3A) regardless of the measure used (Spearman’s Rho = −0.01; Pearson’s r = 0.002; Bray-Curtis similarity = 0.09). The fidelity was slightly higher for molluscs (Fig. 3; Fig. S3A) when using rank correlation (Spearman’s Rho = 0.29) but comparably low using other measures (Pearson’s r = 0.02; Bray-Curtis similarity = 0.08). Similarly to multi-taxic data, the compositional fidelity is low for non-molluscs (Fig. 3; Fig. S3A) with rank correlation Spearman’s Rho of −0.002. However, other measures are somewhat elevated for non-molluscs (Pearson’s r = 0.31; Bray-Curtis similarity = 0.10). Nevertheless, regardless of grouping and measures used, the LA and DA samples are mostly or entirely discordant in terms of their species composition. Similar patterns are observed when compositional fidelity is analyzed at the locality level, with Spearman’s rank correlation (Fig. 3; Fig. S3A) notably higher for molluscs when compared to non-mollusc and multi-taxic data.

Because measures of correlation and similarity can be affected by sample size and unequal sample sizes, and both those issues affect the data analyzed here, it is useful to assess whether sampling biases can explain the observed low compositional fidelity. However, the resampling model of perfect fidelity indicates that the sampling effects are relatively minor and do not vary notably across datasets (Fig. 3). That is, if fidelity were perfect, the estimates of Spearman’s Rho for regional data are expected to be around 0.9 for all three datasets (Fig. 3) and cannot account for the very low Spearman’s Rho observed for all datasets or differences in observed Rho values between molluscs and the other two datasets. Similarly, at the locality level, the expected Rho values for perfect fidelity range between 0.6 and 0.8 and cannot account for differences between mollusks and non-molluscs or low correlation values.

At the locality level, tests for Homogeneity in Multivariate Dispersions (HMD) indicate that DAs do not occupy the same region of the multivariate space as LAs, and are thus overdispersed relative to the centroid of LAs (Fig. 4). This pattern is observed consistently for the multi-taxic assemblage, mollusc and non-mollusc assemblages. However, multivariate dispersion was less pronounced in mollusc assemblages relative to the multi-taxic and non-mollusc assemblages when using Jaccard dissimilarity. All three DAs were significantly over-dispersed, however, the magnitude of overdispersion was somewhat low (0.15−0.27), and total live-dead variation was larger than premortem variation using either presence-absence or proportional abundance (Table S2). The proportion of total variation explained by premortem variation was similar in all three assemblages using both presence-absence (Jaccard) and standardized relative abundance (Bray-Curtis; Table S2). However, the magnitude of the total live-dead variation explained by the premortem component was marginally larger for mollusc assemblages relative to the multi-taxic assemblage when using Jaccard dissimilarity (0.63 vs. 0.57) and smaller when using Bray-Curtis (0.45 vs. 0.57; Table S2). Conversely, the non-mollusc assemblage was lower using Jaccard (0.47) and marginally higher using Bray-Curtis (0.61).