Bivalve body-size distribution through the Late Triassic mass extinction event

Abstract The synergic relationship between physiology, ecology, and evolutionary process makes the body-size distribution (BSD) an essential component of the community ecology. Body size is highly susceptible to environmental change, and extreme upheavals, such as during a mass extinction event, could exert drastic changes on a taxon's BSD. It has been hypothesized that the Late Triassic mass extinction event (LTE) was triggered by intense global warming, linked to massive volcanic activity associated with the Central Atlantic Magmatic Province. We test the effects of the LTE on the BSD of fossil bivalve assemblages from three study sites spanning the Triassic/Jurassic boundary in the United Kingdom. Our results show that the effects of the LTE were rapid and synchronous across sites, and the BSDs of the bivalves record drastic changes associated with species turnover. No phylogenetic signal of size selectivity was recorded, although semi-infaunal species were apparently most susceptible to change. Each size class had the same likelihood of extinction during the LTE, which resulted in a platykurtic BSD with negative skew. The immediate postextinction assemblage exhibits a leptokurtic BSD, although with negative skew, wherein surviving species and newly appearing small-sized colonizers exhibit body sizes near the modal size. Recovery was relatively rapid (~100 kyr), and larger bivalves began to appear during the pre-Planorbis Zone, despite recurrent dysoxic/anoxic conditions. This study demonstrates how a mass extinction acts across the size spectrum in bivalves and shows how BSDs emerge from evolutionary and ecological processes.


Introduction
Body size is a key trait of all organisms and correlates with almost every aspect of their physiology, ecology, and evolution (Roy 2008). An organism's size reflects how it acquires and processes energy and provides valuable information regarding the resource partitioning between species along a range size (Peters 1983;Ernest 2005;Anderson-Teixeira et al. 2009). The body-size distribution (BSD) of species in a guild, clade, or assemblage reflects their niche requirements and is the macroecological and macroevolutionary outcome of the many drivers that constrain body size (Calder 1984;Roy et al. 2000). The "null" shape of the BSD fits to a lognormal distribution (McKinney 1990;McShea 1994;Marquet and Taper 1998;Maurer 1998Maurer , 2003Kozłowski and Gawelczyk 2002;Maurer and Marquet 2013), in which the minimum and maximum sizes are constrained by energetic requirements, biotic interactions, and environmental drivers (Brown et al. 1978(Brown et al. , 1993Schmidt-Nielsen 1984;Marquet and Taper 1998;Allen et al. 2006;Maurer and Marquet 2013;Smith et al. 2016) and the height (or modal size) is the most common bin size used by different species or individuals (Maurer and Marquet 2013). BSDs have been described for several marine assemblages (Warwick 1984;Kendall et al. 1997;Roy et al. 2000;Hildrew et al. 2007;Baltanas and Danielopol 2013;Gearty et al. 2018) and tend to change in response to nutrient availability, the type of ecosystem, and the level of disturbance (Warwick 1984;McClain et al. 2012). For instance, lower levels of available energy result in assemblages of smaller-sized organisms (e.g., right-skewed BSD), with strategies that enable them to survive with lower levels of foraging and/or with mechanisms to cope with starvation (McClain 2004;McClain et al. 2012). Assemblages that occur under higher levels of productivity have diverse taxa of larger-sized individuals and expensive life strategies (e.g., left-skewed BSD) (Rex and Etter 1998;Roy et al. 2000;McClain 2004;McClain et al. 2012McClain et al. , 2018. Therefore, the BSD may provide valuable insights into the determinants of species structure through episodes of environmental change, such as past biotic crises. Body size is highly susceptible to environmental change (Gardner et al. 2011). In marine animals, global warming is expected to lead to a reduction in body size due to the impacts of temperature rise, lower concentrations of dissolved oxygen, lower pH, and changes in productivity (Sheridan and Bickford 2011). Body-size reduction has also been demonstrated in a range of marine taxa during several past mass extinction events (Urbanek 1993;Twitchett and Barras 2004;He et al. 2010;Martinez-Diaz et al. 2016). The temporary reduction in the body size of surviving species in the immediate aftermath of an extinction event was termed the "Lilliput Effect" by Urbanek (1993), and although there are several potential reasons why such a phenomenon could be recorded in the fossil record (Twitchett 2007;Harries and Knorr 2009), in most cases, the smallest sizes are associated with peak temperatures and global warming, and so likely reflect responses to environmental change. During the recovery from mass extinction events, the sizes of surviving taxa and newly originating taxa increase as environmental conditions ameliorate (Payne 2005; Barras and Twitchett 2007;Novack-Gottshall 2008;Metcalfe et al. 2011;He et al. 2016;Sigurdsen and Hammer 2016;Brom et al. 2018;Atkinson and Wignall 2020;Foster et al. 2020).
Theoretical models would suggest that BSD would emerge throughout a cladogenetic diffusion model of species body size (McKinney 1990(McKinney , 1997McShea 1994;Clauset and Erwin 2008), in which small-bodied species are expected to be common relative to larger-bodied species (Stanley 1973(Stanley , 1979Hayami 1978;Gould 1988;McShea 1994;Monroe and Bokma 2013) (Fig. 1). However, the size diffusion process can be punctuated by episodes of biotic crisis (Gould 1988;McKinney 1990), and two possible BSD shapes may result depending on the selective nature of the disturbance (Fig. 1). If the extinction is size selective, then a curtailed BSD would result due to extinction of larger taxa, whereas extinction that is random with respect to size could lead to a platykurtic size distribution (McKinney 1990) (Fig. 1). Although a number of previous studies have addressed size change during extinction events, most of the work has focused on changes within individual taxa. Hitherto, little attention has been paid to the relationship between the dynamics of within-lineage size change and the effect on the evolution of the overall BSD of species arranged along the entire size spectrum, or how the overall BSD is modified by a mass extinction event.
Here, we study how the BSD of bivalves was affected by the Late Triassic mass extinction event (LTE) (∼201.5 Ma; Wotzlaw et al. 2014), the second biggest biodiversity loss of the Phanerozoic (Kocsis et al. 2019). The LTE is directly linked to an intense greenhouse event (McElwain et al. 1999;Berner and Beerling 2007;Kidder and Worsley 2010;Davies et al. 2017;Dunhill et al. 2018) triggered by massive injection of CO 2 into the ocean-atmosphere system (van de Schootbrugge et al. 2008;Blackburn et al. 2013) during eruption of the Central Atlantic Magmatic Province (CAMP) (Fig. 1) (Zaffani et al. 2018). Secondary environmental effects resulting from global warming have also been inferred, including ocean acidification (Crne et al. 2011), widespread anoxia (Jaraula et al. 2013;He et al. 2020), and rapid sea-level change (Hallam and Wignall 1999). Evidence of body-size changes has been recorded in a variety of taxa at different spatial and temporal scales during the LTE (Dommergues et al. 2002; Barras and Twitchett 2007;Atkinson and Wignall 2020), but no studies have examined how the ecological collapse associated with the LTE modified the BSD of marine organisms at local and regional scales.
Using data collected from three wellcorrelated study sites within the United Kingdom, we assess (1) whether the body sizes of individual bivalve taxa and the BSD as a whole changed significantly through the LTE, (2) whether the patterns and timing of changes varies between sites, and (3) whether the LTE was a size-selective extinction, by testing the null hypothesis that the BSD of bivalves does not exhibit measurable changes through the LTE.

Sample Sites and Stratigraphy
Fossil bivalves measured in this study came from three Triassic-Jurassic marine successions in the United Kingdom ( Supplementary  Fig. S1). The first study site was St. Audrie's Bay, Somerset, southwest England (51°1 0 ′ 0.36 ′′ N, 3°7 ′ 8.2 ′′ W) (Fig. 2), which has one of the most complete, well-known, and wellpreserved Triassic-Jurassic marine fossil records (Warrington et al. 2008). The second site sampled was Larne, County Antrim, Northern Ireland (54°51 ′ 26 ′′ N, 5°48 ′ 18 ′′ W), which has a ∼115-m-thick stratigraphically expanded section that spans from the upper part of the Mercia Mudstone Group (Triassic, Norian) through to the Bucklandi Zone of the Lias Group (Jurassic, Sinemurian) and is cut by minor faults (Simms and Jeram 2007). Only 62 m were sampled, from the base of Westbury Formation to the Liasicus Zone of the Blue Lias Formation ( Supplementary Fig. 1 McKinney (1990), the scheme shows how cladogenetic diffusion processes give rise to BSDs. Small sizes are more common and tend to accumulate around the modal size, while clades of larger sizes arriving gradually into the community expand the tail of the distribution. An extinction event punctuates the diffusion process, generating two outcomes depending on whether the extinction is selective or not. Larger sizes are eliminated under selective extinctions, while a nonselective event generates a flatter BSD due to a random probability of extinction.
( Supplementary Fig. 1, Supplementary Files). One sample of 1.5 ± 0.2 kg was collected from each major lithology (limestone and mudstone) on average every 1.98 ± 0.5 m. In addition, to increase the sample size of each stratigraphic horizon, individual measurements were made in the field from bedding plane exposures using a tape measure (SE ± 0.5 mm).
Samples from the different lithologies were processed in different ways. Limestone rock samples of 1.5 ± 0.2 kg (with dimensions of ca. 20 cm length × 20 cm width ×15 cm height) were cut perpendicular to the bedding plane surface, generating slabs approximately 8 cm wide (see Supplementary Fig. 2A). Using a chisel-tipped hammer, each slab (20 × 15 × 8 cm) was then split parallel to the bedding plane, which generated sampling plates of approximately 2 × 20 × 8 cm with specimens exposed on the surfaces of each slab for potential measurement (see Supplementary Fig. 2A). This technique minimized damage to the specimens in order to preserve entire individuals for measurement.
In contrast, sandstone and mudstone samples were broken into ∼50 cm 2 chips, and FIGURE 2. Stratigraphy, location, and paleogeography of study sites sampled in the United Kingdom. A, Sampled localities in this study (white squares). B, Global view of Pangea showing the paleogeographic location of United Kingdom during the Rhaetian-Hettangian (white square). C, Generalized stratigraphic and biostratigraphic scheme of the Triassic-Jurassic of the United Kingdom, with a carbon isotope curve from the section at St. Audrie's Bay, UK. The shaded gray area indicates the extinction horizon, which coincides with the first negative carbon isotope excursion (CIE) (Hesselbo et al. 2004). Line with arrowheads indicates the position of the main negative carbon isotope excursion (Main CIE) (Hesselbo et al. 2004); the durations of each unit were obtained from Ruhl et al. (2010), Bonis et al. (2010), Wotzlaw et al. (2014), andWeedon et al. (2019). CAMP, Central Atlantic Magmatic Province.
water was used to soak the rock and separate the fossils without damaging them. All the bivalves were classified from published information (Swift and Martill 1999;Hodges 2000;Lord and Davis 2010). In all cases, measurements of the length (maximum distance on the anterior-posterior axis) and height (maximum distance on the dorsal-ventral axis) were recorded only from complete, wellpreserved specimens ( Supplementary Fig. 2B, C). All morphometric variables were measured using electronic calipers (SE ± 0.01 mm).

Data Analysis
The body size of each specimen was calculated as the geometric mean (i.e., the square root of the product of length and width) and log 2 transformed to reduce the heteroscedasticity and to facilitate graphic representation (Piazza et al. 2019). Changes in the body sizes of bivalves through the LTE were evaluated in four different ways.
Body-Size Distribution across Space and Time.-To evaluate body-size changes through the LTE and to test whether any changes occurred synchronously at the different localities, mean body sizes per species were pooled within each stratigraphic unit. The shapes of the BSDs of each stratigraphic unit were examined graphically by plotting standard statistical metrics such as the coefficient of variation, skewness, maximum and minimum values, and median. A Shapiro-Wilk test was performed to test the normality of each distribution. A Kruskal-Wallis test was performed to detect significant differences in median values between each stratigraphic unit. A test of homogeneity of variance (Bartlett's test) was performed on each stratigraphic unit to evaluate whether the changes in the variance between localities were synchronous. Finally, the residuals of the body sizes of each stratigraphic unit were calculated and plotted to graphically assess changes in BSD. Each residual value was calculated as the difference between the grand mean ( X = 12.63 mm; i.e., the average size of entire dataset) minus the average body size of each species present in each stratigraphic unit. To test significant differences in the residual distribution of each stratigraphic unit, 95% confidence intervals of the residuals were estimated by an iterative random resampling (9999 iterations). Differences were considered significant when the residual values fell outside the confidence intervals.
Taxonomic and Size Selectivity.-To test whether the LTE was size selective, each bivalve taxon was assigned to one of three extirpation categories: (1) survivors, that is, those recorded in at least one locality and stratigraphic unit from both the pre-extinction and the recovery intervals; (2) colonizers, that is, those taxa that only occur after the extinction event (i.e., between the Cotham Member of the Lilstock Formation and Angulata Zone of the Blue Lias Formation). The extinction horizon coincides with the onset of a negative δ 13 C org excursion that spans from the upper Cotham Member to the lower Langport Member (20-40 kyr)  and is characterized by high species turnover and low diversity Ruhl et al. 2010;Thibodeau et al. 2016;) (see Fig. 2C and Supplementary  Fig. S1). Finally, (3) extirpated taxa (hereafter referred to as "victims") are those that recorded their last appearance either before, during, or after the LTE (i.e., from the upper Westbury Formation to the Langport Member of the Lilstock Formation). To detect significant changes in the BSD across the LTE, the coefficient of variation (i.e., standard deviation divided by the mean) was calculated for each category in each stratigraphic unit. In addition, residual values were calculated for each stratigraphic unit to track and assess the changes in BSD through time. For all dispersion parameters, the statistical differences between different stratigraphic units were evaluated through the estimation of 95% confidence intervals estimated by iterative random resampling (9999 iterations). Differences between successive bins were considered significant when their confidence intervals did not overlap.
Studies of the fossil record have demonstrated that selective extinction against specific lineages or ecological traits would drive important ecological and evolutionary modifications (Hallam and Wignall 1997;Jablonski 2008;Janevski and Baumiller 2009;Harnik 2011;Orzechowski et al. 2015;Payne et al. 2016a,b). Substantial differences with respect to trait between taxonomic groups could change the distribution of the species within genera, families, and orders (Jablonski and Raup 1995;Aberhan and Baumiller 2003;Smith and Roy 2006;Rivadeneira and Marquet 2007;Finnegan et al. 2012;Ros and Echevarria 2012;Harnik et al. 2014;Dishon et al. 2020). To test existence of phylogenetic effect in size selectivity in bivalves across the LTE, two methodological steps based on Rivadeneira and Marquet (2007) were implemented. First, a Moran's autocorrelogram analysis was used to evaluate whether the differences in body size among species are related to their taxonomic distance (Gittleman and Kot 1990;Rivadeneira and Marquet 2007). It has been suggested that this method is a good analytical strategy to assess whether a phylogenetic control may exist on quantitative traits when phylogenetic information is not available (Gittleman and Kot 1990). In this case, a typical phylogenetic signature is evident by the decay in the rescaled Moran's I toward higher taxonomic levels (Gittleman and Kot 1990;Rivadeneira and Marquet 2007). High phylogenetic inertia is interpreted as a high level of similarity among taxa or low variation of the trait between lineages (Gittleman and Kot 1990). If phylogenetic inertia is identified at a taxonomic level, the second step is to identify the selectivity. For that, the body size of each species in each extirpation category is regressed against the mean body size of the taxonomic level with the most variation. If the number of extirpated taxa, for instance, increases significantly along with larger order, then we would be recording a selective extirpation. Linear square regressions were fit to each trait value, and a linear regression with a 0 intercept was fit to all data points. Finally, all the slopes and intercepts were compared. Moran's autocorrelogram analysis was estimated using the package APE (Paradis and Schliep 2019), implemented in R v. 3.6.0 (R Core Team 2019).
Density Distribution and Generalized Linear Model of Different Tiers.-All the recorded bivalves were classified according to their tiering habit (i.e., semi-infaunal, infaunal, and epifaunal) from information available in the Paleobiology Database (PBDB; https://paleobiodb.org) and .
The body sizes of all individuals of each species were pooled to capture the complete size distribution of each tier in each stratigraphic unit. The modality of the BSDs in each stratigraphic unit was examined by plotting the frequency distribution of each mode of life with a bin size of 0.25 on a log 2 scale. Subsequently, a continuous distribution to each mode of life was fit by using a combination of kernel density estimation and smoothed bootstrap resampling (based on 1000 randomizations) (Manly 1996). This procedure approximates the probability density function of a random variable by determining the optimum modality of the data, testing whether a distribution with k + 1 modes fits significantly better than a distribution with k modes (Silverman 1981;Manly 1996). In addition, measures of skewness and kurtosis were estimated for each distribution. A Lilliefors (Kolmogorov-Smirnov) test was performed to test the normality of each distribution, and a Kolmogorov-Smirnov two-sample test (K-S test) was used to compare BSDs between pairs of consecutive stratigraphic units. In this case, the probability values ( p) were adjusted by the Bonferroni correction.
Generalized linear models (GLMs), assuming a Gaussian distribution with a log-link function, were used to examine differences in mean body size between the different tiers in each stratigraphic level. A full factorial design with two factors was implemented: factor 1 (tiering) had three levels (i.e., infaunal, semiinfaunal, and epifaunal); factor 2 (stratigraphic unit) had seven levels (horizons); while body size was used as the response variable. A Gaussian distribution was fit by Akaike information criterion (AIC = −1420.28), using the package propagate (Spiess 2018) implemented in R v. 3.6.0 (R Core Team 2019). Finally, a Newman-Keuls test was performed as an a posteriori test to determine which groups generated the differences.
Null Model.-An unweighted null model was constructed to assess whether any changes in recorded body size could be explained purely by random assembly. A null model is based on randomization of ecological data to evaluate patterns of some ecological or evolutionary process of interest. Here, the randomization was designed to produce a pattern that would be expected in the absence of a particular ecological mechanism, with the aim of emphasizing the potential importance of stochastic mechanisms in producing natural patterns (Gotelli and Graves 1996). The null model was built by unweighted randomization of the rows and columns (n × m) of a matrix of individual sizes arranged by species (n) and samples (m). For each randomization, an average value was produced for each sample (or sampling horizon), and after 9999 iterations, a grand mean and 2.5% and 97.5% percentiles were computed for each sample. Differences in mean body size were considered significant when confidence intervals did not overlap. This procedure was performed for each locality. All the analyses were performed using R v. 3.6.0 (R Core Team 2019). R code for the null model is described in the Supplementary Text.

Taxonomic Richness
A total of 2247 individuals belonging to 39 species and 29 genera of bivalves were sampled in this study (see Supplementary Dataset), with an average of 15 (±2) species per stratigraphic unit. Lowest richness was recorded in the Cotham Member and in the Angulata Zone, with 9 and 10 species respectively (Fig. 3A,B). The first drop in richness coincides with the extinction horizon (Cotham Member; Fig. 2B). In contrast, highest richness was recorded immediately after the extinction event (Langport Member); here, the taxonomic richness of bivalves increased by 90% (19 species) in less than ∼100 kyr . Species richness in the Westbury Formation and the pre-Planorbis, Planorbis, and Liasicus Zones fluctuated between 17 and 18 species, with no significant differences (Fig. 3A,B).

Body-Size Distribution across Time and Space
The BSD of all bivalve species (n = 39) fits to a lognormal distribution (r 2 = 0.87, p < 0.05) (Fig. 3C), with a modal size ∼8 mm, significantly right-skewed with size distributions that range from minimum and maximum sizes of 3.40 and 30.8 mm, respectively. The relationship between species richness and body size (Fig. 3D) exhibits a platykurtic distribution during the extinction event (Cotham Member). Assemblages from the immediate aftermath of the LTE (Langport Member) record a leptokurtic distribution with a modal size of ∼3.75 log 2 mm (Fig. 3D). Most assemblages from the pre-extinction and recovery intervals record mesokurtic BSDs (mode ∼3.1 log 2 mm), except for the pre-Planorbis Zone, which exhibits a bimodal distribution (Fig. 3D). The Angulata Zone recorded a low species number but reports a higher mode (∼4 log 2 mm) and mesokurtic shape (Fig. 3D).
Bivalve assemblages from the Cotham Member record the lowest median size of all stratigraphic units in this study (Fig. 4A). Despite this, the Kruskal-Wallis test did not detect significant differences in median body size across the extinction event (χ 2 = 5.51, p > 0.48; Fig. 4A). However, the dispersal parameters and shape of the BSD change significantly through the Triassic-Jurassic interval ( Fig. 4B-D).
In the Westbury Formation, the BSD of bivalves is centered around the grand mean and lies within the confidence intervals, despite its slightly right-skewed distribution (Fig. 4D, Table 1). However, during the Cotham Member, the coefficient of variation significantly increases ( Fig. 4B) due to the occurrence of some smaller species that skew the distribution toward negative values ( Fig. 4C,D, Table 1). The Langport Member assemblage records the highest leptokurtic BSD, with little size variation and a significant negative skew ( Fig. 4B, D, Table 1, Supplementary Table 1), due to all species being clustered around the grand mean size in the extinction aftermath. The pre-Planorbis Zone records an expansion of the size spectrum, shown by an increase in the coefficient of variation, due to the occurrence of relatively larger taxa, but it retains a negative skew due to the overall high dominance of small taxa in the assemblage (Fig. 4, Table 1). The Planorbis and Liasicus Zones record a further significant expansion in variance, due to the high occurrence of even larger taxa in the assemblages ( Fig. 4B,C, and Table 1). The Angulata Zone records the largest average sizes with a significant right-skewed distribution and a significant reduction in the coefficient of variation (Fig. 4, Table 1). Analysis of the homogeneity of body-size variance showed no significant differences between localities, indicating that the BSDs of bivalves expanded or contracted synchronously at the different study sites (Supplementary Table 2).

Selectivity
The coefficient of variation in the BSDs of victims and survivors is slightly higher in the Cotham Member than in the underlying Westbury Formation (∼62%; Fig. 5A). Subsequently, during the Langport Member, variation of victims, survivors, and colonizers decreased (∼40%). During the overlying Blue Lias Formation, the colonizers and survivors record gradual increases in their dispersion to a peak in the Liasicus Zone (∼83%), followed by a decrease of 37% in the Angulata Zone ( Fig 5A).
The global extinction pulse, which occurred between the Westbury Formation and the Cotham Member, was not size selective and removed a total of seven species of different FIGURE 3. Species richness and body-size patterns of bivalves through the Triassic-Jurassic of the UK sections. A, Raw species richness recorded for each stratigraphic unit and locality. B, Rarefied richness. Average values (±95% confidence intervals) of species richness estimated as sampling size increases. Significant differences were assumed if 95% confidence intervals did not overlap. WF, Westbury Formation; CM, Cotham Member; LM, Langport Member; PPZ, pre-Planorbis Zone, PZ, Planorbis Zone, LZ, Liasicus Zone, AZ, Angulata Zone. C, Body-size distribution of all bivalves recorded across the Triassic-Jurassic sections. A lognormal distribution was fit to the observed BSD. D, Bivariate relationships between body size and raw species richness plotted for each stratigraphic unit. Each line represents the kernel density fit to the data. sizes from along the BSD (Fig. 5B). During the Cotham Member, the average size of the surviving taxa (nine species) declined by ∼13%, with sizes ranging from X = 4.2 mm to X = 19.2 mm (Supplementary Table 1). Both Modiolus sp. and Mytilus cloacinus underwent significant decreases in mean size (by ∼54%), which increased the coefficient of variation of surviving taxa (Fig. 5B). The positive side of the residual distribution is occupied by larger taxa such as Cardinia regularis, Chlamys valoniensis, and Protocardia rhaetica (Fig. 5B). During the Cotham Member, three species disappeared. Of those, Placunopsis alpina and M. cloacinus record significant reductions in size (of ∼63%) before disappearing, while Rhaetavicula contorta did not record any change in size before its extinction. The first colonizer, Pleuromya sp. (∼19.2 ± 8.3 mm), is recorded in the Cotham Member assemblage (Fig. 5B).
The overlying Langport Member records the lowest variance in size and the first pulse of colonizers (Fig. 5A). Nineteen species were recorded in the Langport Member, of which four subsequently disappeared: Astarte sp. and Mesomiltha sp. represent single occurrences and local extinctions, while Liostrea bristovi and Isocyprina concentricum increased in mean body size by ∼34% before becoming globally extinct (Fig. 5B). Eighty percent of surviving taxa (Modiolus sp., P. rhaetica, C. valoniensis, and Pteromya crowcombeia) increased in mean size by >35% from the Cotham Member, while C. regularis recorded no change (Fig. 5B). The high numbers of colonizing species in the Langport Member assemblages and the lowest coefficient of variation recorded (∼40%) suggest that all species present in the aftermath of the LTE exhibited similar sizes (Fig. 5B).
Newly arriving/colonizing taxa gradually appear in the community, leading to an expansion in BSD toward both the positive and negative sides from the Langport Member to the Liasicus Zone (Fig. 5A,B). The first major expansion in body size happened between the Langport Member and the pre-Planorbis Zone, with 11 taxa contributing to a ∼54% expansion in size variance (Fig. 5A,B) and shifting the skewness toward positive values ( Fig  4D). The upper boundary of the BSD is driven by taxa such as Plagiostoma giganteum and Pholadomya sp., which record their first occurrences in the Langport Member, with mean sizes of 19.1 ± 5.2 mm and 19.1 ± 7.8 mm, respectively. Subsequently, both taxa expanded their body size by ∼66% and 80% into the pre-Planorbis FIGURE 4. Temporal changes in the dispersal parameters of the body-size distributions (BSDs) of 37 Triassic-Jurassic bivalve taxa recorded in this study. A, Residual values of the log 2 body size grouped by stratigraphic unit; each point represents one species, with each symbol indicating the locality from which it was sampled (St. Audrie's Bay; Larne, Northern Ireland; and Pinhay Bay). Log 2 mean body size ± SE is indicated for each stratigraphic unit. Zero on the y-axis represents the grand mean estimated from all species. Dashed lines are 95% confidence intervals estimated by bootstrapping with replacement (9999 iterations). B, Mean coefficient of variation of body size for each stratigraphic unit (-•-). Thin lines are 95% confidence intervals estimated by the bootstrapping procedure. C, Minimum (Min.), mean (Mean: -•-), and maximum (Max.) size for each stratigraphic unit. D, Skewness of the BSDs for each stratigraphic unit. Positive values indicate a right skew with a long tail of the data to the right; negative values indicate a left skew to the data. Size kernel-density plots (pale gray) show changes in the BSD between each stratigraphic unit. Key to stratigraphic units as in Fig. 3. and the Planorbis Zones, respectively. The lower boundary of the BSD is defined by the smallest taxa, such as Gervillella ornata (in the Langport Member), Modiolus minimus (in the pre-Planorbis Zone), and extremely small taxa such as Palaeoneilo elliptica (4.5 mm) and Oxytoma fallax (2.2 mm), which expanded the variance of the BSD by ∼34% (Fig. 5A,B) in the Planorbis Zone.

Phylogenetic Signal
Moran's autocorrelogram for body size shows a typical phylogenetic pattern that decays from generic to ordinal level (Fig. 5C), although the results were not statistically significant (Supplementary Table 3). This indicates that the variation in size is distributed relatively uniformly throughout the different taxonomic levels and there is no difference in body size between taxonomic levels. In addition, no signal of selectivity was observed when the ordinal size and species size of each ecological category were regressed (Fig. 5D, Table 2); only the relationship between average ordinal body size and the size of colonizers recorded a significant slope (Table 2).

Body-Size Distribution and Tiering
BSDs were highly variable between tiering categories and stratigraphic units, and only a third of BSDs were normally distributed ( Fig. 6A-G, Supplementary Table 4). The BSDs of infaunal and semi-infaunal bivalves differed significantly in all stratigraphic units (K-S test: d = 0.51, p c < 0.01; Supplementary Table 5). The BSDs of semi-infaunal and epifaunal taxa were not significantly different from each other during the pre-extinction Westbury Formation and later postextinction Angulata Zone (K-S test: d = 0.27, p c > 0.01), whereas the infaunal and epifaunal BSDs had the same shapes during the Cotham Member, pre-Planorbis Zone, and Angulata Zone (K-S test: d = 0.209, p c > 0.01) (Supplementary Table 5). Only during the Langport Member and Planorbis Zone did all tiers have significantly different BSDs from one another (K-S test: d = 0.42, p c < 0.01; Supplementary Table 5).
Dramatic changes in the shapes of the BSDs occurred during and soon after the extinction event (i.e., in the Cotham and Langport Members). For instance, between the Westbury Formation and Cotham Member, the BSDs of all tiers change from being slightly leptokurtic, right skewed, and not-normally distributed, to being mesokurtic with high dispersion (of ∼36%), while epifaunal and semi-infaunal  Table 4).
In the lower Blue Lias Formation, from the pre-Planorbis Zone to the Liasicus Zone, BSDs tend to be mesokurtic, with a switch toward positive skews and with a progressive expansion of variance in all three tiers (Table 1). The Angulata Zone, however, records a contraction in the BSDs of all tiers by around ∼18%. The BSDs of infaunal and epifaunal taxa record left-skewed distributions with high kurtosis (with mode ∼4 log 2 mm), while the BSD of semi-infaunal species is mesokurtic with a positive skew (Supplementary Table 4).
Significant differences were detected by GLMs between the mean body sizes of each tier (Wald χ 2 (2) = 53.29, p < 0.001). Epifaunal taxa recorded the largest mean size (15.7 ± 0.4 mm), infaunal taxa recorded intermediate mean sizes (12.4 ± 0.3 mm), and semi-infaunal bivalves recorded the smallest mean sizes (Fig. 6H, Supplementary Fig. 3). Although the GLMs did not find significant differences between the mean sizes of all bivalves from the Cotham Member (10.3 ± 0.2 mm) and Westbury Formation (8.5 ± 0.4 mm), both of those stratigraphic units record the smallest mean sizes in each section studied (Wald χ 2 (6) = 52.93, p < 0.05) ( Supplementary Fig. 3, Supplementary Table 6). Bivalves from the Langport Member recorded a significant increase in size of ∼30% (13.7 ± 0.2 mm) compared with those from the Cotham Member and Westbury Formation. Mean body size gradually increases from the Rhaetian Langport Member through the overlying Hettangian Blue Lias Formation, reaching an average size of 15.7 ± 0.3 mm in the Angulata Zone ( Supplementary Fig. 3, Supplementary Table 6).
The mean body size of each mode of life significantly increases through the Triassic-Jurassic interval (Wald χ 2 (12) =153.96, p < 0.05) (Fig. 6H, Supplementary Table 6). However, during the extinction event, semi-infaunal bivalves are the only tier to record a significant decrease in mean body size; from 10.4 ± 0.4 mm in the Westbury Formation to 4.0 ± 0.4 mm in the Cotham Member. This decrease is probably due to there being only a single taxon, Modiolus sp., representing this mode of life in the extinction horizon.

Comparison with the Null Model
Comparison of the mean body size of each species per sample with the average curves from the null model shows that the bivalves record a significant directional trend toward larger body size through the study interval (Fig. 7, Table 3). This pattern is broadly consistent across the region, although there appears to be some variation at the local scale during the recovery: mean size in the St. Audrie's Bay section increases significantly in the late Liasicus Zone (35 m above base of measured section); at Larne, larger sizes are recorded from the early Liasicus Zone (at 45 m); while at Pinhay Bay, a significant expansion of body size is recorded during the pre-Planorbis Zone (at 7.5 m). In all cases, the expansion toward larger  Fig. 7A-C). Finally, although sizes were not significantly smaller than expected by chance during the extinction event, the recurrent overlap of the randomized and sampled data during the Cotham and Langport members does indicates a high occurrence of smallersized taxa during those units.

Patterns of Change in Body Size during the Extinction and Recovery Event
Results indicate that the LTE simultaneously affected different localities, within the temporal resolution of the stratigraphic units analyzed. The BSDs of recorded bivalve species changed through the event and recovery, with a flattening of the BSD immediately after the LTE (Fig. 5A,B), followed by a narrowing of BSD and later rapid expansion of the body-size range during the recovery phase. The patterns of body-size changes recorded in this study resemble those recorded in studies of other extinction events (Payne 2005 (2020) reported changes in species composition between the Westbury Formation and Cotham Member (Lilstock Formation) and recorded an apparent reduction of body size of new arrivals and surviving taxa during the early recovery interval (i.e., the Langport Member). This study confirms the same changes during the extinction interval, although the changes in body size are complex.
The first extinction pulse, between the Westbury Formation and Cotham Member, was not size selective and so generated gaps across the BSD (Fig. 5B). In addition, the BSD of bivalve species in the Cotham Member has a negative skew due to the presence of surviving taxa that record a decrease in body size below the grand mean. These taxa include Isocyprina concentricum, Modiolus sp., Placunopsis alpina, and Mytilus cloacinus (Fig. 5B). Immediately after the extinction event (i.e., in the Langport Member), the variance of the BSD decreases significantly, because all recorded taxa cluster around the grand mean (Figs. 4B-D and 5A). In the overlying Blue Lias Formation, an expansion in size is recorded from the pre-Planorbis Zone, and the BSDs shift toward a positive skew due the occurrence of larger taxa during the recovery (Figs. 4D, 5A,B).
Similar increases in body size during the Early Jurassic recovery are also recorded from other studies in the United Kingdom and elsewhere. The trace fossils Thalassinoides, Diplocraterion, and Arenicolites show a substantial size increase from the Langport Member to the Angulata Zone in the southwest United Kingdom (Twitchett and Barras 2004;Twitchett 2007, 2016). Ammonites that survived the LTE expanded their body size by more than 10% toward the top of the Hettangian (Dommergues et al. 2002). Studies at community level from sites within the United Kingdom confirm a progressive accumulation of larger bivalve taxa through the Blue Lias Formation , and recent studies (Atkinson andWignall 2019, 2020) confirm our observations, demonstrating that body-size increase is driven mainly by newly originated clades.
Similar size increases in marine animals during postextinction recovery intervals have been recorded for other extinction events. Following the end-Devonian mass extinction, crinoids FIGURE 6. Changes in the size-frequency distributions of the body sizes of infaunal, semi-infaunal, and epifaunal bivalves through the Triassic-Jurassic of the United Kingdom. A, Westbury Formation (WF); B, Cotham Member (CM; Lilstock Formation); C, Langport Member (LM; Lilstock Formation); D, pre-Planorbis Zone (PPZ); E, Planorbis Zone (PZ); F, Liasicus Zone (LZ); G, Angulata Zone (AZ). The fitting of each curve is based on a combination of kernel density estimation and smoothed bootstrap resampling using bin sizes of 0.25 on a log 2 scale. All distributions are unimodal but with highly variable skewed frequency distributions. H, Changes in the weighted marginal geometric mean values (±SE) for each tier through time. exhibited significant body-size increase (Brom et al. 2018). Conodonts exhibit significant increase in size and high disparity level after the Frasnian-Famennian crisis (Girard and Renaud 2012). Brachiopods, bivalves, gastropods, and trace fossils show a directional trend to larger sizes during the recovery after the late Permian mass extinction (Twitchett 2007;Metcalfe et al. 2011;Pietsch et al. 2016;Foster et al. 2017Foster et al. , 2018Foster et al. , 2020. Similar results were recorded by Rego et al. (2012) with foraminifera from the late Permian to the Late Triassic. Cephalopods and bivalves recorded an increase in body size after the early Toarcian extinction event (Morten and Twitchett 2009;Rita et al. 2019), and venerid bivalves recorded increases in body size after the Cretaceous/ Paleogene (K/Pg) boundary (Lockwood 2005).
The temporal trajectory of the body size of each species, expressed as residuals in this study, resembles the diffusion model in which all species move away from the lower bound or "small mean size" (see Figs. 1, 5B). This could represent a physiological optimum (Brown et al. 1993;McShea 1994) wherein macroevolutionary forces such as extinction risk and selective advantage determine the length of the right and the left tail of the BSD, respectively (Clauset and Erwin 2008), and constrain the upper and lower limits of a taxon's body size around a physiological optimum size.
Testing the Lilliput Effect Model Twitchett (2007) and Harries and Knorr (2009) described four potential and nonexclusive scenarios to explain why small-sized animals dominate fossil assemblages after extinction events: (1) the preferential survival of small taxa and/or selective extinction of large taxa; (2) the postextinction origination of many newly evolved lineages at small initial body size; (3) a potential bias in the preservation of larger taxa in the immediate extinction aftermath; and (4) the Lilliput effect sensu stricto of Urbanek (1993), whereby surviving taxa undergo temporary size reduction in the immediate aftermath of the extinction event.
In the case of Triassic-Jurassic bivalves, our study shows that there was no size selection during the LTE, at least in the UK sections we studied. It appears that the LTE did not selectively eliminate larger or smaller taxa; rather, each body-size class had the same probability of extinction (Fig. 5B). During the extinction event (Cotham Member), the mean body size of the entire community did not change significantly compared with the pre-extinction Westbury Formation. However, the BSD remained highly platykurtic, because taxa were extirpated randomly with respect to body size, which resulted in a redistribution of species along the size axis (Fig. 5B). For instance, M. cloacinus, P. alpina, Modiolus sp., and I. concentricum record apparent dwarfism (i.e., the Lilliput Effect sensu stricto), three taxa record size increases (Pteromya crowcombeia, Cardinia regularis, and Chlamys valoniensis), two species do not record any changes (Rhaetavicula contorta and Protocardia rhaetica), and only one taxon (Pleuromya sp.) emerges at the upper side of the BSD (Fig. 5B).
Langport Member assemblages record a highly leptokurtic BSD with lower coefficient of variation immediately after the extinction, because the sizes of the survivors and victims FIGURE 7. Contrasting null models of Triassic-Jurassic bivalve body size from St. Audrie's Bay (A), Larne (B), and Pinhay Bay (C). Each point represents the mean size of a specific genus recorded at the stratigraphic height shown (see Supplementary Data for detailed stratigraphic logs). Genus-level body sizes were plotted to avoid overcrowding of data points. Thick black lines show the average size of all individual bivalves from each sample, estimated from raw data, with 2.5% and 97.5% percentiles. The gray shaded area with thin black lines represents the null model, with 2.5% and 97.5% percentiles, calculated by row permutation (number of iterations = 9999) of the geometric mean of each individual by species through the samples. Key to stratigraphic units as in Fig. 3. were centered around the grand mean size, with a skew close to zero (see Figs. 4D, 5A). The clustering of species around the mean size immediately after the LTE (Langport Member) is influenced by the arrival of comparatively smaller species belonging to new clades, which could indicate the presence of environmental constraints on size. Subsequently, newly arriving taxa expand through morphospace by increasing their size, and larger taxa become more common, which extends the right tail of the BSD. Thus, the size changes recorded in bivalves through the LTE in the UK sections studied appear to be caused by a combination of the Lilliput effect sensu stricto and the appearance of new, small-sized taxa in the extinction aftermath, and not by selective extinction. The BSD patterns observed during the Cotham and Langport members are also consistent with McKinney's (1990) suggestion that after a nonselective extinction event, the BSD changes to a flatter version of the pre-extinction BSD as result of an equal chance of extinction of each size class. Each size bin had the same probability of extinction, and so tails of the BSD may be more susceptible to disappear, because they contain proportionally fewer taxa. Finally, even though the mechanisms of body-size reduction were not evaluated here, several physiological responses to cope with warming events have been identified in some organisms. For instance, shifts in cell sizes, reduction of the metabolic rate, growth, and/ or modification of life cycles (Atkinson 1994;Millien et al. 2006;Daufresne et al. 2009;Forster et al. 2012;Hermaniuk et al. 2021;Verberk et al. 2021). Recently, Rita et al. (2019) suggested that the high occurrence of juvenile belemnites during the Pliensbachian-Toarcian extinction event could be an ecophysiological response to cope with high temperatures. The faster maturity at smaller sizes could reduce the resource-dependent mortality and intraspecific competition in warm conditions (Sheridan and Bickford 2011;Forster et al. 2012;Le Bris et al. 2017;Verberk et al. 2020Verberk et al. , 2021. In this sense, we could expect a similar ecological signal for bivalves across the end-Triassic extinction, considering the similarity of the environmental drivers. However, further data are required to unravel the mechanisms that drove size changes through the end-Triassic hyperthermal event.

Potential Abiotic Factors and Timing of Recovery
Two questions arise regarding the trajectory of the body-size changes: (1) What drove the changes in bivalve body size during the recovery stage? (2) How does the rate of recovery recorded in our study compare with results from previous studies of the LTE? Regarding the potential drivers of size change, a number of studies have identified abiotic factors that could have driven body-size changes in marine organisms during past extinction events (Twitchett 2007;Zhang and He 2008;Wade and Twitchett 2009;He et al. 2016) and in relation to present-day climate change (e.g., Daufresne et al. 2009;Stige and Kvile 2017). The dissolved oxygen content of aquatic systems has been invoked as one of the main factors in structuring body size and species diversity at different temporal scales (Payne et al. 2011;Riding et al. 2019;Rubalcaba et al. 2020).
Increases in oxygen levels through the Phanerozoic have had a strong control upon the ecological complexity of marine systems and the body size of marine taxa (Novack-Gottshall 2008; Payne et al. 2011), especially during environmental stress. For example, by using the relative abundance of the rare earth element cerium, which is a proxy for paleo-oxygenation of marine environments, He et al. (2010) demonstrated that during the late Permian mass extinction event, the body size of brachiopods decreased strongly with increasing hypoxia. Similar size responses have been described for other taxa during extinction events that are associated with oceanic anoxic episodes: for example, Pliensbachian-Toarcian ammonites (Morten and Twitchett 2009) and belemnites (Rita et al. 2019) and Permian-Triassic gastropods and bivalves (Metcalfe et al. 2011). The burrow diameters of infaunal benthic marine macroinvertebrates also record similar changes through many past extinction events (Twitchett and Barras 2004), including through the LTE (Barras and Twitchett 2007).
The LTE was ultimately caused by a major episode of global warming due to CAMP eruptions, with a rise in atmospheric CO 2 from ca. 750 to ca. 2000 ppmv and an increase of 3°C -4°C in mean seawater temperature (McElwain et al. 1999;Greene et al. 2012). Widespread hypoxia and anoxia are a predicted consequence of global warming due to lower dissolved oxygen concentrations in warmer seawater, weakness of thermohaline ocean circulation, increased likelihood of stratification, and shallowing of euxinic conditions (Kidder and Worsley 2010;Breitburg et al. 2018). Jaraula et al. (2013) used biomarker evidence (e.g., the abundance of isorenieratane) to suggest that fluctuations between dysoxic, anoxic, and euxinic conditions could be one of the main factors constraining the post-LTE recovery phase. Oxygen concentration is only one of the major marine biogeochemical parameters predicted to change during warming events (Rubalcaba et al. 2020) but has the potential to trigger drastic physiological modifications in marine ectotherms (Sheridan and Bickford 2011;Rubalcaba et al. 2020;Verberk et al. 2021). Regarding the timescale, our study indicates that the extinction episode and the appearance of the first colonizers was faster than previously proposed (Thibodeau et al. 2016). However, recent studies performed by Weedon et al. (2019) would assign to the Hettangian a longer duration (4.1 Myr) compared with previous observations (1.7-2 Myr) (Schaltegger et al. 2008;Schoene et al. 2010;Wotzlaw et al. 2014) and cyclostratigraphy (Ruhl et al. , 2016Hüsing et al. 2014). This would suggest that recovery during the Hettangian was slower than previously indicated (Thibodeau et al. 2016;Lindström et al. 2017), whereby the ecological impact of the LTE requires a reevaluation.
The most current chronology of the biotic and abiotic events associated with the LTE establishes that CAMP activity began around the last occurrence of Choristoceras marshi (∼201.52 ± 0.14 Ma) (Lindström et al. 2017), which coincides with the initial negative δ 13 C org excursion, which had a duration of ∼20-40 kyr . Thibodeau et al. (2016) divides the composition changes of the marine fauna with respect to CAMP activity into two stages: first, the "extinction interval" that spans from C. crickmayi to the first occurrence of Psiloceras spelae, that is, the base of the Jurassic; and second, the "depauperate interval" that spans from P. spelae to the Kammerkarites Zone. In the UK succession, the first stage would equate to the interval spanning the Cotham Member to the Triassic/Jurassic boundary (i.e., in the pre-Planorbis Zone), and the second stage would span from the base of the Hettangian to the base of the Liasicus Zone (Fig. 2C).
Our data show that newly colonizing taxa begin to arrive during the Langport Member, and the BSDs increase significantly in variance in all localities (Fig. 5) before the Triassic/Jurassic boundary. Using the timescale of Ruhl et al. (2010) and , this suggests that initial recovery in the United Kingdom had already begun by ca. 100-120 kyr after the extinction, that is, before peak warming and the depauperate interval of Thibodeau et al. (2016) and while CAMP volcanism was still active. Several studies have shown that after the LTE, marine communities along the western coast of Pangea had lower diversity than those in Tethys (Damborenea et al. 2013;Ritterbush et al. 2014), although these apparent differences in recovery rates have not been explored in detail.
Our results agree with previous studies that have suggested that final recovery in the United Kingdom is reached in the late Hettangian (e.g., Barras and Twitchett 2007;Pugh et al. 2014;Atkinson andWignall 2019, 2020;. Although not sampled in this study, previous analyses have also demonstrated that the sizes of benthic marine invertebrates continued to increase into the early Sinemurian (e.g., Barras and Twitchett 2007). The absolute rates of recovery and size change depend on estimates of the duration of the Hettangian stage. Most studies have concluded that the Hettangian lasted between 1.7 and 2 Myr (e.g., Ruhl et al. 2010Ruhl et al. , 2016Hüsing et al. 2014;Gradstein et al. 2020); although Weedon et al. (2019) suggest a longer duration of 4.1 Myr. If correct, this would suggest that recovery during the Hettangian was significantly slower in absolute terms than previously suggested (Thibodeau et al. 2016;Lindström et al. 2017).

Selectivity
Although there is no clearcut phylogenetic signal associated with the body sizes of survivors or extirpated taxa, our data show that semi-infaunal bivalves were most strongly affected during the extinction event (recorded as a decrease in body size) (Fig. 6H, Supplementary Tables S4, S6). These results suggest that ecology was more important in determining the body size in surviving taxa than historical (phylogenetic) aspects. The lack of a phylogenetic signal in the sizes of surviving or extirpated taxa is consistent with the premise that there is no size-related selectivity during mass extinction, meaning that the loss of taxa from any hierarchical level is no different from what would be expected by chance (Jablonski and Raup 1995;McRoberts and Newton 1995). Variability in size-selective extinction has been reported commonly in bivalves (Lockwood 2005;Smith and Roy 2006;Rivadeneira and Marquet 2007); however, the lack a phylogenetic signal in our data could be explained by the low variance in size between taxonomic levels and/or simply because the LTE was not selective at different extirpation categories and taxonomic levels.
Selective extinction against infaunal taxa has also been reported from previous studies of the LTE (Kiessling et al. 2007) and other extinction events (Jablonski and Raup 1995;Rivadeneira and Marquet 2007), and this study found that 55% (6 out of 11) of the infaunal taxa went extinct during the LTE compared with 43% (3 out of 7) of the epifaunal forms. However, as demonstrated by , infaunal suspension feeders are relatively poorly preserved through the extinction interval, particularly in the southwest United Kingdom, and so rock record bias may explain part of this apparent selectivity . When we decomposed the infaunal group into semi-infaunal and shallow-infaunal forms, we found that semiinfaunal bivalves appear to be more strongly affected by the extinction event (80%, 4 out of 5 genera), with Modiolus sp. the only survivor. This may explain why semi-infaunal bivalves record a platykurtic distribution and smallest size during the Cotham Member (Fig. 5B).
Of the 13 taxa that disappear between Cotham Member and Langport Member, 5 have completely aragonitic skeletons, 7 are aragonitic with low Mg-calcite, and just 1 taxon has purely calcitic skeletal mineralogy. There is, however, a strong link between life mode and mineralogy, with almost all infaunal and semi-infaunal bivalves recorded in the Triassic-Jurassic of the United Kingdom having aragonitic skeletons . While this likely explains the taphonomic bias in the infaunal record, it also means that it is not possible to determine whether mineralogy or life habit was most important in determining selectivity.
Several studies have hypothesized that the LTE was accompanied by a biocalcification crisis due to ocean acidification, as evidenced by the extinction of hypercalcified taxa, occurrences of skeletal structures with malformations, decreases in shell thickness, and the high abundance of noncalcareous disaster taxa with organic-walled structures (Kiessling et al. 2007;Hautmann et al. 2008;Crne et al. 2011). However, in their multivariate analysis of all available global data, Dunhill et al. (2018) recently showed that calcification is not a strong predictor of extinction selectivity during the LTE and concluded that ocean acidification was not a major factor in the extinctions. Indeed, benthic mollusks are able to survive in surprisingly acidic conditions, and their shells record characteristic dissolution and repair features when affected by in vivo acidification (Garilli et al. 2015). The presence of such features in fossil mollusks would provide unequivocal evidence of acidified conditions, but thus far has not been reported from the Triassic-Jurassic record.

Taphonomic Effects and Biases within the Dataset
Taphonomic biases are known to affect the fossil record, and it is possible that changes in BSDs through time could be affected by several potential biases, such as lithology, preservation, and changing water depth (Behrensmeyer et al. 2000;Holland 2000;Hunt 2004).  documented in detail the preservational and facies changes affecting the bivalve record of the southwest UK Triassic-Jurassic successions. They also undertook a gap analysis at the genus level to assess the relative completeness of the same stratigraphic units that were sampled in this study and concluded that, in most cases, between 54% and 61% of the genera expected in each unit actually occurred as fossils, and only the pre-Planorbis Zone was substantially worse (recording only 41% of expected genera; . In their analysis of different tiers, Mander and  indicates that infaunal suspension feeders have a relatively poor fossil record overall, but particularly through the extinction and in the pre-Planorbis Zone. In addition, epifaunal bivalves are relatively poorly represented in the Blue Lias Formation compared with the underlying Westbury and Lilstock formations; and semi-infaunal bivalves are better represented from the Planorbis Zone onward than in the underlying stratigraphic units . Certainly, the pre-Planorbis Zone does not appear to be anomalous, in terms of body size, or in any of our comparative analyses of the entire fauna (Figs. 2-4, Supplementary  Fig. S4). Indeed, if preservation acted as a significant factor constraining our observations, then we would not expect to have recorded BSDs with similar shapes to those reported in other marine communities (cf. Warwick 1984;Labra et al. 2015) (see Fig. 3D). In addition, our data have filled in some of the gaps in ranges of several genera reported by . Thus, we are confident that our size data are robust, with respect to taphonomy, and the differences found between the null model and the raw body-size data through the Triassic-Jurassic (Fig. 6) capture a true ecological signal and are not artifacts of preservation bias.
In terms of changing lithofacies, the Cotham Member is much sandier than the other units, representing a much shallower depositional setting, and the Langport Member is more carbonate dominated than other units, representing a lagoonal or ramp setting (Hesselbo et al. 2004;. Although it is unlikely that lithofacies is biasing preservation, and hence diversity, through most of the study interval, it is possible that the depositional setting controls recorded body size in the Lilstock Formation. The relatively small size of the Cotham Member assemblage may be at least partially explained by the environment if, for example, salinity was lower due to greater freshwater influence in those more proximal settings (Wignall 2001;Hesselbo et al. 2004;Warrington et al. 2008). If salinity was normal, then we might expect bivalve species in shallower marine settings to be larger (Hauton 2016). The presence of localized, small, patchy stromatolites within the Cotham Member at some locations in southwest United Kingdom, though not those sampled in this study, also likely indicates abnormal, schizohaline conditions or may reflect the postextinction expansion of stromatolites (Radley et al. 2008). Evidence of fluctuating salinities has also been recorded at some horizons in the overlying Langport Member (Boomer et al. 1999). Thus, the small size of Cotham Member bivalves (Figs. 4D, 5B) and the low coefficient of variation in Langport Member assemblages (Figs. 4B, 5A) seem to be due to an environmental (salinity?) control.
Relative water depths also vary between the study sites. The sections at Pinhay Bay, St. Audrie's Bay, and Larne represent a depth gradient, from inshore to offshore (Hesselbo et al. 2004). Despite this gradient, however, changes in BSDs happen simultaneously at the different localities (Fig. 4A), suggesting that relative position on the shelf does not significantly affect the timing of recorded bodysize changes and overall trends.
Local and global sea levels, and hence water depth, did not remain constant through the studied interval. If fossil richness and abundance of the fossil record are correlated with sea level Holland 2016), then covariation between body size and sea level could be expected. High-quality fossils (abundant and well preserved) would be expected during transgressive to earlyhighstand phases (e.g., the lower Westbury Formation, Langport Member, and lower pre-Planorbis Zone; Hesselbo et al. 2004;. Moderate quality fossils would be expected during the highstand phases (upper Westbury Formation), and poor fossil quality would likely be observed during latehighstand and regressive phases (i.e., the Cotham Member; . The raw species richness and rarefaction analyses of the Cotham Member are potentially consistent with this (Fig. 3A,B). The significant increase in the mean body size of the bivalve community from the Langport Member to the Angulata Zone apparently occurred without significant change in either sea level or lithofacies ( Supplementary Fig. S4), suggesting that other factors were shaping the BSD during the postextinction interval. Smaller sizes were recorded in the Westbury Formation and, in particular, in the Cotham Member ( Supplementary  Fig. S4), and a narrowing of the BSD was recorded during transgressive to earlyhighstand phases (Langport Member) (Fig. 3D, Supplementary Fig. S4).
Alternatively, the nature of the preservation processes, as time averaging and/or transport, could affect the shape of BSD (Fagerstrom 1964;Mancini 1978). Preburial transportation of loose shells by strong currents, could constrain the BSD toward the accumulation of larger sizes or generation of leptokurtic distributions (Fagerstrom 1964), as we see in the Langport Member, for instance. On the other hand, increased degree of time averaging/condensation could lead to a more platykurtic BSD (e.g., Hunt 2004), as was recorded at the Cotham Member. However, the relatively lowenergy Bahamian type that prevailed during the deposition of the Cotham Member and the Langport Member (Hallam 1960;Radley 2002;Hesselbo et al. 2004) preclude the accumulation of large size and/or development of leptokurtic BSD in bivalves. Alternatively, the platykurtic BSD observed in the Cotham Member resulted in very short period of time (∼20 kyr).  in the disappearance of species and the size expansion-reduction of some survivor taxa, whereby the potential effect of time averaging does not seem to fit from the pattern observed in the Cotham Member. Thus, the pattern of change in the BSD recorded in this study represents a true paleoecological signal.
One clear sampling bias is apparent in our data from the Angulata Zone, which are evidently undersampled (Fig. 3B, Supplementary  Fig. S4). The lowest richness and lowest variance in body size recorded from this stratigraphic unit could be related to the relatively low number of samples compared with other units. However, despite that bias, this study still records similar trends in bivalve body size as those recorded by . Thus, the sustained increase in body size recorded after the extinction event in the whole assemblage and in individual taxa is unlikely to be the result of sampling or taphonomic biases, especially given the nature of the Blue Lias Formation rock record  and the fact that similar increases are recorded in other taxa, such as trace fossils (Barras and Twitchett 2007).
By sampling three locations in one region (northwest Europe), this study better represents the true changes in body size through the Triassic-Jurassic extinction and recovery interval than previous studies that have relied on single sites, and it also reduces the effects of variation between localities, species, and stratigraphic units. As noted by Bambach (1977) and Novack-Gottshall (2008), comparison of size trends among phylogenetically and ecologically unrelated taxa living in the same habitat is an important and powerful test of ecosystem-level processes, because environmental changes will often affect entire communities of organisms in a similar manner. Our data show that taphonomic heterogeneity does not strongly affect the BSD patterns found here. Despite this, it is noted that random effects can be more influential when sample sizes are low (Novack-Gottshall 2008).
Our study agrees with previous empirical and theoretical studies of body size (McKinney 1990(McKinney , 1997 and provides new knowledge of the dynamics of body-size change through the LTE and subsequent Early Jurassic recovery. Our results suggest that the LTE acted rapidly and synchronously on the marine bivalve fauna across multiple sites, decreasing taxic richness and strongly modifying the shape of body-size distributions. No significant phylogenetic, ecological, and/or size-selective signals were found. This nonselective extinction resulted in a negatively skewed and platykurtic BSD in the shallow-marine, nearshore Cotham Member, probably controlled at least in part by the local environment. In the overlying Langport Member, the BSD still shows a negative skew, with the body sizes of surviving and extinct taxa accumulated around the median size, recording a low coefficient of variation. Newly colonizing taxa were also small, with body sizes near that of the median ("rebound size"). Thereafter, recovery of body size was relatively rapid, and later in the Hettangian, when the environmental conditions improved, larger individuals were recorded, and so the right tail of the BSD extended significantly.