# Morphological response accompanying size reduction of belemmites during an Early Jurassic hyperthermal event modulated by life history

Paulina S. Nätscher, Guillaume Dera, Carl J. Reddin, Patrícia Rita, Kenneth de Baets

► To cite this version:

Paulina S. Nätscher, Guillaume Dera, Carl J. Reddin, Patrícia Rita, Kenneth de Baets. Morphological response accompanying size reduction of belemmites during an Early Jurassic hyperthermal event modulated by life history. Scientific Reports, 2021, 11, 10.1038/s41598-021-93850-0. insu-03661478

HAL Id: insu-03661478

https://insu.hal.science/insu-03661478

Submitted on 7 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Morphological response accompanying size reduction of belemnites during an Early Jurassic hyperthermal event modulated by life history

Paulina S. Nätscher<sup>1,✉</sup>, Guillaume Dera<sup>1,✉</sup>, Carl J. Reddin<sup>1,✉</sup>, Patrícia Rita<sup>1,✉</sup> & Kenneth De Baets<sup>1</sup>

One of the most common responses of marine ectotherms to rapid warming is a reduction in body size, but the underlying reasons are unclear. Body size reductions have been documented alongside rapid warming events in the fossil record, such as across the Pliensbachian-Toarcian boundary (PToB) event (~183 Mya). As individuals grow, parallel changes in morphology can indicate details of their ecological response to environmental crises, such as changes in resource acquisition, which may anticipate future climate impacts. Here we show that the morphological growth of a marine predator belemnite species (extinct coleoid cephalopods) changed significantly over the PToB warming event. Increasing robustness at different ontogenetic stages likely results from indirect consequences of warming, like resource scarcity or hypercalcification, pointing toward varying ecological tolerances among species. The results of this study stress the importance of taking life history into account as well as phylogeny when studying impacts of environmental stressors on marine organisms.

The most common responses of recent and fossil marine invertebrates to global warming are extinctions<fnref n="1" />, changes in geographical distribution², phenology³,⁴, and decreases in body size⁵⁻¹⁰. However, many biological mechanisms may underlie body size reductions⁸,¹⁰⁻¹³. Generally, individuals in colder environments grow more slowly but become larger as adults¹⁴. This widespread pattern is embodied by two well-established hypotheses: Bergmann's rule, which describes the negative association between temperature and body size in natural environments¹⁵, and the temperature-size rule (TSR)¹⁶, which describes physiological reaction norms relating temperature to body size in laboratory experiments⁸,¹¹,¹⁷,¹⁸. These hypotheses predict that organisms should grow to be larger in colder environments, because growth efficiency decreases with increasing environmental temperature¹⁴. Whether through metabolic processes¹⁹,²⁰ or reduced resource availability, exposure to higher temperatures and accompanying factors like increased *pCO₂* during early development have been shown to lead to altered growth processes²¹⁻²³. Possible outcomes include abnormal growth during the juvenile life stage, or precocious maturity¹³,²⁴, which might lead to paedomorphosis, a heterochronic change marked by retention of juvenile features in adult organisms⁵. While a decrease in body size has been reported from various past hyperthermal events²⁵⁻³³, morphological changes of different life stages are not commonly studied³⁴ although they are crucial for understanding long-term impacts on life history and their ecological implications.

Here we tackle this issue by focusing on one of the most severe extinction crises of the Jurassic: the Pliensbachian-Toarcian crisis. This event was a global, multi-phased, hyperthermal event triggered by volcanic episodes in the Karoo-Ferrar Large Igneous Province³⁵⁻³⁸, which led to the extinction of around 15–20% of marine species³⁹. This disturbance was characterized by an early abrupt warming of at least 5° C at the Pliensbachian-Toarcian boundary event (PToB)⁴⁰⁻⁴³, which is considered a precursor crisis to the early Toarcian warming and oceanic anoxic event (T-OAE)⁴⁴,⁴⁵.

<footnote n="1">¹Geozentrum Nordbayern, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany. ²GET, Université Paul Sabatier, CNRS UMR 5563, IRD, Toulouse, France. ³Museum für Naturkunde, Leibniz Institute for Evolution and Biodiversity Science, Berlin, Germany. ⁴MARE (Marine and Environmental Sciences Centre), 3004-517 Coimbra, Portugal. ✉email: paulina.naetscher@fau.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Palaeoenvironmental proxies in the Peniche GSSP section fluctuate across two hyperthermal pulses during the Pliensbachian-Toarcian crisis (PToB and T-OAE). The proxies from left to right are: 2nd order transgressive-regressive cycles to visualise changing bathymetry⁸⁷, δ¹³C as a proxy for carbon cycle disturbances⁴², δ¹⁸O as palaeotemperature proxy⁴², and δ¹¹B as a proxy for ocean acidification⁴². The grey lines in the background of the isotope curves show the original values⁴², whereas the thicker lines overlaying them represent the lower-resolution curves with data points only at the positions of the stratigraphic beds sampled in this study (see stratigraphic log on the left). The second-to-right curves show median GPA Centroid size as a body size proxy (grey line shows full assemblage pattern, with shaded area representing median absolute deviation (MAD), C. *longiforma* and *P. bisulcata* are shown in green and yellow, respectively, with error bars representing MAD). On the right-hand-side PC1 is displayed as proxy for robustness in purple (for colouration see Centroid Size description). The two belemnite shapes bracketing the robustness curves show typical cross section appearances for C. *longiforma* (in green on the left) and *P. bisulcata* (in yellow on the right).</figcaption></figure>

The Peniche GSSP section (Portugal) offers an excellent geochemical and fossil record of the PToB crisis (Fig. 1). Trace fossils are present throughout the section, indicating consistently oxygenated bottom waters⁴⁶. Sedimentological and geochemical proxies show the rapid warming event⁴⁰⁻⁴³ and sea-level transgression at the PToB⁴⁵ with evidence for ocean acidification⁴², increased weathering⁴¹, collapsed productivity and resource scarcity⁴⁰,⁴⁷, which was followed in the early Toarcian by a second, more severe environmental crisis (the T-OAE) (Fig. 1). The temporary disappearance of macroinvertebrates like belemnites (Fig. 1), additionally indicates harsh conditions during the T-OAE⁴⁸. Belemnites are among the organisms that suffered considerably during this crisis⁴⁸⁻⁵². Their abundance in the fossil record⁴⁸,⁵³⁻⁵⁵, and rapid response to environmental changes, makes these predators suitable for studying responses to past climate changes.

A significant size decrease in adults of the dominant belemnite species was evidenced across the PToB in Peniche⁴⁸, while the extinction of forms with slender belemnite rostra but consistent diversity of robust forms is interpreted as a selective extinction against deep-water dwelling slender forms which were sensitive to decreased oxygen concentrations⁵¹.

Here we use 3D landmark analysis⁵⁶ to quantify morphological changes of belemnites across the PToB crisis in Peniche. Basing our hypotheses on prior studies, we investigate, whether the size reduction across the boundary was accompanied by allometric changes and thus affected certain life stages disproportionally. Firstly, we expect the largest shift in morphospace (see results and methods) across the PToB as an increase in robustness⁴⁸,⁵¹. If this was driven by accelerated growth, we expect a within-species morphological change driven by paedomorphosis of the best represented species⁴⁸. We test for morphological changes on the (1) assemblage and (2) species level and (3) by life stage. Furthermore, we expect these changes to be correlated with increased temperature, rather than with severe shifts in depositional environments.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Subspace occupation within the two-dimensional morphospace (PC1, PC2) of the taxa present in this study. Points represent specimens and are coloured by taxon. The convex hulls connect the outermost values on PC1 and PC2 of each individual taxon. The black shapes around the figure represent the maximum and minimum shapes of the sagittal section through the belemnite rostrum of PC1 and PC2. (see sample sizes in Table 1)</figcaption></figure>

## Results

**Belemnite morphospace structuration.** The consensus shape produced through the Generalized Procrustes Analysis (GPA) is a ventral to dorsal profile through an almost cylindrical belemnite rostrum of medium robustness. The alveolar walls at the reference point on transversal sections are of average thickness and the protoconch is located close to the centre of the rostrum with almost equal distance to the dorsal and ventral sides⁵⁷. The first axis of the Principal Components Analysis (PCA) represents the ratio of width to length (robustness) of the cross-section and accounts for 82% of morphological variation in the data. PC2 explains 7% of the variance and the morphological variation is only nuanced, contrasting thin from thick alveolar walls in the transversal section representing another aspect of rostrum robustness. Based on the percentage of explained variance, we use PC1 as our main morphology proxy, and use 'PC1' and 'robustness' interchangeably.

The studied species occupy different regions of the morphospace (Fig. 2, supplementary table S1). More slender species include *Bairstowius amaliae*, *Catateuthis longiforma*, *Acrocoelites* sp. and indeterminable *Hastitidae* specimens. Out of these, the robustness of *C. longiforma* and *B. amaliae* differs significantly (*p* = 0.014). *Passaloteuthis milleri* and *Passaloteuthis bisulcata* show a wide range in terms of robustness from moderately slender to very robust shapes. Significant differences between both of these more robust species and *B. amaliae* are found (*P. milleri*, *p* < 0.001, *P. bisulcata*, *p* < 0.001). Additionally, the robustness of *C. longiforma* is significantly lower than both, *P. milleri* and *P. bisulcata*, too (*P. milleri*, *p* = 0.003, *P. bisulcata*, *p* < 0.001). *Parapassaloteuthis aff.* *zieteni* is, on average, the most robust species, differing significantly from all other species (*P. bisulcata*, *p* = 0.006, *P. milleri*, *p* < 0.001, *B. amaliae*, *p* < 0.001, *C. longiforma*, *p* < 0.001) (supplementary table S1).

**Morphology through time.** The full assemblage shows both a small but significant size decrease across the Pliensbachian-Toarcian boundary (Hedges' *g* = -0.022, *p* = 0.018), as well as a large, significant size increase from the *Dactylioceras semicelatum* to the *Elegantuliceras elegantulum* subzone after the T-OAE (Hedges' *g* = -6.302, *p* = 0.006), in which only two specimens of the species *Acrocoelites* sp. appear after a comprehensive loss of belemnite abundance and diversity⁴⁸,⁵⁸. The body size pattern of the full assemblage is not accompanied by any significant changes in robustness (Fig. 1, supplementary table S3).

The bootstrapping procedure (see methods) shows that the change in species abundances through the subzones is within the realm of random sampling, with the only exception being the *Elegantulum* subzone (supplementary Fig. S2). Because this sample is from after the T-OAE and the changes in community composition across the boundary appear to be random, we reject the hypothesis that community composition markedly influenced the morphological pattern across the PToB.

For species-level analyses, we focus on *P. bisulcata* and *C. longiforma*, which dominate the assemblages across the boundary and are reasonably represented in our morphological subset. *Catateuthis longiforma*, and to a lesser extent *P. bisulcata* have been shown to decrease in body size in response to the PToB event⁴⁸. While the sample of
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Linear regression lines of GPA Centroid Size with PC1 for cooler (blue line; empty points) and warmer (red line; filled points) times in *C. longiforma* (a) and *P. bisulcata* (b). The shade of the points indicates ontogenetic stage (light grey: juvenile, medium grey: neanic, black: adult). Significance levels of the regressions are indicated next to the legend.</figcaption></figure>

*P. bisulcata* in our dataset does not change in size (supplementary table S2), *C. longiforma* shrinks significantly across the boundary (*p* < 0.001) and increases in body size again in the following subzone (*p* = 0.015) (Fig. 1, supplementary table S2). As assumed (see methods), both species become more robust across the boundary (*P. bisulcata*: Hedge's *g* = -1.054, *p* = 0.043, *C. longiforma*: Hedge's *g* = -0.844, *p* = 0.02) (Fig. 1, supplementary table S3). After this initial increase in robustness, *C. longiforma* returns to a slenderer shape in the *Dactylioceras semicelatum* subzone (Hedges' *g* = 0.229), while *P. bisulcata* becomes even more robust (Hedges' *g* = -1.165) (Fig. 1, supplementary table S3). However, these changes from the *Dactylioceras mirabile* subzone to the *Semicelatum* subzone are not significant (supplementary table S3).

In both species, smaller specimens are generally more robust (Fig. 3). The slope of robustness at size becomes steeper in higher temperatures in both species (*P. bisulcata*—warm: slope = -0.003, adj. R² = 0.778, *p* < 0.001, *P. bisulcata*—cold: slope = -0.001, adj. R² = -0.047, *p* = 0.954, *C. longiforma*—warm: slope = -0.004, adj. R² = 0.309, *p* = 0.004, *C. longiforma*—cold: slope = -0.002, adj. R² = 0.458, *p* < 0.001), meaning that the differences in
robustness among life stages are less noticeable in lower temperature settings, while they become more pro-
nounced in high temperature habitats (Fig. 3). The difference in growth trajectories between warm and cold
settings is insignificant in *C. longiforma* (*p* = 0.138), but significant for *P. bisulcata* (*p* < 0.001).

The linear regression model produced significant evidence for developmental stage being a good predictor of morphological change in *C. longiforma* (*p* = 0.024) (supplementary table S4). *Catateuthis longiforma* shows a negative allometric growth during both Late Pliensbachian subzones, defined as background times, with older life stages being more slender than juveniles (Fig. 4). The older life stages (neanics and adults) become more robust across the boundary (Hedges' g = -0.904, *p* = 0.054) (supplementary table S5), making it look like the growth pattern switches to isometric growth during the hyperthermal event (Fig. 4). In *P. bisulcata* adults basically disappear after the boundary and juveniles become significantly more robust over time (*p* = 0.005) (Fig. 4). However, this seems to be happening gradually, as differences between consecutive subzones are insignificant (supplementary table S5).

The independent effects of temperature reflected by oxygen isotopes (p = 0.002), seawater pH, reflected by
boron isotopes (p = 0.028), and lithology (p = 0.001) on overall morphological variation are significant, out of
which lithology plays the biggest role with 1.4% (supplementary Fig. S4, supplementary table S6). However, this
cannot explain any changes across the boundary as the lithology variable in our samples stays constant there.
When correcting for the influence of lithology to identify correlations between robustness and environmental
variables, the null model is the most parsimonious for the full assemblage (supplementary table S7). Therefore,
no environmental factor seems to consistently influence the overall robustness pattern. For *P. bisulcata* the model
including boron isotopes as proxy for acidification is the most parsimonious (supplementary table S7) and per-
forms significantly better than the null model (p = 0.014) (supplementary table S8). In *C. longiforma* all variables
together, δ¹⁸O, δ¹³C and δ¹¹B, best explain the robustness pattern through time (supplementary table S7). In this
model, all environmental variables have a significant correlation with robustness (δ¹⁸O: p = 0.001, δ¹³C: p = 0.01,
δ¹¹B: p = 0.002) (supplementary Table S8). These results do not change drastically, when using palaeoenviron-
mental proxies (δ¹⁸O, δ¹³C, and Hg/TOC) from different publications³⁸,⁴³,⁵⁹, except for the significant correlation
of the morphological pattern of *P. bisulcata* with the acidification proxy (δ¹¹B), which is not present in the older
geochemical publications (supplementary tables S9 & S10).

Discussion

Belemnites were likely dominant in epipelagic habitats (0–150 m)⁵⁵ and fast-swimming nektonic predators⁶⁰,
and consequently must have had a demanding metabolism⁶¹. Even though we find no significant change in
robustness of the full assemblage across the PToB, the two dominant species (*C. longiforma* and *P. bisulcata*)
both show a significant increase in robustness at the event, which can be attributed to direct and indirect conse-
quences of the warming of the environment, like changes in continental run-off, ocean acidification and resource
sparsity⁴⁰–⁴²,⁴⁶,⁶². At the onset of the T-OAE, a shift towards more robustly shaped taxa (often interpreted as more
nektobenthic, shallow water inhabitants) was interpreted as an effect of anoxia in deeper-water habitats, inhabited
by more slender, pelagic belemnite taxa⁵¹,⁵⁵,⁶¹,⁶³,⁶⁴. However, in our material, the random change in relative spe-
cies abundances across the boundary would be unlikely to impact the increase in assemblage-scale robustness,
which leads us to discard the possibility of bathymetric changes or taphonomic processes having a large impact
on the morphological signal. If these processes had been substantial, they would have resulted in preferential
preservation of more robust species⁶⁵,⁶⁶, or an environmentally triggered shift of assemblage composition, which
is not the case in our data across the PToB.

We interpret the change in morphology of adult specimens that occurs parallel to the decrease in size in *C. longiforma* as an indication for a paedomorphosis by progenesis, i.e. the retention of juvenile morphology in stunted adult specimens<sup>5,34,67–70</sup>. A body size decrease in adults of *C. longiforma* has been observed in the full data set of belemnites in Peniche before<sup>48</sup>. This conclusion is further supported by a correlation of the morphological pattern with environmental changes in palaeotemperature, seawater pH and carbon cycle and the fast recovery after the crisis. In this species juveniles show a lower sensitivity to climate warming than adults. Differential sensitivities of life stages to climate-related stressors (CRS) have been found in many recent metazoan groups<sup>71–74</sup>. It has been argued that particularly very early life stages and spawning adults of marine ectotherms might be particularly susceptible to CRS<sup>8,71,75</sup>, while larger juvenile and young adults should be more resistant. If juvenile growth and reproduction is hindered, this likely will impact their population success<sup>24</sup>, which, if happening over larger parts of their geographic range, augments their extinction risk<sup>76,77</sup>. *Catateuthis longiforma* appeared to respond to the warming event by developmental changes resulting in smaller and more robust adults, and then recovered in morphology and population structure following the initial warming event. This suggests its changes in life history strategy allowed this species to more successfully cope with the initial warming event, which is reminiscent of observations on modern ectotherms<sup>11,13</sup>.

Although the sample sizes of the different *P. bisulcata* ontogenetic stages are too small for sound ecological
inferences, juveniles of this species sampled after the PToB appeared to show abnormal robustness patterns. If
future morphological studies validated this pattern, this could be related to a change in growth in the early devel-
opment. Many experimental studies on modern squid show that several climate-related stressors can change the
morphology of squid's internal skeletons<sup>24,78–81</sup>. In particular, decreased pH in seawater is associated with abnor-
mal early development of statoliths<sup>81</sup> and cuttlebones in decabrachians as well as hypercalcification<sup>23,78–80,82</sup>. All
of these factors would affect buoyancy and locomotion, affecting the ability of young decabrachians to thrive as
normal<sup>23,24,81</sup>. Furthermore, experiments showed that sepiids reared under low feeding conditions develop wider
and shorter (more robust) internal shells<sup>83</sup>, and the PToB event shows evidence for both ocean acidification<sup>42</sup>
and low productivity<sup>62</sup>. Because the robustness variation of *P. bisulcata* correlates significantly with seawater
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Robustness (PC1) of (a) *C. longiforma* and (b) *P. bisulcata* life stages by subzone. The thick grey lines indicate the overall robustness, while the overlaying boxplots are split by ontogenetic stage (left to right: juvenile (light grey), neanic (medium grey), adult (dark grey)). The red background in the first Toarcian subzone (Mirabile) indicates the hyperthermal stress from the PToB event. (see sample sizes in Table 1).</figcaption></figure>

acidification, we suggest that hypercalcification and consequently underfeeding are plausible mechanisms for the morphological changes in *P. bisulcata*, that should be investigated further.

Given that both of the species that exhibit a significant increase in robustness are affected with arguably differing severity, we propose that *P. bisulcata* and *C. longiforma* may have had different thermal and/or bathymetric niches. The thermal tolerance of *C. longiforma* is unknown, but we argue that, because the higher relative abundance of this species at lower latitude localities⁸⁴, and because its robustness recovered when temperature at the
<table><caption>**Table 1.** Table of occurrences of species in the five ammonite subzones in the format: absolute abundance (adult|neanic|juvenile).</caption><thead><tr><th></th><th>P. bisulcata</th><th>C. longiforma</th><th>Bairstowius amaliae</th><th>P. milleri</th><th>Parapassaloteuthis sp.</th><th>Passaloteuthis sp. juv</th><th>Hastitidae indet</th><th>Acrocoelites sp.</th></tr></thead><tbody><tr><td>Elegantulum</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>2(2|0|0)</td></tr><tr><td>Semicelatum</td><td>3 (0|1|2)</td><td>13 (3|6|4)</td><td>-</td><td>3 (1|0|2)</td><td>4 (3|0|1)</td><td>1 (0|0|1)</td><td>-</td><td>-</td></tr><tr><td>Mirabile</td><td>10 (1|5|4)</td><td>17 (2|10|5)</td><td>-</td><td>-</td><td>-</td><td>-</td><td>3 (indet.)</td><td>-</td></tr><tr><td>Elisa</td><td>9 (1|7|1)</td><td>23 (8|12|3)</td><td>8 (1|6|1)</td><td>3 (1|2|0)</td><td>7 (3|0|4)</td><td>2 (0|0|2)</td><td>-</td><td>-</td></tr><tr><td>Solare</td><td>13 (2|3|8)</td><td>-</td><td>10 (2|4|4)</td><td>8 (3|5|0)</td><td>-</td><td>5 (0|0|5)</td><td>-</td><td>-</td></tr></tbody></table>

site decreased after the boundary event, *C. longiforma* appears better to be able to adapt to warmer habitats than *P. bisulcata*. Generally, southern range margins of high latitudinal species are believed to be set by interspecific competition, while northern limits are set by physical tolerance⁸⁵,⁸⁶. We interpret the responses of *P. bisulcata* as indicating that this species was adapted to relative cold to temperate environments, which is consistent with their oxygen isotope signals³⁰,⁵⁵. However, some studies interpreted the family Passaloteuthididae (including both *Catateuthis* and *Passaloteuthis*) as ‘all-rounders’ that could easily migrate vertically or spatially, because of the mix of specialised ecomorphological features they have⁵¹. Geochemical studies on *Passaloteuthis* specimens concluded that they preferentially dwelled in deep water without seasonal influence³⁰,⁵⁵. The epipelagic habitat in Peniche⁸⁷ and the high abundance of this species during the latest Pliensbachian sea-level regressive phase makes this seem unlikely⁴⁸. Additionally, the sea-level rise after the Pliensbachian-Toarcian boundary⁸⁷ likely created accommodation space that should have favoured *P. bisulcata*, aiding their recovery.

Our study found promising evidence for differential susceptibilities to CRS across the life stages of an extinct epipelagic predator, and possibly varying CRS susceptibility between species modulated by their thermal preferences. Further research on thermal niches and varying sensitivities among life stages in the fossil record can help anticipate CRS vulnerability in modern organisms⁸⁸ and inform conservation efforts.

## Material and methods

**Materials.** 930 belemnites were collected during multiple field trips to Peniche (Portugal) in 2016, 2017 and 2018. Specimens were collected from ten consecutive beds in three ammonite zones (*Emaciatum*, *Polymorphum* and *Levisoni*) covering five subzones (*Solare*, *Elisa*, *Mirabile*, *Semicelatum* and *Elegantulum*). Two different sampling techniques were employed to reduce the impact of sampling bias⁴⁸. First, specimens were quantitatively sampled by throwing a 1 m² frame three times per bed and collecting all belemnites and fragments that were found within that frame. Afterwards 30 complete specimens were sampled from each bed to facilitate taxonomic identification. Taxonomy and ontogenetic stages were identified based on traditionally used morphological traits⁴⁸,⁵⁷,⁸⁴,⁸⁹,⁹⁰. Out of the 930 belemnites⁴⁸,⁸⁴, all sufficiently preserved specimens were chosen for this study and represent the full samples’ taxonomic composition and ratio of ontogenetic stages well⁴⁸. The resulting dataset consists of 144 belemnites from seven species in three subfamilies (Table 1, supplementary Fig. S1).

These belemnites were scanned in the micro-CT scanner (v|Tome|X) of the GeoZentrum Nordbayern. After reconstructing the CT scans using VG Studio max, the open source software 3D Slicer⁹¹ was used to virtually clean the belemnite rostra from encrusting organisms, adhering sediment and the sediment filling the alveolus. The rostrum cavum was unevenly fractured at different lengths in most specimens. To correct for this taphonomic bias, we cut the rostrum cavum of all 3D models to 1/6th of the length of the rostrum solidum, as this was a length that was preserved in the available specimens. The specimens are stored at the Museu da Ciência da Universidade de Coimbra (Portugal) and the tomographic data can be accessed on MorphoBank (see Data Availability).

**Methods.** To attain proxies for palaeoenvironmental changes across the PToB, we used WebPlotDigitizer-4⁹² to extract δ¹⁸O (temperature proxy), δ¹³C (carbon cycle proxy) and δ¹¹B (seawater pH proxy) values from a geochemical study of the same outcrop⁴² and assigned them to the corresponding beds in our data⁴⁸ and consequently all belemnite rostra found within each bed.

The ‘geomorph’ R package⁵⁶ was used to manually place 22 3D-landmarks (Fig. 5) on the belemnite rostra. The landmarks that were placed on the dorsal and ventral sides of the rostra were defined as curve sliders. This means they would orientate themselves equidistantly along the flanks using minimum bending energy and can therefore counteract potential inaccuracies in placing them. A 3D-array, produced with the ‘abind’⁹³ package, of the landmarks of all specimens is used for the General Procrustes Analysis (GPA)⁵⁶,⁹⁴. This method rotates and scales the landmarks of all specimens in 3D space to the point of highest similarity. It then calculates the consensus shape (=mean) of all specimens⁹⁵. The ‘Procrustes distances’ between the specimens and the ‘consensus shape’ are used to perform a Principal Component Analysis (PCA).

First, we run a linear model on the relationship between the geometric mean and centroid size (supplementary Fig. S3). The former has been used as a body size proxy successfully before⁴⁸, and centroid size is the size component of a specimen that gets removed during the Procrustes transformation in geometric morphometrics⁵⁶. This is done to test, whether we could use centroid size as a reliable estimate of body size in this study.

To test for changes in morphology on different levels of assemblage organisation, we use Hedge's g from the ‘effsize’⁹⁶ package to quantify changes over time, and assess statistical significance with Mann-Whitney U tests between consecutive samples.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** The positions of the 22 landmarks on the sagittal (left figure) and transversal cross section (right figure) of a schematic belemnite rostrum: The rostrum cavum is cut to 1/6th of the length of the rostrum solidum, landmarks 11–22 are curve-sliding semi-landmarks that orientate themselves equidistantly along the ventral and dorsal side of the sagittal plane through the rostrum to counteract potential inaccuracies through human error.</figcaption></figure>

To test, whether a change in species composition of the assemblages has an impact on the morphological pattern, we first use an ANOVA to test for differences in robustness among the species. We subsequently use a bootstrapping method⁹⁷ to test whether changes in species composition across the subzones are random or not. A vector of the species names was drawn from 29 times (mean sample size on subzone level) with replacement, with the relative abundance of each species in the overall sample defined as the probability for it to be selected. This was repeated 500 times. For each of the 500 randomly calculated samples, the simulated relative abundances for the species were calculated. The range of relative abundances of the individual species is defined as the frame, within which changes in abundance can be considered random.

We then test for a within-species change in robustness, specifically across the boundary. The only two species that cross the boundary are *P. bisulcata* and *C. longiforma*. As both these species are known from the full dataset to have a significantly reduced body size after the boundary, we use Mann-Whitney U tests on those species' centroid sizes between consecutive subzone samples to test whether this works in our subset of specimens. We repeat this to test for changes in robustness (PC1) of these species through time.

To test, whether life stages within the species were affected differently we first employed a linear regression model, using developmental stage as an ordinal predictor variable. The ordinal data structure was used, because it shows trends in smaller data sets more reliably than categorical data structures. We then calculated AICc for both, the null-model and the model containing developmental stage to assess parsimony. Additionally, we run a mood's median test⁹⁸ on the robustness of the life stages of both dominant species among subzones. Because of lower sample sizes in the adults of both species, we group adult and neanic specimens together and compare them to juveniles. A post-hoc pairwise median test⁹⁹ is used to show significance of changes in morphology between consecutive subzones. For all statistical tests we presume that sample size is sufficient, when p-values return at least marginally significant results (<0.1), as insufficient sampling usually leads to type II errors and does not suggest false positive results. We include marginal significance, as lower p-values are less likely to be produced with smaller data sets.

The hypothesis, that accelerated growth at higher temperatures causes the change in morphology in adult specimens, is tested by examining the relationship between size and robustness in both species. We fit two linear regressions for each species. One for specimens found in beds during times of lower temperature (δ¹⁸O > -0.9) and one for belemnites that lived in a higher temperature environment (δ¹⁸O < -0.9) after the boundary. Low temperature and high temperature is defined based on which stratigraphic beds coincide with hyperthermal pulses of the Pliensbachian-Toarcian crisis. The "low temperature" and "high temperature" models are then compared through an ANOVA.

We use variation partitioning from the ‘vegan’ package¹⁰⁰ to calculate the independent effects that environmental changes and lithology have on the morphospace. For this we use four variables: δ¹⁸O, δ¹³C, δ¹¹B as environmental variables, and lithology (categorical: marl, limestone, marly limestone).

To assess the impact of δ¹⁸O, δ¹³C and δ¹¹B on the robustness on the assemblage scale as well as for *P. bisulcata* and *C. longiforma* individually, we use generalised least squares (GLS) modelling from the ‘nlme’ package in R¹⁰¹, correcting for the impact of lithology. This approach was decided on, because it corrects for auto-correlated residuals in the linear model. The chosen environmental variables were proxies for seawater temperature (δ¹⁸O), carbon cycle perturbations (δ¹³C) and seawater acidification (δ¹¹B). The correlation between the variables is <0.7¹⁰² (supplementary Fig. S5). We determine the most parsimonious model using AICc from the ‘nlme’ library in R¹⁰¹, which corrects for uneven sample sizes.

All analyses were conducted in R¹⁰³ (we use the R package ‘viridis’¹⁰⁴ to create figures accessible to colour blind people).

## Data availability

Uncompiled image stacks and mesh files of all belemnites in this study, including specimen identification numbers, have been deposited in a Morphobank repository (http://morphobank.org/per/?P3873). The data and R code used for data analysis have been deposited in a GitHub repository and will be made publicly available upon publication (https://github.com/PauSofN/Belemnite_Morphology_modulated_by_life_stages).
Received: 21 February 2021; Accepted: 30 June 2021

Published online: 14 July 2021

References

1. Reddin, C. J., Kocsis, Á. T., Aberhan, M. & Kiessling, W. Victims of ancient hyperthermal events herald the fates of marine clades and traits under global warming. *Glob. Chang. Biol.* **27**, 868–878 (2021).

2. Reddin, C. J., Kocsis, Á. T. & Kiessling, W. Marine invertebrate migrations trace climate change over 450 million years. *Glob. Ecol. Biogeogr.* **27**, 704–713 (2018).

3. Kordas, R. L., Harley, C. D. G. & O'Connor, M. I. Community ecology in a warming world: The influence of temperature on interspecific interactions in marine systems. *J. Exp. Mar. Bio. Ecol.* **400**, 218–226 (2011).

4. Poloczanska, E. S. et al. Responses of marine organisms to climate change across oceans. *Front. Mar. Sci.* **3**, 1–21 (2016).

5. Hanken, J. & Wake, D. B. Miniaturization of body size: Organismal consequences and evolutionary significance. *Annu. Rev. Ecol. Syst.* **24**, 501–519 (1993).

6. Sheridan, J. A. & Bickford, D. Shrinking body size as an ecological response to climate change. *Nat. Clim. Chang.* **1**, 401–406 (2011).

7. Forster, J., Hirst, A. G. & Atkinson, D. Warming-induced reductions in body size are greater in aquatic than terrestrial species. *Proc. Natl. Acad. Sci. U. S. A.* **109**, 19310–19314 (2012).

8. Ohlberger, J. Climate warming and ectotherm body size - from individual physiology to community ecology. *Funct. Ecol.* **27**, 991–1001 (2013).

9. Garilli, V. et al. Physiological advantages of dwarfing in surviving extinctions in high-CO₂ oceans. *Nat. Clim. Chang.* **5**, 678–682 (2015).

10. Daufresne, M., Lengfellner, K. & Sommer, U. Global warming benefits the small in aquatic ecosystems. *Proc. Natl. Acad. Sci. U. S. A.* https://doi.org/10.1073/pnas.0902080106 (2009).

11. Verberk, W. C. E. P. et al. Shrinking body sizes in response to warming: explanations for the temperature-size rule with special emphasis on the role of oxygen. *Biol. Rev.* **96**, 247–268 (2021).

12. Angilletta, M. J., Steury, T. D. & Sears, M. W. Temperature, growth rate, and body size in ectotherms: Fitting pieces of a life-history puzzle. *Integr. Comp. Biol.* **44**, 498–509 (2004).

13. Hoving, H. J. T. et al. Extreme plasticity in life-history strategy allows a migratory predator (jumbo squid) to cope with a changing climate. *Glob. Chang. Biol.* **19**, 2089–2103 (2013).

14. Angilletta, M. J. & Dunham, A. E. The temperature-size rule in ectotherms: simple evolutionary explanations may not be general. *Am. Nat.* **162**, 332–342 (2003).

15. Vinarski, M. V. On the applicability of Bergmann's rule to ectotherms: the state of the art. *Biol. Bull. Rev.* **4**, 232–242 (2014).

16. Atkinson, D. Temperature and organism size: a biological law for organisms?. *Adv. Ecol. Res.* **25**, 1–58 (1994).

17. Atkinson, D. Effects of temperature on the size of aquatic ectotherms: Exceptions to the general rule. *J. Therm. Biol.* **20**, 61–74 (1995).

18. Forster, J. & Hirst, A. G. The temperature-size rule emerges from ontogenetic differences between growth and development rates. *Funct. Ecol.* **26**, 483–492 (2012).

19. Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. *Science* **293**, 2248–2251 (2001).

20. Schulte, P. M. The effects of temperature on aerobic metabolism: Towards a mechanistic understanding of the responses of ectotherms to a changing environment. *J. Exp. Biol.* **218**, 1856–1866 (2015).

21. Gillooly, J. F., Charnov, E. L., West, G. B., Savage, V. M. & Brown, J. H. Effects of size and temperature on developmental time. *Nature* **417**, 70–73 (2002).

22. Riemer, K., Anderson-Teixeira, K. J., Smith, F. A., Harris, D. J. & Ernest, S. K. M. Body size shifts influence effects of increasing temperatures on ectotherm metabolism. *Glob. Ecol. Biogeogr.* **27**, 958–967 (2018).

23. Rosa, R. et al. Ocean warming enhances malformations, premature hatching, metabolic suppression and oxidative stress in the early life stages of a keystone squid. *PLoS ONE* **7**, e38282 (2012).

24. Pecl, G. T. & Jackson, G. D. The potential impacts of climate change on inshore squid: Biology, ecology and fisheries. *Rev. Fish Biol. Fish.* **18**, 373–385 (2008).

25. Twitchett, R. J. The Lilliput effect in the aftermath of the end-Permian extinction event. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **252**, 132–144 (2007).

26. Harries, P. J. & Knorr, P. O. What does the 'Lilliput Effect' mean?. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **284**, 4–10 (2009).

27. Metcalfe, B., Twitchett, R. J. & Price-Lloyd, N. Changes in size and growth rate of 'Lilliput' animals in the earliest Triassic. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **308**, 171–180 (2011).

28. Chu, D. et al. Lilliput effect in freshwater ostracods during the Permian-Triassic extinction. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **435**, 38–52 (2015).

29. Urbanek, A. Biotic crises in the history of upper silurian graptoloids: a palaeobiological model. *Hist. Biol.* https://doi.org/10.1080/10292389309380442 (1993).

30. Urlichs, M. Stunting in invertebrates from the type area of the Cassian Formation (Early Carnian) of the dolomites (Italy). *GeoAlp* **8**, 164–169 (2011).

31. Morten, S. D. & Twitchett, R. J. Fluctuations in the body size of marine invertebrates through the Pliensbachian-Toarcian extinction event. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* https://doi.org/10.1016/j.palaeo.2009.08.023 (2009).

32. Piazza, V., Ullmann, C. V. & Aberhan, M. Temperature-related body size change of marine benthic macroinvertebrates across the Early Toarcian Anoxic Event. *Sci. Rep.* **10**, 1–13 (2020).

33. Calosi, P., Putnam, H. M., Twitchett, R. J. & Vermandele, F. Marine metazoan modern mass extinction: improving predictions by integrating fossil, modern, and physiological data. *Ann. Rev. Mar. Sci.* **11**, 369–390 (2019).

34. Gerber, S. Comparing the differential filling of morphospace and allometric space through time: the morphological and developmental dynamics of Early Jurassic ammonoids. *Paleobiology* **37**, 369–382 (2011).

35. Pálfy, J. & Smith, P. L. Synchrony between Early Jurassic extinction, oceanic anoxic event, and the Karoo-Ferrar flood basalt volcanism. *Geology* **28**, 747–750 (2000).

36. Caruthers, A. H., Smith, P. L. & Gröcke, D. R. The Pliensbachian-Toarcian (Early Jurassic) extinction, a global multi-phased event. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **386**, 104–118 (2013).

37. Wignall, P. B. Large igneous provinces and mass extinctions. *Earth Sci. Rev.* **53**, 1–33 (2001).

38. Percival, L. M. E. et al. Globally enhanced mercury deposition during the end-Pliensbachian extinction and Toarcian OAE: A link to the Karoo-Ferrar Large Igneous Province. *Earth Planet. Sci. Lett.* **428**, 267–280 (2015).

39. Foster, G. L., Hull, P., Lunt, D. J. & Zachos, J. C. Placing our current 'hyperthermal' in the context of rapid climate change in our geological past. *Phil. Trans. R. Soc. A* **376**, 20170086 (2018).

40. Suan, G. et al. Secular environmental precursors to Early Toarcian (Jurassic) extreme climate changes. *Earth Planet. Sci. Lett.* **290**, 448–458 (2010).
41. Fantasia, A. et al. Global versus local processes during the Pliensbachian-Toarcian transition at the Peniche GSSP, Portugal: A multi-proxy record. *Earth Sci. Rev.* **198**, 102932 (2019).

42. Müller, T. et al. Ocean acidification during the early Toarcian extinction event: Evidence from Boron isotopes in brachiopods. *Geology* **48**, 1184–1188 (2020).

43. Suan, G., Mattioli, E., Pittet, B., Mailliot, S. & Lécuyer, C. Evidence for major environmental perturbation prior to and during the Toarcian (Early Jurassic) oceanic anoxic event from the Lusitanian Basin Portugal. *Paleoceanography* **23**, A1202 (2008).

44. Dera, G. et al. High-resolution dynamics of early Jurassic marine extinctions: The case of Pliensbachian-Toarcian ammonites (Cephalopoda). *J. Geol. Soc. London* **167**, 21–33 (2010).

45. Dera, G. et al. Climatic ups and downs in a disturbed Jurassic world. *Geology* **39**, 215–218 (2011).

46. Miguez-Salas, O., Rodriguez-Tovar, F. J. & Duarte, L. V. Selective incidence of the Toarcian oceanic anoxic event on macroinvertebrate marine communities: a case from the Lusitanian basin Portugal. *Lethaia* **50**, 548–560 (2017).

47. Correia, V. F., Riding, J. B., Duarte, L. V., Fernandes, P. & Pereira, Z. The palynological response to the Toarcian Oceanic Anoxic Event (Early Jurassic) at Peniche, Lusitanian Basin, western Portugal. *Mar. Micropaleontol.* **137**, 46–63. https://doi.org/10.1016/j.marmicro.2017.10.004 (2017).

48. Rita, P., Nätscher, P., Duarte, L. V., Weis, R. & De Baets, K. Mechanisms and drivers of belemnite body-size dynamics across the Pliensbachian-Toarcian crisis. *R. Soc. Open Sci.* **6**, 190494 (2019).

49. Caswell, B. A. & Coe, A. L. The impact of anoxia on pelagic macrofauna during the Toarcian Oceanic Anoxic Event (Early Jurassic). *Proc. Geol. Assoc.* **125**(4), 383–391. https://doi.org/10.1016/j.pgeola.2014.06.001 (2014).

50. Ullmann, C. V., Thibault, N., Ruhl, M., Hesselbo, S. P. & Korte, C. Effect of a Jurassic oceanic anoxic event on belemnite ecology and evolution. *Proc. Natl. Acad. Sci. U. S. A.* **111**, 10073–10076 (2014).

51. Dera, G., Toumoulin, A. & de Baets, K. Diversity and morphological evolution of Jurassic belemnites from South Germany. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **457**, 80–97 (2016).

52. Neige, P., Weis, R. & Fara, E. Ups and downs of belemnite diversity in the Early Jurassic of Western Tethys. *Palaeontology* **64**, 263–283 (2021).

53. Rita, P., De Baets, K. & Schlott, M. Rostrum size differences between Toarcian belemnite battlefields. *Foss. Rec.* **21**, 171–182 (2018).

54. Rita, P. et al. Biogeographic patterns of belemnite body size responses to episodes of environmental crisis. *PeerJ Prepr.* (2019).

55. Hoffmann, R. & Stevens, K. The palaeobiology of belemnites - foundation for the interpretation of rostrum geochemistry. *Biol. Rev.* **95**, 94–123 (2020).

56. Adams, D. C. & Otárola-Castillo, E. Geomorph: An r package for the collection and analysis of geometric morphometric shape data. *Methods Ecol. Evol.* **4**, 393–399 (2013).

57. Schlegelmilch, R. Die Belemniten des süddeutschen Jura. *Die Belemniten des süddeutschen Jura* https://doi.org/10.1007/978-3-8274-3083-0 (1998).

58. McArthur, J. M. et al. Sr-isotope stratigraphy (87Sr/86Sr) of the lowermost Toarcian of Peniche, Portugal, and its relation to ammonite zonations. *Newsletters Stratigr.* **53**, 297–312 (2020).

59. Hesselbo, S. P., Jenkyns, H. C., Duarte, L. V. & Oliveira, L. C. V. Carbon-isotope record of the Early Jurassic (Toarcian) Oceanic Anoxic Event from fossil wood and marine carbonate (Lusitanian Basin, Portugal). *Earth Planet. Sci. Lett.* **253**, 455–470 (2007).

60. Klug, C., Schweigert, G., Fuchs, D., Kruta, I. & Tischlinger, H. Adaptations to squid-style high-speed swimming in Jurassic belemnitids. *Biol. Lett.* **12**, 20150877 (2016).

61. Seibel, B. A., Thuesen, E. V., Childress, J. J. & Gorodezky, L. A. Decline in pelagic cephalopod metabolism with habitat depth reflects differences in locomotory efficiency. *Biol. Bull.* **192**, 262–278 (1997).

62. Mattioli, E., Pittet, B., Petitpierre, L. & Mailliot, S. Dramatic decrease of pelagic carbonate production by nannoplankton across the Early Toarcian anoxic event (T-OAE). *Glob. Planet. Change* **65**, 134–145 (2009).

63. Chamberlain, J. A. Locomotion in ancient seas: Constraint and opportunity in Cephalopod adaptive design. *Geobios* **15**, 49–61 (1993).

64. Rexfort, A. & Mutterlose, J. The role of biogeography and ecology on the isotope signature of cuttlefishes (Cephalopoda, Sepiidae) and the impact on belemnite studies. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **284**, 153–163 (2009).

65. Holland, S. M. The quality of the fossil record: a sequence stratigraphic perspective. *Paleobiology* **26**, 148–168 (2000).

66. Holland, S. M. The non-uniformity of fossil preservation. *Philos. Trans. R. Soc. B Biol. Sci.* **371**, 2 (2016).

67. Korn, D. Impact of environmental perturbations o heterochronic developments in Palaeozoic ammonoids. *Evol. Chang. Heterochrony* 245–260 (1995).

68. Yacobucci, M. M. Plasticity of developmental timing as the underlying cause of high speciation rates in ammonoids. in *Advancing research on living and fossil cephalopods* 59–76 (Springer, Boston, MA, 1999).

69. Landman, N. H. & Gyssant, J. R. Heterochrony and ecology in Jurassic and Cretaceous ammonites. *Geobios* **26**, 247–255 (1993).

70. McNamara, K. J. Heterochrony: the evolution of development. *Evol. Educ. Outreach* **5**, 203–218 (2012).

71. Dahlke, F. T., Wohlrab, S., Butzin, M. & Pörtner, H. O. Thermal bottlenecks in the life cycle define climate vulnerability of fish. *Science* **369**, 65–70 (2020).

72. Pörtner, H. O. & Farrell, A. P. Ecology: Physiology and climate change. *Science* https://doi.org/10.1126/science.1163156 (2008).

73. Pimentel, M. S. et al. Impact of ocean warming on the early ontogeny of cephalopods: A metabolic approach. *Mar. Biol.* **159**, 2051–2059 (2012).

74. Komoroske, L. M. et al. Ontogeny influences sensitivity to climate change stressors in an endangered fish. *Conserv. Physiol.* **2**, 1–13 (2014).

75. Pörtner, H. O., Bock, C. & Mark, F. C. Oxygen- & capacity-limited thermal tolerance: Bridging ecology & physiology. *J. Exp. Biol.* **220**, 2685–2696 (2017).

76. Harnik, P. G., Simpson, C. & Payne, J. L. Long-term differences in extinction risk among the seven forms of rarity. *Proc. R. Soc. B Biol. Sci.* **279**, 4969–4976 (2012).

77. Reddin, C. J., Kocsis, Á. T. & Kiessling, W. Climate change and the latitudinal selectivity of ancient marine extinctions. *Paleobiology* **45**, 70–84 (2019).

78. Dorey, N. et al. Ocean acidification and temperature rise: Effects on calcification during early development of the cuttlefish *Sepia officinalis*. *Mar. Biol.* **160**, 2007–2022 (2013).

79. Sigwart, J. D. et al. Elevated pCO2 drives lower growth and yet increased calcification in the early life history of the cuttlefish *Sepia officinalis* (Mollusca: Cephalopoda) Julia. *ICES J. Mar. Sci.* **73**, 970–980 (2016).

80. Gutowska, M. A., Melzner, F., Pörtner, H. O. & Meier, S. Cuttlebone calcification increases during exposure to elevated seawater pCO2 in the cephalopod *Sepia officinalis*. *Mar. Biol.* **157**, 1653–1663 (2010).

81. Kaplan, M. B., Mooney, T. A., McCorkle, D. C. & Cohen, A. L. Adverse effects of ocean acidification on early development of squid (*Doryteuthis pealeii*). *PLoS ONE* **8**, e63714 (2013).

82. Otjacques, E. et al. Cuttlefish buoyancy in response to food availability and ocean acidification. *Biology* (Basel). https://doi.org/10.3390/biology9070147 (2020).

83. Neige, P. & Boletzky, S. Morphometrics of the shell of three *Sepia* species (Mollusca: Cephalopoda): Intra- and interspecific variation. *Zool. Beitraege*. **38**, 137–156 (1997).
84. Rita, P., Weis, R., Duarte, L. V. & De Baets, K. Taxonomical diversity and palaeobiogeographical affinity of belemmites from the Pliensbachian-Toarcian GSSP (Lusitanian Basin, Portugal). *Pap. Palaeontol.* https://doi.org/10.1002/spp2.1343 (2020).

85. MacArthur, R. H. *Geographical ecology: patterns in the distribution of species*. (Princeton University Press, 1972).

86. Gaston, K. J. *The structure and dynamics of geographic ranges*. (Oxford University Press on Demand, 2003).

87. Duarte, L. Sequence stratigraphy and depositional setting of the Pliensbachian and Toarcian marly limestones in the Lusitanian Basin Portugal. *Ciências da Terra* **16**, 17-23 (2007).

88. Reddin, C. J., Nätscher, P. S., Kocsis, Á. T., Pörtner, H.-O. & Kiessling, W. Marine clade sensitivities to climate change conform across timescales. *Nat. Clim. Chang.* **10**, 249–253 (2020).

89. Doyle, P. & Kelly, R. A. *The Jurassic and Cretaceous belemmites of Kong Karls Land, Svalbard*. (Norsk Polarinstitutt Oslo, 1988).

90. Doyle, P. New records of dimitobelid belemmites from the cretaceous of james ross island Antarctica. *Alcheringa* **14**, 159–175 (1990).

91. Fedorov, A. et al. 3D Slicer as an image computing platform for the Quantitative Imaging Network. *Magn. Reson. Imaging* https://doi.org/10.1016/j.mri.2012.05.001 (2012).

92. Rohatgi, A., Rehberg, S. & Stanojevic, Z. Webplotdigitizer: Version 4.1 of Webplotdigitizer. (2018). https://doi.org/10.5281/zenodo.1137880.

93. Plate, T. & Heiberger, R. Package 'abind'. (2016).

94. Collyer, M. L. & Adams, D. C. RRPP: An r package for fitting linear models to high-dimensional data using residual randomization. *Methods Ecol. Evol.* **9**, 1772–1779 (2018).

95. Sherratt, E. Quick Guide to Geomorph v. 2.0. public.iastate.edu (2014).

96. Torchiano, M. Package 'effsize'. (2020).

97. Gotelli, N. J., Dorazio, R. M., Ellison, A. M. & Grossman, G. D. Detecting temporal trends in species assemblages with bootstrapping procedures and hierarchical models. *Philos. Trans.* **365**, 3621–3631 (2010).

98. Hervé, M. Package 'RVAideMemoire'. (2021).

99. Mangiafico, S. Package 'rcompanion'. (2021).

100. Oksanen, J. et al. Package 'vegan' title community ecology package. *Commun. Ecol. Packag.* **2**, 1–297 (2019).

101. Pinheiro, J. et al. Package 'nlme'. (2021).

102. Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. *Ecography* (Cop.) **36**, 27–46 (2013).

103. R Core Team. R: A Language and Environment for Statistical Computing. *Vienna, Austria* (2019).

104. Garnier, S., Ross, N., Rudis, B. & Sciaini, M. Package 'viridis'. (2021).

## Acknowledgements

This research was performed in the framework of the DFG Research Unit TERSANE (FOR 2332; grant numbers BA 5148/1-1, BA 5148/1-2 to PD Dr. Kenneth De Baets, and AB 109/11-1 to PD Dr. Martin Aberhan) and IGCP 655 (IUGS-UNESCO). We want to thank Alejandra Gomez and Christian Schulbert (FAU) for CT-scanning belemmites and Dr. Emilia Jarochowska for insightful discussions. We are greatly indebted to Luís Vitor Duarte who guided us through the nice outcrops during the expeditions to collect the belemnite specimens, Birgit Leipner-Mata for helping to prepare our specimens and the Museu da Ciência da Universidade de Coimbra (Portugal) for safeguarding the specimens.

## Author contributions

K.D.B. and P.R. collected the specimens in the field. P.S.N. and P.R. CT-scanned the specimens. P.S.N. and K.D.B. designed the study and formed the hypotheses. P.S.N. compiled and analysed the data with support from C.J.R. P.S.N. and K.D.B. produced the figures. P.S.N. wrote the manuscript with support from K.D.B., C.J.R. and G.D.. All authors worked on drafts and approved the final version of the manuscripts.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-93850-0.

**Correspondence** and requests for materials should be addressed to P.S.N.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.