(C) PLOS One [1]. This unaltered content originally appeared in journals.plosone.org.
Licensed under Creative Commons Attribution (CC BY) license.
url:
https://journals.plos.org/plosone/s/licenses-and-copyright
------------
Tooth morphology elucidates shark evolution across the end-Cretaceous mass extinction
['Mohamad Bazzi', 'Subdepartment Of Evolution', 'Development', 'Department Of Organismal Biology', 'Uppsala University', 'Uppsala', 'Nicolás E. Campione', 'Palaeoscience Research Centre', 'School Of Environmental', 'Rural Science']
Date: 2021-08
Sharks (Selachimorpha) are iconic marine predators that have survived multiple mass extinctions over geologic time. Their prolific fossil record is represented mainly by isolated shed teeth, which provide the basis for reconstructing deep time diversity changes affecting different selachimorph clades. By contrast, corresponding shifts in shark ecology, as measured through morphological disparity, have received comparatively limited analytical attention. Here, we use a geometric morphometric approach to comprehensively examine tooth morphologies in multiple shark lineages traversing the catastrophic end-Cretaceous mass extinction—this event terminated the Mesozoic Era 66 million years ago. Our results show that selachimorphs maintained virtually static levels of dental disparity in most of their constituent clades across the Cretaceous–Paleogene interval. Nevertheless, selective extinctions did impact apex predator species characterized by triangular blade-like teeth. This is particularly evident among lamniforms, which included the dominant Cretaceous anacoracids. Conversely, other groups, such as carcharhiniforms and orectolobiforms, experienced disparity modifications, while heterodontiforms, hexanchiforms, squaliforms, squatiniforms, and †synechodontiforms were not overtly affected. Finally, while some lamniform lineages disappeared, others underwent postextinction disparity increases, especially odontaspidids, which are typified by narrow-cusped teeth adapted for feeding on fishes. Notably, this increase coincides with the early Paleogene radiation of teleosts as a possible prey source, and the geographic relocation of disparity sampling “hotspots,” perhaps indicating a regionally disjunct extinction recovery. Ultimately, our study reveals a complex morphological response to the end-Cretaceous mass extinction and highlights an event that influenced the evolution of modern sharks.
Funding: This work was supported by the Royal Swedish Academy of Sciences (GS2017-0018) to M.B., and a Wallenberg Scholarship from the Knut and Alice Wallenberg Foundation to P.E.A. B.P.K. also acknowledges funding from a Swedish Research Council Project Grant (2020-3423), and N.E.C. is funded by an Australian Research Council Discovery Early Career Research Grant (DE190101423). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Copyright: © 2021 Bazzi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
The end-Cretaceous mass extinction (approximately 66 Ma), which marks the Cretaceous/Paleogene (K/Pg) chronostratigraphic boundary, is especially pertinent because it profoundly disrupted marine ecosystems but has disputed implications for shark species diversity and morphological disparity. Indeed, contrasting interpretations advocate either limited [ 17 ] or complex interrelationships of biotic and abiotic drivers that seemingly influenced shark evolution from before, during, and after the K/Pg mass extinction event [ 18 – 21 ]. Here, we explore these contentions via a comprehensive assessment of shark dental morphology across the end-Cretaceous mass extinction. Our study expands on previous studies that targeted either geographically localized [ 22 , 23 ] or clade-specific [ 21 ] assemblages. We use a dataset of 1,239 fossil shark teeth, representing 9 major selachimorph clades sampled at global and regional scales. These groups include the following: the Galeomorphii orders Carcharhiniformes, Heterodontiformes, Lamniformes, Orectolobiformes; the Squalomorphii orders Echinorhiniformes, Hexanchiformes, Squaliformes, Squatiniformes; and the extinct [†] Synechodontiformes. Our approach uses geometric morphometrics to compare disparity and morphospace patterns across a constrained 27.6-million-year interval spanning the Campanian and Maastrichtian ages of the Late Cretaceous (83.6 to 66 Ma) to the Danian, Selandian, and Thanetian ages (= Paleocene epoch) of the early Paleogene (66 to 56 Ma). We test the following hypotheses that: (1) selachimorph disparity was in decline before marine ecosystem disruption during the end-Cretaceous mass extinction, resulting from major marine transgressions during the Maastrichtian [ 23 , 24 ]; (2) selachimorph taxonomic richness depletion at the K/Pg boundary [ 18 , 19 ] was coupled with eco-morphological turnover; (3) apex predator shark lineages were disproportionately impacted, consistent with similar losses in marine tetrapods and osteichthyans [ 18 , 21 , 25 ]; and (4) extinction and recovery patterns were consistent at global and regional scales [ 21 – 23 ].
Sharks constitute one such group because their dental remains are abundant in Mesozoic and Cenozoic marine deposits—a timeframe covering approximately 250 million years (Ma) [ 8 , 9 ]. Extant shark species are also ecologically disparate, encompassing a spectrum of macrophagous to microphagous predators that account for nearly half (42%) of all the currently documented chondrichthyan biodiversity (N = 1,193 species) [ 10 , 11 ]. Nevertheless, the various biological and environmental factors that have shaped shark evolution remain obscure. In particular, their capacity to survive mass extinctions is relevant for understanding the dramatic decline of shark populations observed in our modern oceans [ 11 – 16 ].
Fossils provide the only direct evidence of interplay between organisms and their environments over vast evolutionary timescales [ 1 – 3 ]. They are, therefore, crucial for exploring the drivers of past biodiversity change and can shed light on the origins of modern ecosystems [ 3 ]. However, the analytical challenge is to discern a genuine biological signal from the combined obfuscations of geologic, taphonomic, sampling, taxonomic, analytical, and interpretive biases [ 4 – 7 ]. While these may be impossible to overcome in entirety, the fossil records of some widely distributed and chronostratigraphically extended clades provide exceptional opportunities to characterize macroevolutionary processes through deep time.
Because developmental [ 70 , 71 ] and ontogenetic factors [ 70 , 72 – 74 ] also affect selachimorph tooth morphology, our analyses are presented with the caveat that adequate intraspecific coverage was assumed for each order-level clade.
Sharks are known to exhibit both monognathic (variation along the tooth row) and dignathic (variation between the upper and lower jaws) heterodonty [ 8 , 68 ], although this can be difficult to discriminate from isolated fossil teeth [ 69 ]. To mitigate, we relied upon established diagnostic criteria [ 21 , 70 ] to categorize our specimens as representing parasymphyseal, anterior, lateroposterior, or posterior tooth positions (N MH = 897), and deriving from either the upper or lower tooth rows (N DH = 334) (see S1 Text ) (8) where the aligned Procrustes coordinates (shape) are described as a function of monognathic or dignathic Tooth Position at age t.
We accommodated for sample size biases inherent in the fossil record [ 4 , 62 – 64 ] via rarefaction comparisons of time-scaled PV [ 21 ]. These involved subsampling (999 iterations) of all time-bins to a minimum size (Tables A, C, and D in S1 Text ), after which 95% prediction intervals were calculated. Geographic subsampling focused on the UNESCO World Heritage fossil locality at Stevns Klint in Denmark (Fig H in S1 Text ), which preserves exceptionally rich selachimorph assemblages [ 65 , 66 ] spanning the K/Pg succession [ 67 ]. However, we also calculated partial disparities for each time-bin based on marine depositional basins, designated i in (6). Although the precise boundaries of these basins are ambiguous, they do provide a convenient proxy for comparing regional versus global disparity signals across a broader subsampled series.
We also applied nonparametric bootstrap resampling to estimate confidence intervals around disparity. All post hoc pairwise comparisons of group means were subject to false discovery rate (FDR) adjustments of p-values to mitigate the increased risk of Type I errors associated with multiple comparisons [ 61 ].
H 0 assumes that pairwise absolute differences between PVs across 2 given time-bins (e.g., t1 and t2) will be zero. H A alternatively stipulates that the difference will be greater than zero.
We used a residual randomization permutation procedure (RRPP) with 1,000 permutations [ 55 ] to test null hypotheses for our multivariate shape data as: (7)
Computationally, these equations are solved in geomorph [ 48 ], with the expectation that additive partial disparities [ 60 ] for sampling within each time-bin approximate the total PV given t.
The partial PV for the group i in a specific time-bin t is ∴ (6) where n i is the sample size of group i and N is the total sample size within t [ 48 ].
Disparity within each time-bin was partitioned according to their taxonomic order-level classifications, which determined the clade-specific contributions to the overall disparity. This calculation equates to (3), but with distances measured relative to the mean of the group i ( ), rather than the overall mean [ 48 ]: (5)
All SSW n values in a given time-bin t are then summed and divided by the sample size at that time (N t ) to measure PV across all observations [ 48 ]. Note that the following equation was erroneously presented in Bazzi and colleagues [ 21 ].
We calculated Procrustes variance (PV) as a measure of disparity [ 41 , 59 ] based on the 2D (k × m × N) Procrustes-aligned landmark dataset: (2) where SSW n is the sum of the square distances between the coordinates (x k and y k ) of observation n and their associated mean ( and ) [ 48 ]. This equation can be alternatively depicted as: (3)
Multivariate normality was assessed using a Henze–Zirkler test [ 58 ] (HZ = 4,956, p = 0) (Fig G in S1 Text ). Statistical comparisons between time-bins used a nonparametric Procrustes analysis of variance implemented in the RRPP package [ 55 ].
Morphospace was depicted as time-bin box plots that incorporate the arithmetic mean, median, and modal values. Confidence intervals were calculated using nonparametric bootstrapping with 1,000 resamples. Comparisons between multiple central tendency values accommodate for differences in frequency distributions. Modal shape configurations correspond to the region of maximum frequency calculated as: (1) where is the modal shape configuration of the jth PC; represents the mean shape configuration of the whole sample; is the mode of the jth PC axis; and, Γ j is the rotation matrix corresponding to the jth PC. and Γ j take the form of k × m matrices and were plotted as TPS deformation grids; “k” is the number of landmarks and “m” is the dimensions, in our case two.
The aligned Procrustes coordinates were ordinated via a principal components analysis (PCA) based on the singular value decomposition of the variance–covariance matrix. Shape variation was depicted as both TPS deformation grids and deformation isolines to generate a concentration “heat map” [ 52 ]. Morphospace was visualized using back-transformation [ 53 , 54 ]. All analyses were carried out in R v. 4.0.5 [ 44 ] with the geomorph v. 4.0.0 [ 48 ] and RRPP v. 1.0.0 [ 55 ] packages; visualization used the ggplot2 package [ 56 ]. All data and R scripts are available from Data Dryad under DOI:
https://doi.org/10.5061/dryad.c866t1g5n [ 57 ].
To standardize our digitized specimens for unit size, position, and rotation, we used a generalized Procrustes analysis (GPA) [ 49 , 50 ] that minimizes the bending energies to optimize the positions of the sliding semilandmarks [ 50 , 51 ]. Because large numbers of semilandmarks can impinge on GPA convergence [ 26 , 38 ], we varied the iteration frequency by arbitrarily increasing the max.iter argument in gpagan to compare convergence criteria Q-values (= Procrustes sum of squares). The resulting consensus shape configurations were then inspected (Fig F and Table H in S1 Text ).
Lastly, to screen for possible digitization errors, we extracted a random subsample of 30 tooth images (Table G in S1 Text ) and used a one-way ANOVA to calculate the intraclass correlation coefficient (R) [ 26 , 46 , 47 ] (also see S1 Text ). We also ran a 2-block partial least-squares (2B-PLS) analysis to infer covariation between the landmark datasets. Other exploratory procedures were employed to assess measurement error: (1) a manual survey of digitized images to confirm landmark placement accuracy; (2) screening of outliers visualized in the morphospace plots and associated thin plate spline (TPS) deformation grids, as well as being identified using the plotOutliers search function in geomorph v. 4.0.0 [ 48 ].
Landmark digitization was carried out in tpsDig2 v. 2.31 [ 43 ] with resampling to a standard number of equidistant semilandmarks using customized code in R v. 4.0.5 [ 44 ]. The resulting scheme comprised 2 open curves defined by semilandmarks but anchored by 3 fixed landmarks: Type 1 landmarks delimited the mesial and distal crown–root junctions, and a Type 2 landmark pinpointed the tooth apex (sensu [ 37 ]; Fig B in S1 Text , panel B, Table F in S1 Text ). The number of semilandmarks (k) was determined using resampling of the mesial and distal curves at equal spacings of k = 40, 60, 80, 100, 120, 140, and 160. Qualitative observations (Figs D and E in S1 Text ) found that k = 160 best-captured tooth shape complexity and included distal and mesial curves of 78 and 79 sliding semilandmarks, respectively. Tooth serrations were not digitized because of limited image resolution, but we acknowledge that these structures are functionally important [ 45 ].
Landmark-based geometric morphometrics quantifies biological shapes as a series of evolutionary homologous points in Cartesian space [ 26 , 27 , 37 – 42 ]. However, evolutionary homology [ 37 , 38 ] between landmarks cannot be assumed because of the inherent morphological variability in shark teeth (e.g., the location and number of cusplets). As a result, we designated our landmark placements based on topological rather than evolutionary homology [ 41 ].
Time-binning encompassed 5 geochronological ages spanning the immediate K/Pg interval: Campanian and Maastrichtian/Danian, Selandian, and Thanetian. However, we also implemented an alternative 4-age time-binning scheme that pooled temporally ambiguous specimens assigned to the Danian and Selandian (Tables C and D in S1 Text ). In addition, we carried out analyses using a subsample of the dataset (N = 659) for which “early” and “late” sub-ages could be defined for the Campanian and Maastrichtian and “early,” “middle,” and “late” for the Danian (Table E in S1 Text ). Numerical age values, in Ma, were taken from the International Chronostratigraphic Chart v2020/03 [ 36 ] and used as references for plotting.
Our dataset was compiled using photographs and graphic images derived from first-hand observations or the literature ( S1 Data , Fig A and Tables A and B in S1 Text ). Following recommended best practices [ 21 , 26 , 27 ], we screened the raw image data to include only those depicting complete tooth crowns with adequate resolution to determine the crown–root junction. The global sampling includes nearly all major selachimorph orders (Fig A in S1 Text , Fig B in S1 Text , panel A), except for Pristiophoriformes (Sawsharks), which are sparsely documented [ 8 , 28 ]. We also elevated Echinorhinidae to Echinorhiniformes [ 29 – 31 ] and employed a sensitivity analysis to test the morphospace occupation and disparity effects of †Synechodontiformes, which has been classified as either a clade within Galeomorphii or a neoselachian sister lineage (Fig B in S1 Text , panel A) [ 32 – 35 ]. Finally, most images were captured in the labial aspect, unless only a lingual view was available, and with the tooth apex directed to the left (Fig C in S1 Text ). The equivalency of labial and lingual views [ 21 ] was tested using ordinary least-squares linear models based on a subsample of the labial and lingual view data (N = 866).
Results
Digitization measurement error Visual comparison of the computed consensus (mean) tooth shapes (Fig I in S1 Text) indicates consistent digitization across landmark datasets. Accordingly, an intraclass correlation coefficient (R) of 2% (or 1 in 50) was calculated based on the aligned Procrustes coordinates and their error replicate counterparts (N = 30) (Table I in S1 Text). Pearson product–moment correlation (t = 1,874.7, df = 9598, p-value << 0.001, R = 0.99) and a 2B-PLS test (r-pls = 0.997, p-value = 0.001, Z = 7.136) unambiguously demonstrate dataset compatibility (Fig J in S1 Text).
Geographic distribution of disparity Samples from the North American Western Interior (WIB) and European epicontinental (EEB) basins accounted for most of the disparity during the Campanian (Fig 6). Alternatively, the Maastrichtian incorporated a substantial contribution from the Ouled Abdoun and Tarfaya basins (= Eastern Atlantic rim [EAR]) of Morocco (Fig 6). In general, geographic regions with small representative subsamples yielded comparatively low disparities during the Maastrichtian, including the North African Mediterranean Tethys and North American Western Atlantic rim (WAR) (Fig 6). By contrast, intensive sampling of Stevns Klint in Denmark and the Limhamn Quarry in Sweden skewed the Danian-Selandian geographic distribution toward these regions and revealed the inordinate influence of lagerstätten deposits on our global disparity signal (Fig 6). Lastly, the Thanetian showed a return to more geographically widespread sampling across a continuous trans-Atlantic belt, spanning the African epicontinental basins and EAR to the EEB and WAR (Fig 6). PPT PowerPoint slide
PowerPoint slide PNG larger image
larger image TIFF original image Download: Fig 6. Geographic distribution of disparity. Mesozoic marine depositional basins are based on Wretman and colleagues [75]. Basins with sample sizes of N < 10 were omitted. Pie graphs depict sample proportions. The data used in this analysis can be accessed online at
https://doi.org/10.5061/dryad.c866t1g5n.
https://doi.org/10.1371/journal.pbio.3001108.g006
[END]
[1] Url:
https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001108
(C) Plos One. "Accelerating the publication of peer-reviewed science."
Licensed under Creative Commons Attribution (CC BY 4.0)
URL:
https://creativecommons.org/licenses/by/4.0/
via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/