The extinction and survival of sharks across the end-Cretaceous mass extinction

Sharks (Selachimorpha) are iconic marine predators that have survived multiple mass extinctions over geologic time. Their fossil record is represented by an abundance of teeth, which traditionally formed the basis for reconstructing large-scale diversity changes among different selachimorph clades. By contrast, corresponding patterns in shark ecology, as measured through morphological disparity, have received comparatively limited analytical attention. Here, we use a geometric morphometric approach to comprehensively examine the dental morphology of 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 during the Cretaceous/Paleogene transition. Nevertheless, selective extinctions did impact on apex predator lineages characterized by triangular blade-like teeth, and in particular, lamniforms including the dominant Cretaceous anacoracids. Other groups, such as, triakid carcharhiniforms, squalids, and hexanchids, were seemingly unaffected. Finally, while some lamniform lineages experienced morphological depletion, others underwent a post-extinction disparity increase, especially odontaspidids, which are typified by narrow-cusped teeth adapted for feeding on fishes. This disparity shift coincides with the early Paleogene radiation of teleosts, a possible prey source, as well as the geographic relocation of shark disparity ‘hotspots’, perhaps indicating a regionally disjunct pattern of extinction recovery. Ultimately, our study reveals a complex morphological response to the end-Cretaceous mass extinction event, the dynamics of which we are only just beginning to understand.

such as, triakid carcharhiniforms, squalids, and hexanchids, were seemingly unaffected. 29 Finally, while some lamniform lineages experienced morphological depletion, others 30 underwent a post-extinction disparity increase, especially odontaspidids, which are typified 31 by narrow-cusped teeth adapted for feeding on fishes. This disparity shift coincides with the 32 early Paleogene radiation of teleosts, a possible prey source, as well as the geographic 33 relocation of shark disparity 'hotspots', perhaps indicating a regionally disjunct pattern of 34 extinction recovery. Ultimately, our study reveals a complex morphological response to the 35 end-Cretaceous mass extinction event, the dynamics of which we are only just beginning to 36 understand. 37

INTRODUCTION 38
Fossils provide the only direct evidence for the interplay between organisms and their 39 environments over vast evolutionary timescales [1][2][3]. They are, therefore, crucial for 40 exploring the drivers of biodiversity and ecosystem change in the past, with potential insights 41 on the origin of present-day ecosystems [3]. However, the analytical challenge is to discern a 42 genuine biological signal relative to geologic, taphonomic, sampling, and philosophical biases 43 [4][5][6][7]. While these may be impossible to overcome in entirety, the fossil records of some 44 widely distributed and chronostratigraphically extended clades offer reasonable proxies for 45 modelling macroevolutionary processes through deep time. 46 Sharks constitute just such a 'model group' because their dental remains are abundant in 47 Mesozoic and Cenozoic marine depositsa timeframe that covers ~250 million years (Ma) 48 [8,9]. Extant shark species are also ecologically ubiquitous, encompassing a spectrum of 49 macrophagous to microphagous predators that account for nearly half (42%) of all the 50 currently documented chondrichthyan biodiversity (n=1193 species) [10,11]. Nevertheless, 51 after nearly 200 years of scientific research [9], the various biological and environmental 52 factors that shaped shark evolution in the distant past remain obscure. In particular, their 53 capacity to survive past mass extinctions is relevant for understanding the dramatic decline of 54 shark populations observed in modern oceans [11][12][13][14][15]. 55 The end-Cretaceous mass extinction (~66 Ma), which marks Cretaceous/Paleogene (K/Pg) 56 chronostratigraphic boundary, is especially pertinent because it profoundly disrupted marine 57 ecosystems, but has had disputed implications for shark species diversity and morphological 58 disparity. Indeed, contrasting interpretations have advocated either limited [16], or complex 59 interrelationships between various biotic and abiotic drivers influencing shark evolution from 60 before, during, and after the K/Pg mass extinction event [17][18][19][20]. Here, we attempt to resolve 61 8 Fig) found that k = 160 best-captured tooth shape complexity (e.g., adequate shape coverage 134 of hexanchiforms) via distal and mesial curves of 78 and 79 sliding semilandmarks, 135 respectively. At this time, tooth serrations were not digitized due to image resolution, but we 136 acknowledge that they are functionally important [45]. 137 Finally, to screen for possible digitization errors, we extracted a random sub-sample of 30 138 images (S7 Table) and used a one-way ANOVA to calculate the intra-class correlation 139 coefficient (R) [25,46,47] (also see Supplementary text). We also ran a two-block partial 140 least-squares (2B-PLS) analysis to infer covariation between the datasets. Other exploratory 141 procedures performed as part of assessing measurement error include: 1) a manual survey of 142 digitized images to confirm accurate landmark placement and 2) screening for outliers. The 143 latter was performed visually, via observation of ordinated shape-space and associated thin-144 plate spline (TPS) deformation grids and analytically, through an outlier search function, 145 plotOutliers, in the geomorph package geomorph v. 3.3.1 [48]. 146

Morphometric analysis 147
To standardize digitized specimens to unit size, position, and rotation, we use a generalized 148 Procrustes analysis (GPA) [49,50]  where SSWn is the sum of the square distances between the coordinates (xk and yk) of 180 observation n, and their associated mean (̅̅̅ and ̅̅̅). This may be alternatively depicted as: 181 (2) All SSWn values in a given time-bin t are then summed and divided by the sample size at that 182 time (Nt) to measure PV, across all observations. 183 Disparity within each time-bin was partitioned according to our taxonomic order-level 184 classifications, which determined the clade-specific contributions to overall disparity. This 185 calculation equates to (2), but with distances measured relative to the mean of the group i 186 The partial PV for the group i over a specific time-bin t is  188 where ni is the sample size within group i and N is the total sample size within t. 189 Computationally, these equations are solved in geomorph [48], with the expectation that 190 additive partial disparities [59] for sampling within each time-bin approximate the total PV 191 given t. 192 11 We used a residual randomization permutation procedure (RRPP) with 1000 permutations 193 [55] to test null hypotheses for our multivariate shape data as: 194 where H0 assumes that pairwise absolute differences between PVs across two given time-bins 195 (e.g. t1 and t2) will be zero. HA alternatively stipulates that the difference will be greater than 196

zero. 197
We also applied non-parametric bootstrap resampling to estimate confidence intervals around bins to a minimum time-bin size (S1, S3 and S4 Table)  However, we also calculated partial disparities for each time-bin based first on depositional 209 basins, and then using countries of origin; designated i in (5). Although geopolitical 210 boundaries are artificial, they provide a convenient proxy for comparing regional versus 211 global disparity signals across a broader sub-sample series. difficult to define from isolated teeth, we conducted our analyses with the caveat that 225 adequate intraspecific coverage was assumed for each order-level clade. 226

Genus-level diversity 270
Raw genus-level occurrences showed a diversity increase in the Campanian to Maastrichtian 271 time-bins, followed by a prolonged decline throughout the Danian-Thanetian (Fig 3B). 272 However, sample-standardized curves using SQS and CR suggest stable diversity across the 273 extinction event (Fig 3B). 274

Global partial disparity 312
Our PV values indicate that selachimorph global disparity during the Campanian and 313 Maastrichtian was predominantly driven by lamniforms (Fig 6A). Not surprisingly, 314 lamniforms were also the most numerically abundant taxa in our dataset (S1 Fig). Lamniform 315 global disparity subsequently decreased in the Danian time-bin ( Fig 6A), but increased again 316 in the Selandian time-bin. Altogether, lamniforms accounted for almost the total global 317 variance observed within Selachimorpha (Fig 6A). Carcharhiniforms and hexanchiforms 318 otherwise contributed towards a global disparity increase across the Danian-Thanetian time-319 bins ( Fig 6A). 320 Although sampling is admittedly uneven, our break-down of family-level clade disparity (Fig.  321 6B) indicates that archaeolamnids, mitsukurinids, and anacoracids are the primary drivers 322 within the Campanian time-bin. Odontaspidid sharks increase their disparity following the 323 extinction of anacoracids in the Danian-Selandian time-bin, and expanded substantially in the 324 Thanetian time-bin (Fig. 6B). In correspondence, triakids, hexanchids, and squalids were the 325 major non-lamniform contributors to disparity in the Danian-Selandian time-bin, with 326 scyliorhinids becoming dominant in the Thanetian time-bin. Other family-level clades, such 327 as ginglymostomatids and hexanchids also increased their disparity from the Danian-328 Selandian to Thanetian time-bins (Fig. 6B). accounted for most of the disparity in our Campanian time-bin (Fig 7A and S14 Fig). Selandian time-bins (Fig 7C). Both the U.S.A. and Morocco sub-samples returned to high 343 disparity in the Thanetian time-bin (Fig 7D).  Pairwise comparisons along all PCs found significant differences in morphospace distribution 356 between time-bins (Fig 8, Table 2). The PC1 distributions are platykurtic (kurtosis < 3) and 357 positively skewed, except during the Thanetian, which was negatively skewed (g1Thanetian= -358 0.0336). Computed distribution-specific interquartile ranges (IQR) show minimal 359 morphospace dispersion on PC2 (S31-S32 Table). A Hartigan's dip test of multimodality did 360 not reject a unimodal distribution on either PC1 or PC2. However, a significant positive 361 morphospace shift occurred between the Maastrichtian-Danian along PC1 (Fig 8A, S16 Fig,  362 S31-S32 Table). Notably, there was no corresponding change in the modal-shape 363 configurations, although a reduction in positively loaded morphospace accompanied range-364 shortening among the minimum values ( Fig 8A, S16 Fig). The average value along PC1 shifted positively from the Selandian to the Thanetian (Fig 8A). 366

Superorder-level clade morphospace 393
Galeomorphs and squalomorphs occupied comparable regions on PC1 throughout the entire 394 Maastrichtian-Thanetian interval (Fig 10A), but exhibit notable differences in their mean and 395 medial values. On average, galeomorphs are characterised by tall, narrow teeth, whereas 396 squalomorphs are low crowned. The most pronounced shift occurs across the K/Pg Boundary 397 along PC1, as well as on PC3 and PC4 (S34 Table), where galeomorphs exhibit a reduction in 398 negative PC1 values (Fig 10A). Squalomorphs underwent a significant positive shift on PC2 399 (p=0.016) from the Maastrichtian to the Danian-Selandian (Fig 10A-B, S35 Table). 400

Disparity dynamics across the K/Pg Boundary 429
Our results indicate that the global disparity of selachimorphs across the K/Pg extinction 430 event was largely stable. Stasis is consistent with previous dental disparity-based studies of 431 sharks [20] and interpretations of the elasmobranch fossil record [16], but contrasts with 432 significant losses in species richness of almost 50% [17,18,77], suggesting shark diversity and 433 disparity were decoupled across the K/Pg Boundary. Importantly, our evidence for static 434 disparity in selachimorphs across the boundary is retained after accounting for heterodonty 435 and differences in order-level representation between the K/Pg time-bins (S15 Fig). Overall, 436 we find that these sources of variation only alter the disparity inferred for the Thanetian that, 437 upon correction, cannot be differentiated from that of the Campanian (see Fig. 3A, S15 Fig,  438 and S29 Table). In either case, morphological stasis about the K/Pg Boundary suggests that 439 the loss of 17% of families [17] was not linked to major losses in ecological diversity and, 440 overall, conforms to a 'non-selective' extinction model [78]. However, at lower taxonomic 441 scales, lamniforms experienced the selective extinction of Cretaceous anacoracids [20] with 442 triangular, blade-like teeth (Fig 11A-B) and, in the aftermath, lamniform morphospace 443 underwent a noticeable Paleocene proliferation of apicobasally tall, laterally cusped teeth, 444 driven by odontaspidids (Fig 11A-B). This dynamic of loss in one area of morphospace and 445 subsequent proliferation in another supports a 'shift' extinction model [20], the result of 446 which is a stable disparity. Concomitantly, carcharhiniforms characterized by mesiodistally 447 broad and low teeth proliferated during the Paleocene (Fig 11A-B), perhaps reflecting a 'non-448 adaptive' radiation along with minimal ecological divergence [20]. In comparison to 449 galeomorphs, squalomorph disparity seemed largely unaffected by the extinction event, 450 despite previously identified moderate taxonomic losses among squalids [17]. Much like the 451 overall pattern, squalids may have undergone a non-selective extinction. However, 452 22 squalomorphs remain poorly sampled across the end-Cretaceous extinction event (S1 Table), 453 and associated patterns should be interpreted cautiously. 454 Despite the proliferation of certain tooth types in the aftermath of the extinction event, static 455 disparity does not support a Paleocene diversification of sharks, which is in contrast to the 456 substantial ecological and taxic diversification of coeval actinopterygian fishes (e.g., from that recovered by our global disparity analysis. However, a slight disparity increase (Fig.  469 3A) was recovered from the early to middle Danian, which conforms to an overall increase 470 (albeit statistically non-significant) in absolute values amongst galeomorphs (Fig. 4A). This 471 regional pattern is indicative of likely short-term recovery in the aftermath of the extinction 472 event, but otherwise posit a direct correlation with distribution trends evidenced globally for 473 the early Paleocene. 474 On a broader level, regional disparity occurred in 'hotspots' that transitioned over the sampled 475 interval. Results indicate that (1) the Western Interior Basin assemblages from the U.S.A. 476 23 contributed to much of the measurable global disparity during the Campanian and 477 Maastrichtian (Fig 7A-B); (2) the Scandinavian Boreal assemblages from Denmark and 478 Sweden likewise predominate disparity for the K/Pg extinction interval (Fig 7C and S14 Fig);  479 and (3) the Mediterranean Tethyan and Atlantic margin basins of Morocco provide most of 480 the disparity for the Thanetian (Fig 7D). While these results are clearly influenced by 481 preservation biases (e.g., spatial variation in sampling across the Cretaceous/Paleogene

Ecological implications 493
The reorganization of selachimorph communities across the K/Pg mass extinction event 494 [17,18,20,24] can almost exclusively be attributed to large body-sized, apex-predatory 495 anacoracid lamniforms (e.g., the reduction of negatively loaded tooth structures) [20,21]. Regardless, the continuity of large-bodied sharks across the boundary, does suggest that 508 extinction susceptibility was not limited to relative body-size. Indeed, studies on living 509 cartilaginous fishes have shown that habitat preferences and reproductive strategies play a 510 more important role in determining survival, whereas, body-size has a modest effect on 511 extinction risk [96]. 512 Finally, for sharks, the K/Pg event should also be viewed in terms of diversification, rather 513 than simply extinction. In particular, lamniforms (odontaspidids) and carcharhiniforms 514 (triakids, scyliorhinids) morphologically proliferated in the early Paleocene, potentially taking 515 advantage of teleost fishes as an newly emergent food resource [16]. These patterns, highlight 516 the importance of biotic factors (i.e., predator-prey relationship) within the context of global 517 mass extinction, and it's contrasting effect among shark groups. 518

CONCLUSIONS 519
Understanding the evolutionary dynamics of sharks during the catastrophic end-Cretaceous 520 mass extinction event [17,18,20,21,77]  geometric morphometric evaluation of selachimorph disparity across multiple clades based on 524 their prolific dental fossil record. The combined observation of virtually static levels of 525 25 morphological disparity and the potential for differential extinction dynamics across the globe 526 suggest that, in a broad sense, the K/Pg extinction event was not as severe for sharks, as it was 527 for other vertebrates [90,91,99,100] Nevertheless, one particularly dominant Mesozoic clade 528 of sharks, anacoracid lamniforms, underwent a selective extinction, which saw to the loss of 529 triangular, blade-like teeth that were perhaps associated with feeding on large-bodied 530 contemporaneous tetrapods [20]. Other coeval selachimorph clades, however, survived and 531 appear of have been largely unaffected by the extinction event, with some groups of 532 carcharhiniforms and lamniforms even proliferating during the Paleocene. We interpret this as 533 an extinction-mediated ecological 'shift', whereby certain sharks underwent significant 534 changes in their morphology, but without substantial modifications to their overall disparity 535 [20]. The potential drivers of this transition include possible shifts in food resources, 536 supported by the Paleocene diversification of teleost fishes, although other mechanisms 537 cannot be rejected at this time. Finally, our application of a morphospace-disparity framework 538 offers a powerful interpretive tool that compliments traditional taxonomy-based assessments 539 of discrete character data [e.g., 101,102]. We therefore advocate similar continuously valued 540 disparity descriptors to reconstruct deep time evolutionary patterns in the fossil record of 541 sharks in the future, while acknowledging that variation sources, such as heterodonty and 542 differential spatiotemporal sampling, need to be accommodated. 543