Bayesian analyses indicate bivalves did not drive the downfall of brachiopods following the Permian-Triassic mass extinction

Domain: earth_science
Published: 2023-09-09
Editorial outcome: minor_revision

Abstract
Certain times of major biotic replacement have often been interpreted as broadly competitive, mediated by innovation in the succeeding clades. A classic example was the switch from brachiopods to bivalves as major seabed organisms following the Permian-Triassic mass extinction (PTME), ~252 million years ago. This was attributed to competitive exclusion of brachiopods by the better adapted bivalves or simply to the fact that brachiopods had been hit especially hard by the PTME. The brachiopod-bivalve switch is emblematic of the global turnover of marine faunas from Palaeozoic-type to Modern-type triggered by the PTME. Here, using Bayesian analyses, we find that unexpectedly the two clades displayed similar large-scale trends of diversification before the Jurassic. Insight from a multivariate birth-death model shows that the extinction of major brachiopod clades during the PTME set the stage for the brachiopod-bivalve switch, with differential responses to high ocean temperatures post-extinction further facilitating their displacement by bivalves. Our study strengthens evidence that brachiopods and bivalves were not competitors over macroevolutionary timescales, with extinction events and environmental stresses shaping their divergent fates.

Introduction
The Phanerozoic evolution of animal diversity has been shaped by the balance of extinction and origination processes underpinned by biotic factors and abiotic drivers, the Red Queen versus Court Jester models, respectively. The interplay of abiotic and biotic factors over different temporal and spatial scales caused the waxing and waning of clades and their replacement by one another. Of these, the decline of brachiopods, coupled with the rise of bivalves in the Phanerozoic is a textbook example of clade replacement in palaeontology (Fig. 1). Both clades originated in the Cambrian, survived the ‘Big Five’ mass extinctions, and thrive in today’s oceans. They have lived in the ocean for >500 million years (Myr) and are two of the most diversified invertebrate clades in marine ecosystems. Nevertheless, brachiopod diversity declined dramatically at the point of the Permian-Triassic mass extinction (PTME) ~252 Myr ago, whereas bivalve diversity increased, demonstrating the brachiopod-bivalve switch in their relative richness (Fig. 1).

The PTME coinciding with the brachiopod-bivalve switch also marks one of the largest events in the history of marine life, the switch from Palaeozoic- to Modern-type evolutionary marine faunas, the most dramatic turning point in the 540 Myr of the Phanerozoic. The rise of the modern marine fauna, dominated by bivalves, gastropods, crustaceans, echinoids, and neopterygian fishes, was also a point of substantial increase in energy capture by marine life mediated by new photosynthesising plankton groups from the Late Triassic onwards. The Mesozoic marine revolution began in the Late Triassic and Jurassic, represented by enhanced arms races between predators and prey and, in turn, the dominant new groups were faster and meatier than their Palaeozoic precursors.

Owing to similarities in body plan (i.e., two shells) and overlaps in ecology (i.e., similar feeding behaviours, modes of life, and living habitats), brachiopods and bivalves have long been regarded as competitors. Gould and Calloway first challenged this view and considered that both groups possessed comparable diversity patterns in the Palaeozoic, and that the PTME reset their initial diversities and altered their relative dominance. In arguing against active competition between members of both clades, they were described as ‘ships that pass in the night’. However, debate on the competition between the two clades and their ecological roles in ecosystems continues. Moreover, the decline of brachiopods after the PTME has also been attributed to intense predation, decreased ability to expand habitat distribution, or increased energy flux. Nevertheless, data supporting the previous scenarios mostly relied on the relative richnesses of brachiopods and bivalves, without considering their underlying diversification dynamics (i.e., origination and extinction rates), nor have their relationships with biotic (e.g., self- and predator diversity dependence) and abiotic (e.g., temperature, seawater chemistry, tectonic regulation) factors received substantial statistical scrutiny (but see Liow et al. and Reitan and Liow).

Here we quantify the diversification dynamics of brachiopods and bivalves, along with their drivers, in a Bayesian framework (PyRate) to investigate the timing and triggers of the brachiopod-bivalve switch. First, we estimate the origination and extinction rates, and diversity of the post-Cambrian brachiopods and bivalves using fossil occurrence data from the Paleobiology Database (PBDB). This analysis shows the long-term pattern of their diversification history and highlights the importance of the Permian–Jurassic interval in assessing the brachiopod-bivalve switch. We then exhaustively revise the taxonomy and stratigraphy of the global fossil record of Permian–Jurassic brachiopods and bivalves and use these bespoke datasets to examine the dynamics of the brachiopod-bivalve switch at high temporal resolution, with emphasis on differences between ecological guilds and geographical regions. Finally, we use a multivariate birth-death model (PyRateMBD) to evaluate the potential factors driving the extinction and recovery of brachiopods and bivalves across the PTME.

Results
Diversification dynamics of brachiopods versus bivalves
The Bayesian analyses of post-Cambrian datasets show that the large-scale diversification dynamic patterns of brachiopods and bivalves were commonly shaped by major biotic and environmental events, such as the Great Ordovician Biodiversification Event, the ‘Big Five’ mass extinctions and their subsequent recovery phases (Fig. 2). A closer examination of the diversification rates, however, indicates very different patterns of the two clades before and after the Triassic-Jurassic boundary (i.e., largely similar trends from the Ordovician to Triassic vs. largely divergent trends from the Jurassic to recent). Despite different volatilities (which may have been resulted from different preservation rates), brachiopods and bivalves share a gradually decreasing trend in extinction and origination rates during the early to middle Palaeozoic, typical of rates of all marine animals. In the Jurassic, brachiopods experienced frequent extinction and origination events but the rates of bivalves were largely stable. After the Jurassic, brachiopods possessed almost stable extinction and origination rates except for a minor surge in extinction rate across the Cretaceous-Paleogene (K-Pg) boundary (Fig. 2d, e), while bivalves had much more volatile extinction and originating rates (Fig. 2a, b). The flat rates of the post-Jurassic brachiopods are unlikely to have been caused by their low diversity because the confidence interval is rather narrow, and PyRate has the ability to detect significant rate shifts in such a situation. Overall, the Bayesian analyses show comparable trends of diversification dynamics of pre-Jurassic brachiopods and bivalves, but less similar trends for post-Triassic taxa. A similar pattern is also found in analyses of the same datasets using the traditional method, i.e., the per-capita rate (Supplementary Figs. 39, 40, 43–48). The correlation between brachiopod and bivalve extinction rates is more prominent than that of their origination rates, indicating that their extinctions were generally caused by similar environmental events. The newly calculated diversities of post-Cambrian brachiopods and bivalves are largely comparable with those previously published and clearly demonstrate the brachiopod-bivalve switch across the PTME (Figs. 1, 2c, f).

Given the different evolutionary trends of brachiopods and bivalves before and after the Triassic-Jurassic boundary and the most pronounced diversity switch of brachiopods and bivalves occurring over the Permian–Triassic transition (Fig. 1), the newly revised datasets of Permian–Jurassic brachiopods and bivalves were analysed in detail. To enhance temporal resolution, the revised Permian–Jurassic datasets consist of stratigraphical occurrences that were assigned to the substage level when possible. The general trends derived from the emended datasets (Fig. 3) remain similar to those derived from the uncorrected data in the PBDB (Fig. 2) and largely agree with the rates estimated for discrete time bins (Supplementary Figs. 41, 42), but show more detailed fluctuations of rates. From the Early to Middle Permian, both brachiopods and bivalves showed decreasing origination rates, while brachiopod extinction rate exhibited a more pulsed, episodic pattern (Fig. 3a, d). The PTME severely impacted both groups, but the extinction rate of brachiopods was almost double that of bivalves (Fig. 3a, d). As a result, brachiopod diversity declined more sharply than that of bivalves across the PTME (Fig. 3a–f).

After the PTME, both groups underwent rapid diversifications in the first 10 Myr of the Triassic. More detailed comparisons indicate that bivalves recovered immediately after the PTME, while brachiopods rebounded ~2 Myr later, in the Olenekian (Fig. 3a, d). In the Ladinian, both groups displayed decreasing origination rates, followed by a pronounced surge in origination rate at the beginning of Late Triassic Period (Fig. 3a, d). Their origination and extinction rates then remained relatively consistent until the end of the Triassic, despite weak fluctuations, where both groups exhibited concordant spikes in origination and extinction rates during the Triassic-Jurassic mass extinction (TJME; Fig. 3a, d).

Jurassic bivalves and brachiopods displayed markedly different extinction and origination rates. Bivalves showed no shifts in origination rate, but their extinction rate curve peaked sharply at the Middle-Late Jurassic boundary and again at the end of the Jurassic (Fig. 3a). In contrast, brachiopods experienced a permanent increase in extinction rate during the Toarcian oceanic anoxic event and a further, albeit temporary spike near the Middle-Late Jurassic boundary (Fig. 3d). Although suffering high extinction rate, brachiopods still achieved a high origination rate in Middle Jurassic (Fig. 3d).

The decline of brachiopods and rise of bivalves coincide with the PTME. However, brachiopods suffered more, with several orders wiped out, whereas the main groups of bivalves were not eliminated. Consequently, it is somewhat inappropriate to compare the diversity of all Palaeozoic brachiopods with that of all Mesozoic bivalves. We therefore subdivided brachiopods into two major groups: the orders that went extinct in the PTME and the orders that survived or originated after the PTME (referred hereafter as PT and PT groups, respectively; see Methods for the detailed order names and assignments). Our results reveal that PT brachiopods did not show a strong trend of long-term diversity decline despite a sharp drop across PTME (Fig. 3i). Instead, their diversity curve frequently fluctuated, and by the Middle–Late Jurassic they returned to levels previously displayed in the Permian (Fig. 3i). Clearly, the superficial decline of the entire brachiopod phylum after the PTME is likely caused by the extinction of some diverse Palaeozoic groups, rather than a failure by the PT brachiopods to re-diversify. PT brachiopods did re-diversify in Middle–Late Triassic, and Early and Middle Jurassic, although their diversity remained lower than that of bivalves (Fig. 1).

Overall, the brachiopod-bivalve diversity replacement was neither the ‘double-wedge’ pattern (i.e., brachiopod diversity gradually declined, while bivalve diversity increased), nor the ‘mass extinction’ pattern (i.e., bivalves could only diversify after the extinction of brachiopods) (Fig. 1). Based on assumptions of ecological saturation and competition between brachiopods and bivalves, Sepkoski used the coupled logistic model to explain the observed diversity pattern of the two clades. That model assumes that the increase of bivalve diversity decreases the net (i.e., origination minus extinction) diversification rate of brachiopods, and the solution to the model requires a decline in brachiopod diversity in the middle–late Palaeozoic and a decline in net diversification rate of bivalves when approaching the present. Obviously, the newly calculated diversification rate patterns do not support it, nor do the high brachiopod diversity in the Permian (the peak is also prominent after sampling standardisation) or the rapid and non-stop increase in Cenozoic bivalve diversity. Conversely, Steele-Petrović hypothesised that the Mesozoic replacement of bivalves was achieved by a series of extinction events. Because of greater abilities of bivalves to resist extinction, to colonise, and to expand distribution, they rapidly invaded vacant habitats left by brachiopods after mass extinctions, and then prevented brachiopods from occupying these ecospaces, essentially as the ‘incumbent replacement model’. The superiority of bivalves as Steele-Petrović argued is shown by the fossil record: they rapidly ‘took over the world’ after the PTME, and had a much higher abundance than brachiopods in the earliest Triassic. However, the presence of bivalves did not prevent brachiopods from recovering in the Olenekian and Anisian as expected by the model (Fig. 3). In addition, the brachiopod radiation in the Middle Jurassic also indicates that brachiopod diversification was not prohibited by the bivalves. All the evidence suggests that brachiopods and bivalves had their own evolutionary paths. Neither the ‘passive replacement model’ nor the ‘active displacement model’ can sufficiently explain their diversification history.

Diversification dynamics between ecological groups
For brachiopods, there is a heterogeneity of rates between the different lifestyles (cemented, reclining, pedicle-attached, infaunal) (Supplementary Figs. 5–9). Infaunal brachiopods did not show significant rate shifts, but the other three groups displayed similar trends to one another in the Permian. Cemented and reclining taxa suffered a higher extinction rate during the PTME than pedicle-attached taxa, as most of the former belong to the PT group. Although a new cemented clade, Thecideida, arose in the Triassic, its fossil record is sparse and thus did not have a strong impact on the rates.

Bivalves have a greater number of ecological lifestyles than brachiopods. According to their mode of life and relative position of the animal and sediment, we categorised them into five groups: cemented, reclining, epibyssate (including semi-infaunal endobyssate; see Methods), shallow infaunal, and deep infaunal (Fig. 4, Supplementary Figs. 12–18). The large-scale patterns show that both epifaunal and infaunal bivalves possessed elevated extinction rates in the PTME and TJME, decreasing origination rates through the Permian, and elevated origination rates through Early Triassic to Carnian (Fig. 4a, c). Both groups rebounded rapidly from the PTME, with similar origination rates. Nevertheless, the extinction rate of epifaunal bivalves remained relatively high in the Triassic like brachiopods (Fig. 2d, g), while that of infaunal bivalves sharply dropped after the PTME, and remained unchanged except at the ends of the Triassic and Jurassic (Fig. 4a, c). After the Norian, epifaunal taxa maintained a constant origination rate (mean λ = 0.020), while infaunal taxa showed a slightly higher absolute value (mean λ = 0.027) coupled with a further origination episode after the TJME.

In the Triassic, the extinction rate of epifaunal bivalves was generally higher than that of infaunal bivalves, as demonstrated previously. In the Jurassic, the origination rate of epifaunal bivalves was low or very slightly higher than their extinction rate, while infaunal bivalves showed higher origination rates than extinction rates (Fig. 4a, c). The differences in net diversification rate between the two ecological categories therefore underpin their respective evolutionary trends in diversity: the relative proportion of infaunal bivalves increased gradually while epifaunal taxa correspondingly decreased (Supplementary Fig. 27), supporting an infaunalization process observed by previous authors which continued to the present. This infaunalization has been linked to an early onset of the Mesozoic marine revolution by some authors, but the lack of abundant predators in the Triassic suggests an alternative interpretation. Although some predators presumably originated in the Late Triassic to Jurassic, according to our results, very few shifts in epifaunal and infaunal origination rate occurred during this period. Some groups developed novel adaptations that may have been helpful in the arms race between predator and prey, but these innovations did not cause either group to diversify rapidly and substantially. Besides, shallow infaunal taxa, the key group of infaunal bivalves, had a higher extinction rate than epibyssate ones (mean μ = 0.025 vs mean μ = 0.016, respectively; Fig. 4e, f, i, j), indicating that an infaunal lifestyle did not protect them from extinction. Alternatively, the slightly higher extinction rate and even higher net diversification rate of infaunal bivalves may exhibit their intrinsic ecological success independent of any response to predation pressure. As a result, our results demonstrate that Triassic–Jurassic bivalves show no obvious sign of the escalation-driven diversification that left a mark on their evolutionary rates. Further studies on smaller temporal and spatial scales may illuminate when and how predators affected the evolution of bivalves (and brachiopods).

Diversification dynamics between geographic regions
Diversification rates and diversities of brachiopods and bivalves show distinct regional heterogeneity, highlighting the importance of biogeographic nuance in broader patterns of diversity (Supplementary Figs. 19–26). Except for in south-western Tethys, brachiopods and bivalves in all other regions displayed an increase in net diversification rate after the PTME. Both north-western Tethys and eastern Tethys possess the most abundant fossil occurrences and bear a higher similarity in rates to the global trends, corroborating previous findings that the signal from the ‘global’ fossil record may instead be disproportionately driven by data from one or two well-sampled regional records. In these two regions, brachiopods showed a slightly delayed increase in origination rate after the PTME compared to bivalves. Although the multivariate birth-death (MBD) model was not applied to regional data to explicitly infer the drivers of origination and extinction, the estimated rates indicate that brachiopod origination was neither suppressed by bivalve diversification, nor stimulated by their extinction. Further, bivalves did not necessarily have a higher net diversification rate than brachiopods.

The role of competition
We conducted multivariate birth-death (MBD) analyses on the diversification dynamics for the Permian–Jurassic period as a long-time window analysis, which is further subdivided into four smaller time windows: Asselian–Wordian, Capitanian–Ladinian, Carnian–Toarcian, and Aalenian–Tithonian. Both long- and short-time window analyses have advantages and disadvantages. Long-time window analysis is more useful in revealing factors that contribute to large-scale and long-term processes that shape the basic pattern of a clade’s diversification history, but if the effects of predictors are highly variable through time, the combined effects of predictors may not explain the observed diversification rates. Short window analysis reveals time-varying relationships, however, it risks over-parameterisation, and the noise of a largely stable factor may be magnified in a short-time window analysis, affecting the evaluation of other factors.

We first performed the MBD analysis on brachiopods and all bivalves, and found that many factors correlated with their diversification dynamics (Fig. 5a–d). Of these, their diversities are one of the most outstanding factors shaping their rates. If brachiopods and bivalves competed with each other, increased diversity of one clade should reduce the origination rate and/or increase the extinction rate (= reduce the net diversification rate) of the other. Bivalve diversities in some intervals (Permian–Jurassic, Capitanian–Ladinian, and Carnian–Toarcian) show strong negative correlations with origination rates of brachiopods (Fig. 5a; Supplementary Tables 1, 3, 4), suggesting competitive pressure as conventionally argued. However, such an interpretation requires caution as the negative relationship between diversity and origination rate also affected bivalves themselves (Fig. 5c). In fact, this negative correlation is largely caused by the rebounded origination after the extinction when diversity is low. For brachiopods and bivalves, this ‘diversity-dependent’ pattern is only prominent in the aftermaths of mass extinctions. During ‘normal’ times without mass extinction (e.g. Asselian–Wordian, Aalenian–Tithonian), this ‘diversity-dependence’ was not uncovered, and brachiopods were not ‘suppressed’ by bivalves (Supplementary Tables 2, 5, 7, 10). The reason why bivalve diversity rather than brachiopod diversity played a role (Fig. 5) is simply because the bivalve diversity trajectory (Fig. 3c) is closer to the shape of origination rate curves (Fig. 3a, d). Bivalve diversity rapidly decreased then increased across the PTME and TJME. As a result, the sharp declines in diversity correspond well to the spikes of extinction and origination rates (Fig. 3). Accordingly, it is likely that diversity loss during mass extinction stimulated the increase in origination rates of both brachiopods and bivalves, but brachiopods were not necessarily competitors of bivalves (i.e., the common cause effect).

Of the various bivalve life-modes, epifaunal taxa have larger overlaps with brachiopods (which are mostly epifaunal) in food and habitat than infaunal bivalves. Therefore, epifaunal bivalves should be the direct competitor of brachiopods if bivalves and brachiopods do have this relationship. Infaunal bivalves, in contrast, have more diverse feeding behaviours and occupy different niches from brachiopods. Although the activities of infaunal bivalves may disrupt substrates and thus affect occupation by brachiopods, it is a quite different mechanism from direct competition. For this reason, the MBD analysis was also conducted for the three ecological groups, epifaunal brachiopods, epifaunal and infaunal bivalves. Similar to the result described above (Fig. 5a–d), the new results also indicate that ‘diversity-dependence’ is an important mechanism that regulates their origination rates (Fig. 5e–j): the infaunal bivalve diversity strongly and negatively correlated with origination rates of all three groups. This ‘diversity-dependence’, or the elevated origination rate after mass extinction, is a ubiquitous phenomenon across multiple marine clades, which has been attributed to weakened competition, extinction of predators, or others. Regardless of specific ecological meaning, these correlations again suggest that the long-term (i.e., Permian–Jurassic) diversification dynamics of the three groups were largely governed by comparable factors, and bivalves (or specifically, epifaunal bivalves) were not competitors of brachiopods over macroevolutionary scales.

Liow et al. estimated diversification and sampling rates of post-Cambrian brachiopods and bivalves using an alternative method, and analysed the ‘causative’ or ‘correlative’ relationship of these rates based on stochastic differential equations. They detected a ‘causative’ link between bivalve extinction and brachiopod origination and concluded that brachiopods and bivalves are more than ‘ships that pass in the night’. Nevertheless, similar to our results (especially the per-capita rates estimated in discrete time bins as their method; Supplementary Figs. 39, 40, 43–45), they found that the diversification rates of brachiopods and bivalves correlated with one another. As a result, brachiopod extinction also ‘caused’ brachiopod origination and therefore, it is difficult to tell whether it was bivalves, or other factors affecting both clades (i.e., extinction events) that triggered the origination of brachiopods. To deal with this problem, Reitan and Liow re-analysed the data in multivariate models. The best model confirmed their previous view that bivalve extinction ‘caused’ brachiopod origination, but some relationships in this model, such as brachiopod origination causing bivalve origination, bivalve extinction causing brachiopod extinction, do not fit well with ecological theories (although they interpreted these relationships as reflecting similar origination and extinction rates of brachiopods and bivalves).

The absence of competition on macroevolutionary scales does not mean that brachiopod-bivalve competition cannot occur in small temporal and spatial units. Their similar ecologies would suggest that such competition could occur. For example, through cage experiments, Thayer observed that mussels can increase the mortality of brachiopods especially when predators are absent. In the fossil record, intermittently segregated distribution of brachiopods and bivalves was also observed in some regions. However, both modern and fossil assemblages suggest that such competitive exclusion rarely occurs when predators and other disturbances are present. Besides, when brachiopods and bivalves coexist, the physiological ‘superiority’ of bivalves does not guarantee their dominance. In today’s ocean, there are some regions where brachiopods are not affected by bivalves, but instead are much more abundant than bivalves. The opportunistic settlement of these animals means the colonisation of a species is largely a matter of chance; the result is many species can coexist, regardless of their competitive abilities.

Extrapolating microevolutionary dynamics to macroevolutionary time scales is potentially problematic. Nevertheless, for competition-driven macroevolutionary patterns, it is intuitively assumed that competition could decrease population size or growth rates of the inferior clade, raising its extinction probability. If the competition between brachiopods and bivalves was intense and pervasive, then should have affected their diversification rates, which is not observed in our results. Instead, while competition exclusion driven by bivalves could have happened sporadically, the cumulative evolutionary pressure of such interactions was negligible.

The role of top-down drivers
Many external factors correlated with the diversification dynamics of brachiopods and bivalves. For origination rates, the predator diversity and the ‘diversity-dependence’ discussed above had the highest absolute value of coefficients and contributed greatly to the long-term pattern. For extinction rates, sulfur and strontium isotope excursions were usually negatively correlated with them, while sea temperature and continental fragmentation generally had a positive correlation parameter (Fig. 5).

A comparison of brachiopod and bivalve responses to predictors over various time windows shows that almost no environmental predictors displayed a clear positive correlation with the bivalve diversification but a negative one with brachiopods (i.e., the selection of bivalves against brachiopods). Sea temperature is an exception. In the long-time window analysis (Permian–Jurassic), brachiopod and bivalve origination rates showed different responses to temperature (Fig. 5). When focusing on the critical time of the brachiopod-bivalve switch (i.e., Capitanian–Ladinian; Supplementary Tables 3, 8, 13, 18, 23), this difference is much more significant as temperature had a very high negative coefficient with brachiopod origination rate, while the correlation coefficient with bivalve origination rate is positive (although insignificant). This indicates that high temperatures in the Early–Middle Triassic possibly inhibited the origination rate of brachiopods but did not affect (or perhaps even stimulated) the bivalve origination rate. This interpretation is supported by the timing of recovery of the two groups. Seawater temperature increased rapidly to exceptionally high values after the PTME and in this hot environment bivalves diversified with a high origination rate (Fig. 3). But for brachiopods, even though rare Mesozoic-type taxa originated in the Induan, their main diversification phase occurred in the Olenekian to Anisian when temperatures cooled (Fig. 3). Greater bivalve tolerance to hyperthermal conditions has also been noted in the Toarcian oceanic anoxic event when brachiopods experienced a more prominent extinction and reduction in body size than bivalves. Among modern taxa, although data are sparse, brachiopods seem to be more sensitive to warming than bivalves. Together with a higher extinction rate of brachiopods at the PTME, these results support the hypothesis that physiological differences between brachiopods and bivalves caused the brachiopod-bivalve switch.

In contrast to the ‘competition’ hypothesis, Stanley suggested that the diversification of advanced predators caused the decline of brachiopods, but our analyses do not fully support this scenario. Predator diversity positively correlated with the extinction rate of Permian–Jurassic brachiopods, which seems to support Stanley’s hypothesis. However, the increasing predatory diversity was also accompanied by an increase in brachiopod origination rate. More importantly, the Cenozoic diversification of predators did not cause a substantial increase in extinction rate or decrease in origination rate for brachiopods because their net rates were even slightly higher in the Cenozoic than Cretaceous (Fig. 2e). This is possibly because brachiopods experienced intensified predation before the Cenozoic and had become marginal components of benthic ecosystems, so they were no longer subject to predators. Alternatively, predation may not be the key force driving the decline of brachiopods, which is reinforced by the observation that brachiopods are typically not the first targets of predators.

The decline of brachiopods was also attributed to their limited ability to expand their distribution after the PTME. It has been shown that in the Triassic and Jurassic, the longitudinal distribution of brachiopod genera was restricted, possibly resulting in a higher extinction risk. However, our MBD analyses suggest that continental fragmentation seems not to have had a noticeable effect on brachiopods (Fig. 5). Rather, it correlated positively with the bivalve extinction rate. This relationship is difficult to explain because it was not seen in all short-time window analyses, and further studies are needed to test it. For brachiopods, restricted geographical distribution did not always negatively affect their diversity, as narrowly distributed taxa also tend to have a higher origination rate. For instance, many Anisian brachiopods were endemic taxa that did not persist into the Late Triassic, but they contributed substantially to the Middle Triassic recovery of brachiopods.

Study limitations
Although we modified some data in the PBDB, there are always problems with taxonomy of fossils, which is inevitable in any palaeontological studies. In some cases, the hinge morphology of bivalves and internal characters of brachiopods are not well preserved, preventing precise identification. Considering the long temporal span and large size of the dataset, these taxonomic issues should not severely affect the estimation of rates. The analyses of a particular group or a region, nevertheless, might be strongly influenced, especially when the dataset is small.

A further limitation is that we were not able to include all possible environmental proxies, such as seawater oxygen content and pH values because of the lack of long-term and continuous records. The hypothesis that the brachiopod decline was caused by habitat loss was not tested because it is difficult to assess the size of habitats preferred by brachiopods. Further, for the geochemical proxies analysed, the temporal density of datum points is highly variable. ‘Hot’ intervals like the Permian-Triassic boundary are well documented, while some intervals such as the Late Triassic are poor-ly studied, hindering the discovery of environmental events and potential relationships.

The analyses of regional datasets extract the spatially heterogenous rates concealed by the global signal, but they could be greatly biased by regional sedimentation, tectonic and other factors due to the relatively small size of the spatial windows. In addition, the current curves of some geochemical proxies are composed of data from all over the world, which can be spatially heterogeneous. As a result, attempting regional MBD analyses over long timescales is challenging. Other confounding factors, such as the precise ages of fossil occurrences and uneven sampling are unavoidable, although these uncertainties and heterogeneities are partially addressed by our Bayesian methodologies. Finally, concerns have been raised about whether diversity and diversification rates can be used as proxies for assessing biotic interactions in the fossil record (e.g., competition, predation-driven evolution) in place of the detailed data on individuals and populations available to modern ecological studies. Resolving this question is beyond the scope of this paper, but our results nonetheless indicate that large-scale diversity patterns do not show any trace of competitive replacement between brachiopods and bivalves. Future studies on finer temporal and spatial scales (for both palaeontological and geochemical data) may resolve some of these problems and enable more precise tests of the hypotheses we present here.

Concluding remarks
Although brachiopods and bivalves had different diversity trends, their diversification dynamics were of broadly similar patterns before the Jurassic. Through detailed study of Permian–Jurassic taxa, we confirmed that the well-known brachiopod-bivalve switch was caused by the PTME, which killed most brachiopods. In addition, the different responses of brachiopods and bivalves to high seawater temperature may have further accelerated this switch. Contrary to conventional wisdom, we find no signs of long-term competition between brachiopods and bivalves. In the Triassic–Jurassic, the surviving brachiopod lineages successfully re-diversified without being suppressed by the bivalves. The vast difference between bivalve and brachiopod diversities today is largely due to great diversity loss of brachiopods in the PTME and rapid bivalve diversification in the Cretaceous and Cenozoic which brachiopods failed to match. Consequently, the post-extinction (PTME) ecological replacement became evident mainly from the Late Jurassic onwards, 100 Myr after the PTME. Why the bivalve origination was rapid and unconstrained and what suppressed brachiopod origination after the Jurassic when origination contributed more to the diversity may also help answer why brachiopods and bivalves are so different in diversity and abundance in the modern oceans. This study shows that brachiopods and bivalves had their own evolutionary paths. Their diversification histories were strongly controlled by rapid changes in environments (e.g., mass extinctions), and their differing physiological abilities to cope with these changes shaped their divergent fates.

Methods
Fossil occurrence data
To investigate the long-term patterns of bivalve and brachiopod diversification dynamics, we first investigated their post-Cambrian (Ordovician–Quaternary) fossil records. Cambrian records were not considered because bivalve fossil occurrences from this time are so rare. We then focused on their rates from the Permian to Jurassic, which is the critical interval when the brachiopod-bivalve switch took place. All analyses were carried out at genus level.

Fossil occurrence data of brachiopods and bivalves were manually downloaded from the Paleobiology Database (PBDB) on 23/11/2022. Occurrences with ambiguous generic assignments were excluded, but those with an unambiguous genus name and an open species name (i.e., the species name is sp.) were retained. Terrestrial records (indicated by the environment field in the PBDB) were excluded. Uncertain records (primary generic resolution noted by cf., aff., ?, informal, ex. gr., “”) with identical primary and accepted generic names were discarded, while those with different names (i.e., the species was re-assigned to another genus) were retained. In addition, absolute ages (maximum and minimum ages) of the occurrences were updated based on ages in the Geological Time Scale 2020 using the chrono_scale function of the fossilbrush R package. Next, occurrences with a high temporal uncertainty (>10 Myr) but not from an international stage were removed from the dataset. For example, a Guadalupian (~15 Myr uncertainty) occurrence was discarded, but a Norian record was kept although the latter has a ~18 Myr uncertainty.

This cleaning procedure was first applied to the post-Cambrian datasets. No taxonomic modifications were applied to the datasets, aside from deletion of bracketed subgenus names from within genus names. After data cleaning, 179,030 and 153,011 occurrences were retained for the bivalve and brachiopod datasets, respectively, including a total of 2484 bivalve genera and 3427 brachiopod genera. For the Permian–Jurassic analysis (~298.9–143.1 Ma), in addition to the above revisions, more detailed curations were made to generate a taxonomically and temporally more precise database. It should be noted that, although we focused on diversification dynamics from the Permian to the Jurassic, all Carboniferous–Cretaceous occurrences were kept during the revision stage; these were useful in analysing the ranges of genera (the pacmacro_ranges function; see below). Moreover, to ensure the rates were accurate near the boundaries of the studied interval (i.e., earliest Permian and latest Jurassic) and not affected by edge effects, the dataset analysed actually included the Kasimovian–Valanginian (~307.0–132.6 Ma) occurrences.

The additional revision process of the Permian–Jurassic dataset includes four main steps.
(1) In addition to the PBDB data, we added 1522 and 4538 occurrences of Permian–Jurassic brachiopods and bivalves, respectively. These occurrences were all reported from China and have not yet been included in the PBDB. The bivalve occurrences added are largely from the Triassic; the brachiopod occurrences added are mainly from the Permian, and a few Jurassic brachiopod occurrences are also included. These records are important especially for the regional analyses, considering that Triassic bivalve records and the Jurassic brachiopod records are relatively sparse in eastern Tethys in the PBDB.
(2) The name of every genus in the dataset was examined. Names that are currently not used and synonyms of other taxa were revised. In contrast to simply deleting the subgenus name as we did for the post-Cambrian data, every subgenus was checked, and whether it should be upgraded to genus level was decided based on recent literature. If it was regarded as a genus, the subgenus name was retained and the genus name was discarded; if not, only the genus name was retained.
(3) The duration of every genus was checked, and doubtful records were discarded. The duration of each genus was based on well-curated databases of taxon stratigraphic ranges such as Treatise on Invertebrate Palaeontology, Sepkoski’s compendium (accessed through the fossilbrush package), and recent literature. The pacmacro_ranges function in the fossilbrush package was used to generate an additional, statistical set of stratigraphic references. This function uses the stratigraphic density of occurrences to highlight doubtful occurrences. We used the default settings of the function to examine the occurrences, and those falling outside the 90% density intervals of their taxon stratigraphic ranges were treated as doubtful. These doubtful occurrences were checked, and highly unreliable ones checked against well-curated databases were discarded. It is true that deleting these data but not emending them will unavoidably omit some useful records, but revising them all is infeasible considering the size of the dataset. The unrevised very old or poorly described records would have negatively affected estimation of diversification rates. Thus, removing these doubtful records is better than retaining them in the dataset.
(4) The ages of the occurrences were updated. Where possible, occurrences were revised to substage level based on lithostratigraphic information (i.e., formation, member) and recent literature that report the ages of the stratigraphical units. For Permian records, each stage was divided equally into early and late parts. The absolute ages of the local Triassic and Jurassic stages of New Zealand were updated according to the correlation of biozones and Geological Time Scale 2020.

After these revisions, 58,319 occurrences in 684 genera, and 71,987 occurrences in 1352 genera were left in the bivalve and brachiopod datasets, respectively.

The PTME caused the extinction of many Palaeozoic brachiopod orders. However, this catastrophe had little effect on high-level clades of bivalves. To compare the diversification patterns between the Mesozoic brachiopods and bivalves, the Permian–Jurassic brachiopods were then categorised into two groups. One group (PT) consists of brachiopod orders that went extinct in the PTME (including the Spiriferida, Productida, Orthotetida, Orthida, and Dictyonellida). The other group (PT) is composed of brachiopod orders that survived the PTME (including the Lingulida, Terebratulida, Athyridida, Rhynchonellida, Spiriferinida, and Craniida) and originated after this event (the Thecideida).

Ecological categories of brachiopods and bivalves
All the Kasimovian–Valanginian brachiopods and bivalves were assigned to an ecological category according to their mode of life and relative position of the animal and sediment. Bivalve lifestyles were categorised into five groups: epibyssate (note: including semi-infaunal endobyssate here), cemented, reclining, shallow infaunal, and deep infaunal. Epibyssate bivalves usually live above the sediment-water surface (epifaunal), while endobyssate forms have parts (semi-infaunal) or all (infaunal) of their shells buried in the sediments. Although some shell forms are typical of epibyssate or endobyssate groups, transitional forms in morphology between the two are also commonly present. Semi-infaunal endobyssate taxa have parts of their shell exposed to the water and are vulnerable to predation, and so were classified together with fully epifaunal elements, but fully infaunal endobyssate genera were placed in the infaunal group. The recliners have a much lower taxonomic consistency compared with other ecological groups. These include bivalves with various types of shell morphology: taxa with a thin and flat shell (e.g., Bositra), a thick and heavy shell (e.g., Neomegalodon), or a highly inequivalve shell (e.g., Exogyra). Some free-living forms that could possibly swim were present during the Permian–Jurassic, but they were very rare and not a big ecological group. We categorised these possible ‘swimmers’ into two groups based on the following criteria: if the byssus is retained in adults, it was classified in the epibyssate group; if the byssal notch is absent and the overall shell morphology is closer to the adept swimmers, it was classified in the reclining group. The epibyssate, cemented, and reclining taxa belong to the epifaunal group. The infaunal bivalves include shallow and deep infaunal groups. Shallow infaunal taxa are usually nonsiphonate; they live close to the water column and are easily exhumed by currents. Deep infaunal taxa are often siphonate and have a greater chance of avoiding predation. Borers that live in hard substrates were classified in the deep infaunal group given that they have a stronger resilience to waves and predators. The assignment to eco-group was based on the lifestyle of the adults inferred from their shell morphology, following Stanley, Ros-Franch et al., Mondal and Harries, and other literature. Since a genus may have more than one mode of life, the major lifestyle of the genus was selected. A question mark (fewer than 1% of taxa) was given if the life strategy is poorly known or highly variable.

Unlike that of bivalves, the life strategy of brachiopods is mostly inferred because many fossil brachiopods do not have living analogues. The ecology of brachiopods given here was categorised entirely based on taxonomic classification rather than detailed morphology-function analysis. We categorised the brachiopods into four ecological types: pedicle-attached, cemented, reclining, and infaunal. The first three were clustered as the epifaunal group. Pedicle-attaching type is the most common lifestyle among the Rhynchonelliformea; we assigned all the orthides, athyridides, rhynchonellides, spiriferides, spiriferinides, and terebratulides to this ecological category. The cemented type includes those attached by the shell surface and/or with spines to support the body, such as the orthotetides, craniides, strophalosiidines, and lyttoniidines. The reclining type consists of the chonetidines and productidines. Although all have a concavoconvex shell, the strophalosiidines usually possess a cicatrix near the ventral umbo and rhizoid spines anchoring them to other objects, while the adult chonetidines and productidines normally lack these features, but instead develop long spines near the hinge line or on the body to stabilise them on the substrate. The infaunal type consists only of the deep infaunal lingulides. There must be exceptions among these ecological types. For example, the pedicles of some rhynchonellides, terebratulides, and spiriferides may lose their function when the shells grow large, and therefore, the animal turns into a recliner. These exceptions were not considered in this study.

Estimation of diversification rate
Previous studies usually used taxonomic diversity to investigate the brachiopod-bivalve switch. However, diversity arises from the interplay of origination and extinction processes and the same diversity pattern can be generated by completely different diversification dynamics. For example, an increase in diversity may have resulted from elevated origination rate or decreased extinction rate or both, and they correspond to different mechanisms. To investigate the drivers of diversification, separating origination and extinction processes is essential. Diversification dynamics of all datasets were analysed in a Bayesian framework implemented in PyRate (v.3). Incorporating both the fossil preservation process and the birth-death process, PyRate jointly estimates the preservation rates (q), the times of origination and extinction of each genus (Ts and Te), and the origination and extinction rates (λ and μ). Differentiating from previous methods that estimate rate changes at discrete time points (usually stage boundaries), PyRate considers the diversification of a given group as a continuous process. Using the reversible-jump Markov Chain Monte Carlo (rjMCMC) algorithm, PyRate can explicitly estimate the number and temporal placement of statistically significant rate changes, minimising the risk of under- and over-parameterisation. Analyses of simulated data indicate that PyRate generates reliable estimations of diversification rates under a variety of preservation models and outperforms traditional methods.

Because of our large datasets that include a great number of taxa (and parameters), it is challenging for the MCMC chain to converge. Therefore, we followed a two-step procedure to estimate diversification rates for our data. First, we used PyRate with the Gibbs MCMC sampling algorithm (in place of the default Metropolis-Hastings algorithm) to estimate the preservation-corrected Ts and Te of each taxon by sampling directly from their approximate posterior distribution. However, the rapid convergence achieved by this algorithm is at the cost of reduced resolution in estimated diversification rates. To solve this problem, the second programme, LiteRate, was employed. LiteRate was developed to estimate the origination and extinction rates of large datasets. It does not consider the preservation process at all but has the same birth-death model as PyRate. Thus, preservation-corrected times (Tss and Tes) generated by PyRate can be supplied to LiteRate to calculate high-resolution origination and extinction rates. In the end, all the preservation and birth-death parameters are properly estimated. Besides, Tss and Tes estimated by PyRate were then used to run the multivariate birth-death analyses in PyRateMBD (see below).

To incorporate the age uncertainties of the fossil occurrences, ten replicates of the original dataset were generated. Occurrences of the same fossil assemblage/site, using the collection number as the indicator, were assigned a coeval age randomly taken from their stratigraphic age range. The analysis was done for every replication, and the parameters estimated by the ten replicates were combined to plot the final result. For the preservation process, we used the -PP model test function in PyRate to select the best fitting model from among a homogeneous Poisson process, a non-homogeneous Poisson process, and a time-variable Poisson process. The result supported the time-variable Poisson processes, which allows the preservation rate to vary as a piece-wise series of constant-rate time bins. The times of preservation rate shift were set to stage boundaries in the Permian–Jurassic analyses and to series boundaries in the post-Cambrian analyses. The preservation rates across different time bins were assigned a gamma distribution; the shape parameter was set to 1.5 and the rate parameter was set to 0 to allow PyRate to estimate the rate directly from the data (-pP 1.5 0).

Unlike previous studies, we did not model potential preservation rate heterogeneity between lineages (the Gamma model, -mG option) arising from differences in environment, ecological style, shell mineralogy, and other such biases. In theory, this heterogeneity is commonly present and should be included in analyses. However, incorporation of the Gamma model increased the complexity of the overall model. When it is used along with the Gibbs model, for some datasets, the posterior values occasionally show sudden jumps even after a great number of iterations (Supplementary Fig. 28). For those replications that seemingly reached ‘convergence’ (the posterior samples are stable), the preservation rates may vary greatly between randomised replications, which indicates that the ‘convergence’ is spurious and that prohibitively many more generations would be needed by our datasets. We did some experimental analyses of our datasets using different models (Gibbs vs. Metropolis-Hastings algorithms, Gamma vs. non-Gamma models) and compared the estimated ps, Tss, and Tes. The results indicate that the Gibbs sampling can generate highly comparable ps, Tss, and Tes to the Metropolis-Hastings algorithm (Supplementary Fig. 29), especially when the Gamma model is not used (the estimates of these parameters are almost identical; Pearson’s r > 0.99990, p < 0.001). As the Metropolis-Hastings algorithm can accurately recover the parameters of the preservation and birth-death processes, the preservation-corrected Tss and Tes generated by the Gibbs-non-Gamma model and the origination and extinction rates subsequently estimated by LiteRate can also achieve this.

PyRate with the Gibbs algorithm was run for 10 million iterations. To reduce the size of log files generated, 1000 posterior samples were saved. All posterior samples were assessed using Tracer (v.1.7.1). Convergence was confirmed if the effective sample size (ESS) was >200 for all parameters, and a proper burn-in percentage was suggested. For replications with ESSs < 200, longer iterations were analysed. Mean times of origination and extinction were generated using the -ginput function in PyRate after a 10% burn-in, and then analysed by LiteRate. LiteRate was run for 200 million generations, saving 1000 posterior samples. Similarly, longer generations were run for replications that did not reach convergence. After a 10% burn-in, the posterior estimates across all replications were combined, then the mean and 95% highest posterior density intervals of origination, extinction, and net (i.e., origination minus extinction) rates were calculated (-combLog and -plotRJ functions in PyRate).

Diversity estimation
We used two methods, the -ltt (lineage through time) function of PyRate and the recently developed Bayesian model-based mcmcDivE, to estimate the diversity of brachiopods and bivalves. The former is calculated based on Tss and Tes generated by PyRate. The latter considers the uneven sampling of the fossil record and uses the preservation rates generated by PyRate and the occurrence data to estimate sampling-corrected diversity. For all datasets (including global, ecological, and regional datasets), mcmcDivE was run for 2 million generations at 1 Myr intervals. Convergence was assessed using Tracer (v.1.7.1). More generations were analysed if ESSs < 200. The first 10% of the posterior samples were discarded as burn-in. The median and 95% highest posterior density intervals of diversity were calculated.

The two methods were developed based on completely different theory. Instead of considering the unsampled data between the first and last appearance of lineages (i.e., Tss and Tes) in the -ltt function, the mcmcDivE was designed to compensate for taxa that are not sampled for a specific time period (e.g., 1 Myr). We compared the results of these methods, finding that mcmcDivE tends to produce noisier diversity trajectories (Supplementary Fig. 30). According to simulated analyses, mcmcDivE was shown to be more accurate than other methods, but whether it still performs well with real data (with high temporal and spatial sampling bias) needs further confirmation. Because the following multivariate analyses use Tss and Tes to decipher relationships between rates and factors, -ltt diversities were adopted in the MBD analyses.

Potential drivers of diversification dynamics
Except for the possible brachiopod-bivalve interactions, other potential factors/predictors may also affect the evolutions of the two clades. Here we selected ten predictors (Supplementary Fig. 31), including biotic and abiotic factors to infer their influences on the diversification dynamics of brachiopods and bivalves.

Biotic factors include three categories, brachiopod and bivalve diversities, and predator diversity. The mean diversities generated by the -ltt function were referred to as brachiopod and bivalve diversities. The intensity of predator-prey interactions is also an important factor that might affect the fates of brachiopods and bivalves. Both brachiopods and bivalves are subject to predation, and the radiations of advanced predators, such as neogastropods, crabs, and teleost fishes were hypothesised to have caused the late-Mesozoic decline of brachiopods. Previous studies showed that potential predators of the Permian–Jurassic shelled animals include echinoderms, arthropods, gastropods, fishes, and reptiles. Here diversities of two benthic groups, arthropods and echinoderms, were selected to represent the diversity of predators. Fishes and reptiles were not considered because their diversity patterns are poorly known and their fossil records are prone to be strongly affected by sampling bias such as the Lagerstätte effect. Besides, the largest group that consumes shelled animals today, the teleosts, were not diverse until the Late Cretaceous. In the class Gastropoda, although some clades are predatory, known carnivorous forms originated later than the Jurassic. The possibly parasitic platyceratid gastropods occurred in the Palaeozoic, but they are very few in the Permian–Jurassic interval. Overall, only arthropods and echinoderms have a relatively continuous and abundant fossil record. Fossil occurrences of arthropods and echinoderms were downloaded from the PBDB on 23/11/2022. Similar to the post-Cambrian brachiopod and bivalve datasets, these data were also cleaned using the same procedure. For arthropods, the Eurypterida, Xiphosura, Decapoda, and Stomatopoda were regarded as potential predators and the data of other groups were discarded; for echinoderms, the occurrences of Asteroidea, Echinoidea, and Ophiuroidea were kept. The data of the two groups were analysed separately using PyRate (the method is identical to that for brachiopods and bivalves), and diversity was calculated using the -ltt function of PyRate. Then, the means of their estimated diversities were summed as a predictor to be used in multivariate analyses (Supplementary Fig. 31). It is difficult to avoid some biases in the estimate of predation intensity due to the uncertainty of trophic structure in ancient oceans and an incomplete fossil record. The curve we generated, although imperfect, nonetheless follows the expectation that predators became increasingly diverse in the Mesozoic.

The seven abiotic factors include seawater temperature, eustatic sea level, marine carbon isotope composition (δ13C), sulfur isotope composition (δ34S), strontium isotope composition (87Sr/86Sr), global continental fragmentation, and proportion of carbonate occurrences (Supplementary Fig. 31). These factors have been widely tested to explain the diversification dynamics of marine organisms including brachiopods and bivalves. The trajectory of seawater temperature came from Scotese et al.; their ‘global average temperature’ estimated from oxygen isotopes (δ18O) was adopted. The Permian eustatic sea level was from Haq and Schutter. The Triassic and Jurassic sea level was from Haq. The carbon isotope data were from Cramer and Jarvis, and these were smoothed using the locally-weighted polynomial regression. The dataset of δ34S was from Present et al. Isotopic values derived from bulk rock carbonate‐associated sulfate were discarded because they show greater variability than isotopic values from other sources and may not reflect the δ34S composition of palaeo-seawater sulfate. Secular change curves of 87Sr/86Sr followed McArthur et al. The original global continental fragmentation index was based on the Seton et al.’s tectonic plate model. Here, we recalculated the index value based on the PALEOMAP tectonic model, to be consistent with that used in the regional analysis. A high value means all plates are not touching, and a low value suggests the presence of a supercontinent. The last factor, the proportion of carbonate occurrences, was selected to describe the habitat (substrate) condition. The vast majority of the post-Permian brachiopods were pedicle-attached. Mesozoic brachiopods became more and more frequent in carbonate lithologies which are on average firmer than siliciclastics, and the loss of habitats because of strengthened bioturbation possibly drove the Mesozoic–Cenozoic decline of brachiopods. To approximate habitat availability, we calculated the proportion of occurrences (of all brachiopods and bivalves) reported from carbonate lithologies at the substage level. This factor also indicates whether the temporal distributions of occurrences were affected by alternations of rock lithology, and these rate variations, therefore, might be attributed to the quality of fossil record. The proportional number of carbonate formations is not employed as the proxy because the number of occurrences varies greatly among formations. It is inappropriate to treat a highly fossiliferous formation and a fossil-poor unit as an equally weighted unit in the analyses, although the proportions of the formations correlate positively with occurrences of fossils (Spearman’s ρ = 0.39, p < 0.01).

Multivariate birth-death analyses
The multivariate birth-death (MBD) model was employed to assess the influence of biotic and abiotic factors on the diversification dynamics of brachiopods and bivalves in the software package PyRateMBD. In this model, origination and extinction rates are correlated with time-continuous variables (i.e., factors, predictors) through an exponential function or linear function, so that the observed rates represent the function of the baseline rates, correlation parameters, and external factors. Using an MCMC algorithm, PyRateMBD can jointly estimate the baseline origination (λ0) and extinction (μ0) rates, and correlation parameters (Gλ and Gμ). A positive G implies that the factor correlates positively with the rate and vice versa. Besides, by using a horseshoe prior, the analysis can assess the support of the assumed relationship, preventing over-parameterisation. A correlation parameter with a high shrinkage weight (w) close to 0 is likely to be noise, while a value close to 1 represents a true signal.

We performed the multivariate analysis for five datasets. The first two are the whole Permian–Jurassic brachiopods and bivalves, and their diversities are included as predictors. This analysis is similar to previous studies in treating bivalves as a whole group. The other three datasets analysed are epifaunal brachiopods (which shows very little difference from the whole brachiopod dataset because infaunal brachiopods are rather low in diversity), epifaunal bivalves, and infaunal bivalves, so the three diversity trajectories were included as predictors. Because of greater overlap in food resource and habitat, epifaunal bivalves were hypothesised as competitors of brachiopods. However, no studies separately compared epifaunal and infaunal bivalve diversities with brachiopod diversities (although Sepkoski noticed this problem). By dividing bivalves into two ecological groups, it will help to distinguish their potentially different relationship with brachiopods.

Before analysis, all factors were rescaled to 0 and 1 to remove biases from the absolute magnitudes of input values. By default, PyRate uses stepwise interpolation to connect the discrete values of the predictors. To precisely track the variation of the values, factors were instead linearly interpolated at 0.1 Myr intervals. We ran the MBD model for 30 million generations and saved 1000 posterior samples. For every dataset, ten replications were analysed separately (the Tss and Tes estimated by PyRate were used as input, avoiding re-modelling the complex preservation process), and the ten log files generated were combined after a 10% burn in. Convergence was assessed using Tracer (v.1.7.1), and longer generations were analysed for replications with ESSs < 200. The 95% highest posterior density intervals of correlation parameters and mean shrinkage weights were calculated. If the shrinkage weight of a correlation is greater than 0.5, the relationship was regarded as significant.

As suggested by other authors, the drivers of diversification may change over time, and the relationships revealed over an extended period may fail in being significant within a short period. Therefore, we first conducted the MBD analyses for the whole Permian–Jurassic dataset. Next, the entire Permian–Jurassic interval was subdivided into four shorter windows: Asselian–Wordian, Capitanian–Ladinian, Carnian–Toarcian, and Aalenian–Tithonian, and each was analysed using the MBD model. These time windows have different durations (from ~27.3 to ~62.3 Myr), but all are composed of six or seven stages/ages. Given that the age of occurrence in the datasets relies largely on stages/substages, subdivision of time bins based on the number of stages rather than the absolute range (e.g., 40 Myr) is more appropriate. For all MBD analyses, both exponential and linear models were run. We calculated Bayes factors to investigate the relative support of the two models following Lehtonen et al. Except for the Asselian–Wordian datasets which did not show clear support for either one model (absolute difference of the log Bayes factors is <2), the exponential model was supported in all cases. Therefore, only the results of the exponential model were reported.

Regional analyses
The global diversity and diversification dynamics represent the summed diversity across a set of geographically and environmentally distinct regions, and the real diversification process occurring in a region may be masked by the global trend. To investigate whether diversification patterns of brachiopods and bivalves are different among geographical regions, we spatially standardised the data and analysed their respective rates. The workflow of spatial standardisation followed Flannery-Sutherland et al. Instead of dividing the regions arbitrarily, we made spatial windows according to the palaeobiogeographical pattern of fossil occurrences. First, the occurrences were spatially binned using a hexagonal grid through the icosa R package. Second, a grid-taxon matrix was made based on the presence/absence of taxa in these grids. Last, we used two methods, network analysis and partitioning around medoids clustering, to identify biogeographic regions. Network analysis was implemented using the igraph R package. Partitioning around medoids clustering was implemented using the fpc R package, and the modified Forbes index was used to calculate the dissimilarity matrix for partitioning around medoids clustering. Then the result of grouping can be used as a reference to construct the spatial windows. This procedure was implemented in an R function, compare_biogeography, and can be directly applied to study the palaeobiogeography of other fossil datasets. Based on the analyses of brachiopod and bivalve occurrences (Supplementary Figs. 32–35), four biogeographical regions having relatively abundant and continuous fossil records (especially the Permian–Triassic interval) were selected to analyse the regional rates, including northern Panthalassa, north-western Tethys, south-western Tethys, and eastern Tethys (Supplementary Fig. 36).

Prior to spatial standardisation, the palaeo-coordinates of the fossil occurrences were updated based on the modern coordinates and revised midpoint ages, using the PALEOMAP tectonic model through the GPlates Web Service (https://gwsdoc.gplates.org). A spatial window was constructed for each of the four regions using the spacetimewind function. The generated spatial window is relatively stable in size, and moves over time to track the movement of the focal plate to guarantee that the analytical results were not biased by the species-area effect and continental drift. Fossil occurrence data were subsampled into each region and then standardised using the spacetimestand function. We applied the minimum spanning tree (MST) length as the metric to standardise the occurrence data. In this procedure, fossil occurrences are binned into hexagonal grids, and an MST is reconstructed from the grid centres that contain fossil data. If the MST length is greater than the target value, the cells with the smallest amount of data are removed from the MST, until the threshold is reached. We calculated the MST length of raw data in each time bin, and the median MST length of all time bins was adopted as the threshold. After spatial standardisation, the data show much less variance in MST length (Supplementary Fig. 37). For rate analysis, the PyRate + LiteRate procedure is inapplicable because many regions have some bins that lack any records of fossil occurrences, and then LiteRate cannot estimate the rates successfully. As a result, only PyRate with the Gibbs algorithm was employed for the analyses. The rates estimated clearly show various patterns among regions despite the loss of some resolution.

Diversification dynamics revealed by traditional methods
In addition to Bayesian model-based methods which estimate continuous rates, we also analysed the diversification dynamics of global brachiopods and bivalves using traditional methods (e.g., per-capita, three-timer, and its improved versions, gap-filler and second-for-third), which are based on discrete time bins. The workflow was identical to that of Kocsis et al. Every stage was regarded as a time bin, except the Neogene (four bins) and Quaternary (one bin). Occurrences were assigned to these time bins based on their ages, and those that could not be safely assigned to a certain bin were discarded. Rates and diversities were calculated using the divDyn R package. Raw diversity and sampling-corrected diversity (using the Shareholder Quorum Subsampling algorithm) were estimated by the two metrics, range-through and sampled-in-bin diversities (corrected by the three-timer sampling completeness). Foote’s per-capita rates and Alroy’s second-for-third rates are reported (Supplementary Figs. 39–42). Compared with the per-capita rates of brachiopods and bivalves, the second-for-third rates are more volatile, especially for the post-Jurassic brachiopods. After checking the output of PyRate and divDyn, we attribute the high volatility of post-Jurassic brachiopod rates to the low preservation rate (sampling rate) during this time. Simulation studies have indicated that when diversity or sampling is low, three-timer and its related methods tend to generate results that are more volatile than the per-capita rate. Generally, the per-capita rate is more accurate than other traditional methods, although it has its own drawbacks. Therefore, we did additional correlation tests on per-capita rates of brachiopods and bivalves. Rates in the first and last two bins were not included to remove the influence of edge-effects.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review, Guo et al., MS# NCOMMS-23-12673-T\nSummary: This is an ambitious paper, trying to tackle the recalcitrant problem of understanding the\ncontrasting Phanerozoic diversity trajectories of bivalves and brachiopods. The paper sets the stage\nbroadly, but while it makes clear and concise reference to the long (and diffuse) past literature on\ntrying to understand the contrasting trajectories, it only incidentally (even if often) touches on these\nvarious previous ideas and data. Instead, it focuses on the authors’ own assessment of the diversity\ndynamics (origination and extinction rates, and the resulting diversity trajectories) of the two\ngroups, their correlations with each other and with a variety of environmental variables, touching on\ndifferences geographically, as well as breaking the analysis up ecologically. This a very ambitious\nresearch agenda. But in the end, this manuscript is too long, tackles too many issues, most of which\nare treated too diffusely, fails (for the most part) to make rich contact with the past literature, as\nwell has having some methodological issues. They cover so much ground that I don’t know what the\nnon-specialist reader is supposed to take away from the paper that is new, and I think the\ndiscussion points are too simple for the specialist reader. Taken together, I feel these all combine to make the\npaper unsuitable for Nature Communications. Below I provide a little more detail on just some of\nthese issues.\n1) Abstract. The abstract reflects the diffuseness of the manuscript, viz (their text in italics): (a)\nTimes of major biotic replacement have traditionally been interpreted in a Neodarwinian\nsense as mediated by innovation in the succeeding clades. I am not sure how useful this\nstatement is, and is incomplete, for the traditional approach (is there any other?) also makes\nreference to changes in the environment (which they do point out in their text). (b) … the\nswitch from brachiopods to bivalves … [has been] attributed to competitive exclusion of\nbrachiopods by the superiorly adapted bivalves or simply to the fact that brachiopods had\nbeen hit especially hard by the PTME. OK, fine, but the idea of long-term competitive\nexclusion I do not think has any support. Then we change gears abruptly with (c) … the\ntwo\nclades displayed highly comparable trends of diversification before the Jurassic (which I can’t\nsee in the data, and even if present, I don’t understand the significance). Then back we go to\nthe Permian-Triassic transition with (d) Insight from a multivariate birth-death model shows\nthat the extinction of major brachiopod clades during the PTME set the stage for the BBS (I\nthink we have known this for 50 years), with differential responses to high ocean\ntemperatures post-extinction further facilitating their displacement by bivalves (ah, an\ninteresting result, ignoring the fact that correlation is not causation). Then, a final\nconclusion, that largely confirms past conclusions (although contradicted by some of the\nliterature which they don’t resolve): (e) Our study indicates [I think ‘confirms’ would be a\nmore accurate statement] that brachiopods and bivalves were not competitors over\nmacroevolutionary timescales, with extinction events and environmental stresses shaping\ntheir divergent fates. This is so broad as to be not so informative. So, I don’t see what non-\nspecialists (or even specialists) are supposed to take away from this paper in its current\nform, although there are many aspects of analysis that will be of interest to specialists if\npresented more thoroughly in the specialist literature.\n2) Similar trends between pre-Jurassic bivalve and brachiopod diversity dynamics. There are\nseveral issues here. First, this conclusion seems to be via visual inspection only, and I cannot\nsee it in the data, which makes we concerned that it is not true (unless the authors are\nsimply referring to shared volatility compared with the post-Triassic records). This needs\nsome quantification. Related to this, are they referring to just the largest changes, i.e., the\nfact that both groups respond to the mass extinctions? I would also note that even if the\ntrajectories are in fact similar, this does not mean that the clades are not in competition (as\nimplied in the abstract): there could have been biotic interaction between the two clades on\nshorter timescales, i.e., persistent competition on ecological timescales, with both\nresponding similarly to larger perturbations (just as competitor companies might also\nrespond similarly to global changes in consumer spending, thus showing similar large scale\ntrajectories even though in competition). So, an analysis of the residuals might be of value\nhere.\n3) Methodological concerns:\na) The brachiopod diversity trajectory (see Figure 1). This does not match my sense of their\ntrajectory, largely due to the huge Permian peak. I suspect this is due to the authors\naddition of the very rich (largely Permian) Chinese data, which makes me suspect the\ncurve is biased by heterogeneous temporal and spatial sampling.\nb) This then raises the general issue of how to deal with sampling biases. It appears that they\nhave let PyrRate deal with preservation issues. It is crucial that the authors explore the\nrobustness of their conclusions by also treating their data with some sort of sample\nstandardization process, for example SQS, and then analyzing them with alternative\napproaches, for example divDyn, to compute origination and extinction rates, etc. That is,\nI simply don’t know if I can trust the analysis – if more than one analytic approach leads to the\nsame conclusions, that would greatly increase the power of the paper.", "reviewer_2": "Overall, the submitted manuscript is strong. The results are entirely predictable: the hypothesis\nthat bivalves outcompeted brachiopods is theoretically reasonable, but has been clearly and\nunambiguously falsified since good global biodiversity data became available starting in the late\n1970s. Nevertheless, there are practitioners who continue to promote (without evidence) this\noutdated hypothesis. In this context, the manuscript is important. Even though, the results are\nexpected, presenting new lines of evidence that bivalves did not outcompete brachiopods after the\nPermian mass extinction (and vice versa before the extinction) is a valuable contribution to the\nliterature.\nThe methods are strong and the writing is clear. I have a few minor suggestions, but I don't have\nany major reservations or concerns with the manuscript.\nThe most serious comment I have is that I think the language regarding the patterns of bivalve\nand brachiopod diversification before and after the Jurassic is greatly overstated. The first\nparagraph of the \"Diversification dynamics of ...\" Section on page 5 states that both clades show\nsimilar dynamics before the Jurassic, but the trends become distinct after the Jurassic. I don't\nthink this claim is supported by Figure. 2. My read on figure 2 is that the the Jurassic is, indeed,\nan inflection point in the macroevolutionary history of the two clades. However, the dynamics flip\naround that point. Before the Jurassic, the bivalves show relatively stable rates (Fig 2a and 2b)\nwhile the Brachiopod rates are volatile (Fig 2d and 2e). After the Jurassic, they flip and the\nbivalves are volatile and the brachiopods are stable. I don't think this negates the overall\ninterpretation of the data, but it does require a bit more explanation. I wonder if it has to do with\nthe overall diversity? Does the birth-death model employed inherently produce more volatile rates\nwhen diversity is greater than 200-300 (genuine question)? Alternatively, perhaps asking reader to\nvisually inspect raw time series is not the best way to evaluate the relative similarity of the trends.\nMight a cross plot of bivalve vs. brachiopod rates would be more effective??\nMinor Comments & Suggestions:\nI don't like the first part of the title \"Red in tooth and gill\". This is a variation on a common\nphrase, but is culturally specific and won't make sense to all readers. Admittedly, though I am\nfamiliar with the original phrase \"red in tooth and claw\", I had to look up the exact meaning.\nEven if this stems primarily from ignorance on my part, I suspect it will be unfamiliar to a\nlarge portion of Nature Communication's readership who do not speak English as their native\nlanguage.\nOn of the things that plagues paleontology is the overuse of acronyms. I think the manuscript uses\na few too many acronyms. Specifically, I feel the use of BBS for brachiopod-bivalve switch is\ngratuitous: I recommend dropping the acronym and just writing out brachiopod-bivalve switch (or\ntransition). Moreover, the methods have turned into a bit of an acronym bonanza. The use of some\nis fine, but I would, again, recommend spelling out some of the acronyms that are only used a\ncouple of times (e.g., ESS, HDPI, PAM). There are just a lot of acronyms to remember and it's a\ndistraction while reading to have to stop and either think about what they mean or go back in the\ntext to find their first use.\nThere are a number of places in the manuscript where two items are listed and then the second\nitem is referred to as \"latter\" in a subsequent clause. As with the acronyms, this can be a\ndistraction to the reader, especially if \"latter\" does not occur directly after the list. I recommend\nwriting out the object being referred to where possible. A good example is on Page 8 Lines 190-\n192. \"...Cretaceous and Cenozoic, which the latter ...\". I was initially confused because I thought\nlatter referred to the Cenozoic. I had to reread the sentence to remind myself that it was\nbrachiopods being referenced.\nPage 5 Line 54: I think \"largest events in the history of marine life\" should should reference the\nend-Permian mass extinction, not the \"BBS\".\nAll time series figures: Please drop the vertical gray bars marking alternating time intervals. These\nare chart junk and make the plots difficult to read without conveying any important information.\nImportant events are directly labeled (e.g., PTME), which render the gray bars superfluous. In\naddition, the yellow lines and labels for the PTME, etc. would be easier to read if the were a bit\ndarker.\nPage 7 Line 170: I think a word other than \"observations\" should be used here. I don't think the\nobservations made in this analysis strengthen or weaker previous observations. Previous\nobservations have already been made and are forever fixed. It's the interpretations of those\nobservations that are strengthened, weakened, supported, etc.\nPage 9 Line 218: I don't think \"greater\" is the best word to use here. I think \"more\" or \"a greater\nnumber\" would be better. Greater implies better or larger in size.\nPage 15 First Paragraph. I think a very short discussion, or at least a citation, of the very\nimportant paper by Adrain and Westrop (2000, Science: An empirical assessment of toxic\npaleobiology) is appropriate here.\nPage 17 Line 429. I humbly request that the authors add the > 6000 fossil occurrences mentioned\nhere to the PBDB. This is, of course, not a requirement for publication and should not be\nconsidered in the editorial decision, but it would be a great service to the community. If the\nauthors contact the database executive committee they can help arrange a bulk upload from\nexcel/csv files so that the data does not need to be keystroked into the database by hand.\nPage 21 Line 544: \"PyRate considers the diversification of clades\". What are the clades here? Are\nthey simply Bivalvia vs. Brachiopoda or are rates being calculated for orders, families, etc. I think\nclarification on this point would be helpful.\nPage 22 Line 587: Should \"modal\" be \"model\"?\nPage 22 Line 591: \"Fake\" is not the right word to use here. Fake implies something that was\ncreated intentionally as a deception. I think spurious is a better word choice: it has a standard\nstatistical usage that makes sense in this context.\nPage 25 Line 663: I as a little confused why Echinoidea and Ophiuroidea were included as\npredators of bivalves and brachiopods. While it's likely that both may consume shelled\nbenthic invertebrates (particularly larval/juvenile stages), I don't think they constitute a major\nsource of food for either group: I am not convinced they should be included. If there is\nevidence/literature contrary to my understanding, please include a justification for their inclusion.\n-Noel Heim"}, "author_response": "RESPONSE TO REVIEWER COMMENTS\nReviewer #1 (Remarks to the Author):\nSummary: This is an ambitious paper, trying to tackle the recalcitrant problem of\nunderstanding the contrasting Phanerozoic diversity trajectories of bivalves and brachiopods.\nThe paper sets the stage broadly, but while it makes clear and concise reference to the long\n(and diffuse) past literature on trying to understand the contrasting trajectories, it only\nincidentally (even if often) touches on these various previous ideas and data. Instead, it\nfocuses on the authors’ own assessment of the diversity dynamics (origination and extinction\nrates, and the resulting diversity trajectories) of the two groups, their correlations with each\nother and with a variety of environmental variables, touching on differences geographically,\nas well as breaking the analysis up ecologically. This a very ambitious research agenda. But\nin the end, this manuscript is too long, tackles too many issues, most of which are treated too\ndiffusely, fails (for the most part) to make rich contact with the past literature, as well has\nhaving some methodological issues. They cover so much ground that I don’t know what the\nnon-specialist reader is supposed to take away from the paper that is new, and I think the\ndiscussion points are too simple for the specialist reader. Taken together, I feel these all\ncombine to make the paper unsuitable for Nature Communications.\nREPLY: Many thanks for your valuable comments on our study. We note that you are\nconcerned about: 1) lack of focus (diffuseness), 2) too many points, and 3) omission of some\nimportant references to the earlier literature. We have considered all of these points carefully\nand revised the manuscript thoroughly.\nWe added some discussions about our new results and previous hypotheses, especially\nthese three important BBS scenarios: Sepkoski (1996; competition hypothesis), Steele-\n         \nseries analyses showing the competition between the two). Additional important references\nhave been cited in the revised version.\nWe re-analysed the same dataset using traditional methods (see responses below). The\nper-capita rates largely agree with the rates derived from the PyRate analysis. Moreover,\ncorrelation tests of these rates support our conclusion that the pre-Jurassic rates of\nbrachiopods and bivalves have similar trends.\nWe also performed additional MBD analyses on the three groups: epifaunal brachiopods,\nepifaunal bivalves, and infaunal bivalves. Because of larger overlap in habitat and feeding\nmechanism, epifaunal bivalves were hypothesized to be competitors of brachiopods. The\nindependent multivariate analysis of epifaunal and infaunal bivalves is an interesting issue\nthat was overlooked by previous studies. These analyses also connect the ‘Rates of different\necological groups’ section and the multivariate section, making the manuscript more coherent.\nBelow I provide a little more detail on just some of these issues.\n1) Abstract. The abstract reflects the diffuseness of the manuscript, viz (their text in italics):\n(a) Times of major biotic replacement have traditionally been interpreted in a Neodarwinian\nsense as mediated by innovation in the succeeding clades. I am not sure how useful this\nstatement is, and is incomplete, for the traditional approach (is there any other?) also makes\nreference to changes in the environment (which they do point out in their text).\nREPLY: Thanks, we have revised the statements.\n(b) ... the switch from brachiopods to bivalves ... [has been] attributed to competitive\nexclusion of brachiopods by the superiorly adapted bivalves or simply to the fact that\nbrachiopods had been hit especially hard by the PTME. OK, fine, but the idea of long-term\ncompetitive exclusion I do not think has any support.\nREPLY: Although the brachiopod-bivalve switch is a classic example of biotic replacement,\nonly a few studies have investigated their long-term diversification dynamics in detail. Gould\nand Calloway (1980) denied the presence of long-term competition between the two groups.\nIn contrast, Sepkoski (1996), Liow et al. (2015), and Reitan and Liow (2017) re-studied this\nissue later on, and these authors found some evidence, and clarified the idea of competitive\nreplacement between brachiopods and bivalves. Whether or not the evidence is reliable, the\nanswer to this issue still remains open to debate. We agree with you that the idea of long-term\ncompetitive exclusion has been rarely supported, but obviously, new studies based on the\nupdated datasets using advanced analytical methods are needed.\nThen we change gears abruptly with (c) ... the two clades displayed highly comparable trends\nof diversification before the Jurassic (which I can’t see in the data, and even if present, I don’t\nunderstand the significance).\nREPLY: Your suggestions are followed in the revised version (See responses below).\nThen back we go to the Permian-Triassic transition with (d) Insight from a multivariate birth-\ndeath model shows that the extinction of major brachiopod clades during the PTME set the\nstage for the BBS (I think we have known this for 50 years), with differential responses to\nhigh ocean temperatures post-extinction further facilitating their displacement by bivalves\n(ah, an interesting result, ignoring the fact that correlation is not causation).\nREPLY: Thanks, we agree with you that correlation is not causation. We have re-worded\nsome sentences in the manuscript.\nThen, a final conclusion, that largely confirms past conclusions (although contradicted by\nsome of the literature which they don’t resolve): (e) Our study indicates [I think ‘confirms’\nwould be a more accurate statement] that brachiopods and bivalves were not competitors over\nmacroevolutionary timescales, with extinction events and environmental stresses shaping\ntheir divergent fates. This is so broad as to be not so informative. So, I don’t see what non-\nspecialists (or even specialists) are supposed to take away from this paper in its current form,\nalthough there are many aspects of analysis that will be of interest to specialists if presented\nmore thoroughly in the specialist literature.\nREPLY: It is true that our results ‘confirm’ the previous hypothesis rather than proposing a\nnew scenario on the brachiopod-bivalve replacement. However, the previous studies (Gould\nand Calloway, 1980; Sepkoski, 1996) usually utilize diversity data, rather than diversification\ndynamics to figure out this issue. The same diversity trajectory can be achieved by\ncompletely different combinations of origination and extinction rates. Thus, disentangling the\ntwo components is important to answer this question. Liow et al. (2015) did so, but they\nconcluded that both brachiopods and bivalves are competitors, denying the conclusion given\nby Gould and Calloway (1980). In addition to competition, other factors have also been\nproposed to explain brachiopod and bivalve extinction or origination, especially across the\nPTME. However, most of these hypotheses were not tested in a multivariate analysis. In\nrecent years, more accurate and advanced methods such as PyRate have been developed. We\nemploy these new methods to re-study the long-debated brachiopod-bivalve switch issue.\nTherefore, although our analyses did not lead to startling or shocking new conclusion,\nwe argue for the merits of our work in two ways: (a) the brachiopod vs. bivalve replacement\nis a textbook example, widely known to many evolutionary biologists as well as\npalaeobiologists, and so our work will be important to those general readers; and (b) we reach\na decisive conclusion based on use of the most detailed data set ever deployed and the widest\narray of state-of-the-art techniques, resolving a long-term but unresolved debate point. So, we\nhave a headline that will startle many analysts and commentators in evolutionary study:\nbrachiopods were not driven to extinction by ‘superior’ bivalves.\n2) Similar trends between pre-Jurassic bivalve and brachiopod diversity dynamics. There are\nseveral issues here. First, this conclusion seems to be via visual inspection only, and I cannot\nsee it in the data, which makes we concerned that it is not true (unless the authors are simply\nreferring to shared volatility compared with the post-Triassic records). This needs some\nquantification. Related to this, are they referring to just the largest changes, i.e., the fact that\nboth groups respond to the mass extinctions? I would also note that even if the\ntrajectories are in fact similar, this does not mean that the clades are not in competition (as\nimplied in the abstract): there could have been biotic interaction between the two clades on\nshorter timescales, i.e., persistent competition on ecological timescales, with both\nresponding similarly to larger perturbations (just as competitor companies might also\nrespond similarly to global changes in consumer spending, thus showing similar large scale\ntrajectories even though in competition). So, an analysis of the residuals might be of value\nhere.\nREPLY: Thank you for your comments and suggestions. It is true that the “similar trend” that\nwe stated before was based on visual inspection. Unlike traditional methods that estimate\nrates for every time bin, PyRate with the rjMCMC algorithm estimates the number and\ntemporal placement of statistically significant rate changes, which reduces the risk of over-\nparameterisation. Therefore, the rates can be “flat” if there is no significant rate shift, and the\nabundant tied data renders direct and effective correlation analysis.\nFollowing your suggestions below, we calculated diversification rates using traditional\nmethods, and we also performed additional analyses to the per-capita rates. Correlation tests\nof these rates agree with our previous conclusion that origination and extinction rates of\nbrachiopods and bivalves correlate positively with one another, and this relationship is more\nprominently demonstrated by the pre-Jurassic data.\nWe agree with you that there could be biotic interactions between the two clades in some\ncases, as Thayer’s experiments showed. However, whether these sporadic, localized events\nwere the main cause of the macroevolutionary decline of brachiopods is questionable. Our\nstudy focuses on large-scale diversification dynamics.\nYour suggestions on the analysis of residuals are also highly appreciated. Following\nyour advice, we performed the reduced major axis regression analysis on their extinction and\norigination rates, respectively, and then calculated the residuals (Fig. S49). Clearly, the\nresiduals do not correlate with brachiopod or bivalve diversities (p > 0.7), i.e., an\nunnoticeable role of competition on macroevolutionary timescales.\n3) Methodological concerns:\na) The brachiopod diversity trajectory (see Figure 1). This does not match my sense of their\ntrajectory, largely due to the huge Permian peak. I suspect this is due to the authors addition\nof the very rich (largely Permian) Chinese data, which makes me suspect the curve is biased\nby heterogeneous temporal and spatial sampling.\nREPLY: Many thanks for your comments. Overall, the diversity trajectory of brachiopods\nshown in Figures 1 and 2 is close to the trajectories that were previously published elsewhere,\nincluding those based on:\n(1) PBDB database and “sample in bin” diversity (e.g., Alroy, 2010, sampling\nstandardised using the SQS algorithm; Close et al., 2020, spatially and SQS standardised);\n(2) PBDB database and “range through” diversity (e.g., Payne et al., 2014);\n(3) well-curated databases and “range through” diversity (e.g., Sepkoski’s Compendium\n[Sepkoski, 1996]; Treatise on Invertebrate Paleontology [Curry and Brunton, 2007]).\nAll these trajectories show high peaks of diversity in the Ordovician–Devonian and\nPermian (see the figure below), although the peaks have various heights.\nDiversity peaks of Permian brachiopods published in previous literature (Sepkoski, 1996;\nCurry and Brunton, 2007; Alroy, 2010; Payne et al., 2014; Close et al., 2020). Diversity\ncurves of Alroy (2010) and Close et al. (2020) are sampling standardised. Colours in the\ndiagram of Close et al. (2020) represent diversity in different paleolatitude zones.\nThe fossil occurrences that we added indeed do not increase the steepness of the\ndiversity peak of Permian brachiopods, because the curves in Figures 1 and 2 were totally\nbased on the raw data downloaded from the PBDB. These curves are used to reflect the long-\nterm diversification dynamics of brachiopods and bivalves. For the subsequent Permian-\nJurassic analysis, we did many revisions including adding additional occurrences. Despite all\nthis, we found that only two of the Permian genera that we added are absent in PBDB.\nWe compared the taxonomic compositions in the raw PBDB dataset, the revised PBDB\ndataset, and those in the Treatise, and found that the high Permian peak of brachiopod\ndiversity in Figures 1 and 2 is mainly caused by two types of taxa. First, the raw PBDB\ndataset contains 89 Permian genera that should be discarded because these taxa are invalid, or\nsynonymous to other taxa, or should not occur in the Permian (i.e., identifications of the\nfossils are problematic). Second, many new Permian genera have been established or re-\nvalidated in recent years after the Treatise of Brachiopoda was published in 2000s. We\ncounted the number of genera that are recorded in the PBDB but were not included in the\nTreatise of Brachiopoda, and found that 114 Permian genera were established or re-validated\nafter the Treatise was published. In contrast, there are 31 and 22 new genera in Devonian and\nCarboniferous, respectively that were not included in the Treatise. As a result, given that so\nmany new genera were established from the Permian fossil records, the Permian peak of\ndiversity should be higher than the diversity peak of the same age shown by Curry and\nBrunton (2007; see figure above). It is noteworthy that when the diversity is calculated for\nSQS subsampled data (Fig. S39), the Permian diversity peak is still prominent.\nIn summary, the pronounced peak of Permian diversity is a character of the PBDB\ndataset itself, and does not appear to be a bias that arose from the various analytical methods\nthat we employed. When the raw dataset is taxonomically emended, the Permian peak is less\nsteep (the highest Permian diversity is <400 in Figure 3). This highlights the importance of\ntaxonomic emendation for the datasets that are used for big data analysis, as we did for the\nPermian-Jurassic dataset. It should be noted that although the raw post-Cambrian diversity is\nless accurate, it does not influence our subsequent analyses. The multivariate analysis is\nbased on the Permian-Jurassic diversities which were calculated from the taxonomically\nrevised datasets.\nRefs cited in this Reply:\nAlroy J. The shifting balance of diversity among major marine animal groups. Science 329,\n1191–1194 (2010).\nCarlson SJ. The evolution of Brachiopoda. Annual Review of Earth and Planetary Sciences\n44, 409–438 (2016).\nClose RA, Benson RBJ, Saupe EE, Clapham ME, Butler RJ. The spatial structure of\nPhanerozoic marine animal diversity. Science368, 420–424 (2020).\nCurry GB, Brunton CHC. Stratigraphic distribution of brachiopods. In: Treatise on\nInvertebrate Paleontology, Part H (Revised), Brachiopoda, Vol. 6. (ed Selden PA).\nGeological Society of America and University of Kansas (2007).\nPayne JL, Heim NA, Knope ML, McClain CR. Metabolic dominance of bivalves predates\nbrachiopod diversity decline by more than 150 million years. Proceedings of the Royal\nSociety B: Biological Sciences281, 20133122 (2014).\nSepkoski Jr, JJ. Competition in macroevolution: the double wedge revisited. In: Evolutionary\nPaleobiology (eds Jablonski D, Erwin DH, Lipps JH). University of Chicago Press\n(1996).\nb) This then raises the general issue of how to deal with sampling biases. It appears that they\nhave let PyrRate deal with preservation issues. It is crucial that the authors explore the\nrobustness of their conclusions by also treating their data with some sort of sample\nstandardization process, for example SQS, and then analyzing them with alternative\napproaches, for example divDyn, to compute origination and extinction rates, etc. That is, I\nsimply don’t know if I can trust the analysis – if more than one analytic approach leads to the\nsame conclusions, that would greatly increase the power of the paper.\nREPLY: It is true that incompleteness of fossil records and sampling biases are often an\nimportant problem and tough issue for fossil record-based statistical analysis and modelling.\nIndeed, the great performance of PyRate has been tested repeatedly (although some papers\nquestioned the accuracy of PyRate; see response and discussion in Flannery-Sutherland et al.\n[2022]). Analyses of simulated datasets indicated that PyRate can unravel the underlying\ndiversification processes (Silvestro et al., 2014) under a variety of preservation scenarios.\nEven if a large fraction of taxa is not sampled, PyRate is also able to reveal some true signals,\nand the confidence intervals cover the real trajectory.\nPyRate has many advantages over traditional methods (Silvestro et al., 2014, 2019), and\nits Bayesian framework allows for testing the significance of rate shifts. Silvestro et al. (2019)\nperformed many simulation tests, and showed that the relative error of PyRate results is much\nlower than those produced by traditional methods (e.g., Alroy’s three-timer, Foote’s\nboundary-crossing or per-capita). Traditional methods based on discrete time slices usually\ngenerate more volatile curves when the true evolutionary rates are constant (Silvestro et al.,\n2019), reflecting that traditional methods may have the problem of over-parameterisation.\nAlthough these uncertainties can be represented by confidence intervals from multiple\nanalyses, their widths are usually greater than those of PyRate, especially when diversity is\nlow.\nRecently, Warnock et al. (2020) proposed a fossilised birth-death (FBD) model to\nestimate the diversification of fossil data. Their analysis of simulated data suggested that the\nFBD model is more accurate than the birth-death (BD) model that PyRate uses. However, in\nmost cases, both BD and FBD models generate comparable results, which are better than that\nderived from traditional methods (especially the three-timer) (Warnock et al., 2020). The\nFBD model has been presented in PyRate, but, due to its complexity, this model is difficult to\nbe applied to our large datasets. Therefore, we believe that PyRate with the BD model is\ncurrently the most suitable method for the analysis of our datasets.\nFollowing your suggestion, we have re-analysed the diversification rates using the\ntraditional methods (those provided in the divDyn package) and diversity. It is expected that\ndifferent methods will generate slightly different results. Here we report the per-capita and\nsecond-for-third (the improved version of three-timer and gap-filler) rates. An obvious\ndifference between the two is, the second-for-third rate of brachiopods is highly volatile after\nthe Jurassic, but the per-capita rate in the same period is low, like our PyRate result. We\nchecked the divDyn and PyRate output files and attribute this difference to the low\npreservation rate of brachiopods at that time (the low sampling rate of Cenozoic brachiopods\nwas also reported by Liow et al. (2015)). Analyses of simulated data have indicated that when\ndiversity or sampling rate (preservation rate) is low, three-timer or its related methods tend to\ngenerate results that are more volatile and less accurate than the per-capita method (Alroy,\n2014; Silvestro et al., 2019; Warnock et al., 2020). Although the per-capita rate (and every\nother method) has its own drawbacks, we used it to make additional comparisons and tests\n(see Supplementary Figures). Generally, the per-capita rates estimated here show largely\ncomparable trends with our PyRate trajectories.\nRefs cited in this Reply:\nAlroy, J. Accurate and precise estimates of origination and extinction rates. Paleobiology 40,\n374–397 (2014).\nFlannery-Sutherland JT, Silvestro D, Benton MJ. Global diversity dynamics in the fossil\nrecord are regionally heterogeneous. Nature Communications 13, 2751 (2022).\nSilvestro D, Schnitzler J, Liow LH, Antonelli A, Salamin N. Bayesian estimation of\nspeciation and extinction from incomplete fossil occurrence data. Systematic Biology 63,\n349–367 (2014).\nSilvestro D, Salamin N, Antonelli A, Meyer X. Improved estimation of macroevolutionary\nrates from fossil data using a Bayesian framework. Paleobiology 45, 546–570 (2019).\nWarnock RCM, Heath TA, Stadler T. Assessing the impact of incomplete species sampling on\nestimates of speciation and extinction rates. Paleobiology46, 137–157 (2020).\nReviewer #2 (Remarks to the Author):\nOverall, the submitted manuscript is strong. The results are entirely predictable: the\nhypothesis that bivalves outcompeted brachiopods is theoretically reasonable, but has been\nclearly and unambiguously falsified since good global biodiversity data became available\nstarting in the late 1970s. Nevertheless, there are practitioners who continue to promote\n(without evidence) this outdated hypothesis. In this context, the manuscript is important.\nEven though, the results are expected, presenting new lines of evidence that bivalves did not\noutcompete brachiopods after the Permian mass extinction (and vice versa before the\nextinction) is a valuable contribution to the literature.\nThe methods are strong and the writing is clear. I have a few minor suggestions, but I don't\nhave any major reservations or concerns with the manuscript.\nREPLY: Many thanks for the positive comments.\nThe most serious comment I have is that I think the language regarding the patterns of\nbivalve and brachiopod diversification before and after the Jurassic is greatly overstated. The\nfirst paragraph of the \"Diversification dynamics of ...\" Section on page 5 states that both\nclades show similar dynamics before the Jurassic, but the trends become distinct after the\nJurassic. I don't think this claim is supported by Figure. 2. My read on figure 2 is that the the\nJurassic is, indeed, an inflection point in the macroevolutionary history of the two clades.\nHowever, the dynamics flip around that point. Before the Jurassic, the bivalves show\nrelatively stable rates (Fig 2a and 2b) while the Brachiopod rates are volatile (Fig 2d and 2e).\nAfter the Jurassic, they flip and the bivalves are volatile and the brachiopods are stable.\nI don't think this negates the overall interpretation of the data, but it does require a bit more\nexplanation. I wonder if it has to do with the overall diversity? Does the birth-death model\nemployed inherently produce more volatile rates when diversity is greater than 200-300\n(genuine question)? Alternatively, perhaps asking reader to visually inspect raw time series is\nnot the best way to evaluate the relative similarity of the trends. Might a cross plot of bivalve\nvs. brachiopod rates would be more effective??\nREPLY: Thanks for your comments. We are sorry that we did not do correlation test of these\nrates. Unlike traditional methods that estimate rates for every time bin, PyRate with the\nrjMCMC algorithm estimates the number and temporal placement of statistically significant\nrate changes, which reduces the risk of over-parameterisation. Therefore, the rates can be\n“flat” if there is no significant rate shift, and the abundant tied data renders direct and\neffective correlation analysis.\nSilvestro et al. (2014, 2019) indicated that PyRate can extract the underlying\ndiversification dynamics even if the preservation rate or diversity is low. When the\npreservation rate is low, although it cannot accurately estimate the absolute value of rates, it\ncan still detect the time of rate shifts. In such a case, the final mean rates may be flatter, but\nthe wide confidence interval still covers the real rates. Accordingly, the flat brachiopod rates\nwith narrow confidence intervals in the Cenozoic should be a true signal. The low diversity is\nalso not a problem. Many previous studies using PyRate focused on clades with low\ndiversities and rate shifts were successfully recovered. Some fossil groups that we studied\n(e.g., cemented bivalves and deep infaunal bivalves shown in Fig. 4) also have very low\ndiversity, but the mean rates have many fluctuations.\nThis is also the reason why we said the trends of brachiopod and bivalve diversification\nrates are “similar” to one another. The bivalve rate is much less volatile than the brachiopod\nrate in the early–middle Paleozoic, but at the same time, the confidence interval is wide. The\nlow volatility of the bivalve rate in this interval is possibly caused by the low bivalve\npreservation rate, but the long-term general trend is similar to brachiopods: both clades had\nrates gradually decreasing in Ordovician–Devonian.\nWe also calculated their diversification rates using traditional methods following\nReviewer 1’s suggestion (see responses above), and the result allows for a correlation test.\nIt is shown that before the Jurassic, brachiopod and bivalve rates are significantly correlated,\nbut the Jurassic–Quaternary rates are less correlated. This provides additional evidence for\nour previous conclusions.\nRefs cited in this Reply:\nSilvestro D, Schnitzler J, Liow LH, Antonelli A, Salamin N. Bayesian estimation of\nspeciation and extinction from incomplete fossil occurrence data. Systematic Biology 63,\n349–367 (2014).\nSilvestro D, Salamin N, Antonelli A, Meyer X. Improved estimation of macroevolutionary\nrates from fossil data using a Bayesian framework. Paleobiology 45, 546–570 (2019).\nMinor Comments & Suggestions:\nI don't like the first part of the title \"Red in tooth and gill\". This is a variation on a common\nphrase, but is culturally specific and won't make sense to all readers. Admittedly, though I am\nfamiliar with the original phrase \"red in tooth and claw\", I had to look up the exact meaning.\nEven if this stems primarily from ignorance on my part, I suspect it will be unfamiliar to a\nlarge portion of Nature Communication's readership who do not speak English as their native\nlanguage.\nREPLY: Thanks, we have removed this part from the title.\nOne of the things that plagues paleontology is the overuse of acronyms. I think the\nmanuscript uses a few too many acronyms. Specifically, I feel the use of BBS for brachiopod-\nbivalve switch is gratuitous: I recommend dropping the acronym and just writing out\nbrachiopod-bivalve switch (or transition). Moreover, the methods have turned into a bit of an\nacronym bonanza. The use of some is fine, but I would, again, recommend spelling out some\nof the acronyms that are only used a couple of times (e.g., ESS, HDPI, PAM). There are just a\nlot of acronyms to remember and it's a distraction while reading to have to stop and either\nthink about what they mean or go back in the text to find their first use.\nREPLY: Many thanks. We have deleted the unnecessary acronyms and spelled out some of\nthe acronyms.\nThere are a number of places in the manuscript where two items are listed and then the\nsecond item is referred to as \"latter\" in a subsequent clause. As with the acronyms, this can be\na distraction to the reader, especially if \"latter\" does not occur directly after the list. I\nrecommend writing out the object being referred to where possible. A good example is on\nPage 8 Lines 190-192. \"...Cretaceous and Cenozoic, which the latter ...\". I was initially\nconfused because I thought latter referred to the Cenozoic. I had to reread the sentence to\nremind myself that it was brachiopods being referenced.\nREPLY: Thanks, we have revised them.\nPage 5 Line 54: I think \"largest events in the history of marine life\" should should reference\nthe end-Permian mass extinction, not the \"BBS\".\nREPLY: Thanks for the comments. We have revised it.\nAll time series figures: Please drop the vertical gray bars marking alternating time intervals.\nThese are chart junk and make the plots difficult to read without conveying any important\ninformation. Important events are directly labeled (e.g., PTME), which render the gray bars\nsuperfluous. In addition, the yellow lines and labels for the PTME, etc. would be easier to\nread if the were a bit darker.\nREPLY: Thanks very much for your suggestions. We have removed all the gray bars in the\nmain figures to make them clearer. But we have kept them in the Supplementary Figures to\nmake it easy for the reader to compare these time series.\nPage 7 Line 170: I think a word other than \"observations\" should be used here. I don't think\nthe observations made in this analysis strengthen or weaker previous observations. Previous\nobservations have already been made and are forever fixed. It's the interpretations of those\nobservations that are strengthened, weakened, supported, etc.\nREPLY: Thanks for your comments. We have deleted this statement.\nPage 9 Line 218: I don't think \"greater\" is the best word to use here. I think \"more\" or \"a\ngreater number\" would be better. Greater implies better or larger in size.\nREPLY: Thanks very much. We have corrected it.\nPage 15 First Paragraph. I think a very short discussion, or at least a citation, of the very\nimportant paper by Adrain and Westrop (2000, Science: An empirical assessment of toxic\npaleobiology) is appropriate here.\nREPLY: Thanks for your suggestions. We have cited this paper.\nPage 17 Line 429. I humbly request that the authors add the > 6000 fossil occurrences\nmentioned here to the PBDB. This is, of course, not a requirement for publication and should\nnot be considered in the editorial decision, but it would be a great service to the community.\nIf the authors contact the database executive committee they can help arrange a bulk upload\nfrom excel/csv files so that the data does not need to be keystroked into the database by hand.\nREPLY: Thanks for your suggestions. The additional occurrence data that we added have\nbeen provided in Supplementary Material. They have identical column names with PBDB\ndata and can be easily read and manipulated in R. We will try to upload these occurrences to\nPBDB so it will be good for future studies.\nPage 21 Line 544: \"PyRate considers the diversification of clades\". What are the clades here?\nAre they simply Bivalvia vs. Brachiopoda or are rates being calculated for orders, families,\netc. I think clarification on this point would be helpful.\nREPLY: Apologies for this confusion. It should be “a group”. The diversification of any\ngroups is regarded as a continuous process by PyRate.\nPage 22 Line 587: Should \"modal\" be \"model\"?\nREPLY: Thanks. Corrected.\nPage 22 Line 591: \"Fake\" is not the right word to use here. Fake implies something that was\ncreated intentionally as a deception. I think spurious is a better word choice: it has a standard\nstatistical usage that makes sense in this context.\nREPLY: Thanks. Corrected.\nPage 25 Line 663: I as a little confused why Echinoidea and Ophiuroidea were included as\npredators of bivalves and brachiopods. While it's likely that both may consume shelled\nbenthic invertebrates (particularly larval/juvenile stages), I don't think they constitute a major\nsource of food for either group: I am not convinced they should be included. If there is\nevidence/literature contrary to my understanding, please include a justification for their\ninclusion.\nREPLY: Thanks for your comments. We agree that bivalves and brachiopods both are\nprobably not major foods of the Echinoidea and Ophiuroidea. But it is clear that Echinoidea\ncan consume bivalves when they lived together (Penchaszadeh et al., 2004). Bivalves or\nbrachiopods may not be a big component of their food resource, but Echinoidea can also\nincrease brachiopod mortality by grazing (Tomašových, 2008), which is also an important\nmechanism of the Mesozoic marine revolution (Vermeij, 1977). For the Ophiuroidea, it seems\nthat both brachiopods and bivalves are not major sources of food. However, their predation\non brachiopods is repeatedly mentioned in papers and the Treatise (Foukeand LaBarbera,\n1986), and therefore, we included them in the analyses. The Ophiuroidea is only a small part\nof the predators. Its diversity doesn’t show rapid changes. The rise of the predator diversity\ntrajectory in the Middle–Late Jurassic is largely contributed by crustaceans and echinoids. In\nfact, it is very difficult to estimate the predation intensity in ancient ecosystems, and the\nproportions of contributions from different groups are difficult to estimate, so we only\nroughly estimate the predation pressure by adding diversities of multiple predatory groups.\nWe admit that this curve is far from perfect, but we think that it should be more or less useful\nin a long-window analysis.\nFouke, B.W., and LaBarbera, M. Ecology of shallow water brachiopods in Jamaica: Proof of\npredation and its implications. 4th North American Paleontological Convention\nAbstracts 4, A16 (1986).\nPenchaszadeh, P. E., Bigatti, G. & Miloslavich, P. Feeding of Pseudechinus magellanicus\n(Philippi, 1857) (Echinoidea: Temnopleuridae) in the SW Atlantic coast (Argentina).\nOphelia 58, 91–99 (2004).\nTomašových, A. Substrate exploitation and resistance to biotic disturbance in the brachiopod\nTerebratalia transversa and the bivalve Pododesmus macrochisma. Marine Ecology\nProgress Series 363, 157–170 (2008).\nVermeij, G. J. The Mesozoic marine revolution: Evidence from snails, predators and grazers.\nPaleobiology 3, 245–258 (1977).\n-Noel Heim"}, {"round": 2, "reviewer_comments": {"reviewer_1": "REVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nSummary: This is the second time I have seen this paper. I felt the first version was too diffuse,\ndid not make proper contact with the literature, and coupled with my feeling that this issue was\nalready settled, to warrant publication in this venue. However, up reading this revision, I am\nconvinced that this is indeed should be published in Nature Communications – it is a very good\npaper, with several sophisticated yet clearly expressed arguments relevant to establishing\ncausality of diversity change in the fossil record, not just for bivalves versus brachiopods, but in\ngeneral. This is an important paper, well executed.\nBelow are a series of largely minor comments that I think will help make the paper more readable\n(frankly, I don't know how anyone can write a sophisticated paper in another language, but these\nauthors have done a fine job). In the interests of full disclosure, I did not scrutinize the methods\nsection.\nSuggestions for improvement\nLines 21-22: The authors’ say that “we find that unexpectedly the two clades displayed highly\ncomparable trends of diversification before the Jurassic”. It seems to me that the details are\nactually quite different (I don't see the ‘highly comparable trends’), but they DO share similar\nsame large-scale trends, so perhaps this generality (‘similar large-scale trends’) is the phrase that\nbelongs in the abstract (not ‘highly comparable trends’).\nLine 38: The reference to the ‘Big Five’ (Raup and Sepkoski 1982) is over 40 years old; maybe the\nauthors’ could cite the 2023 review of the Big Five in the new journal “Cambridge Prisms:\nExtinction”.\nLine 79: ‘dependence’ mis-spelled.\nLine 92: Replace ‘any difference arising between’ with ‘differences between’.\nLine 102: I still can’t see the ‘comparable trends of diversification dynamics, especially\nbefore the Jurassic’ in the Figure (my eye is going to the fine detail of the wiggles). But if the next\nsentence explains what you mean, then I agree. So maybe just insert ‘viz.’ between the sentences to\nindicate that you are going to tell is in what way they are similar – then I am OK with the\nstatement.\nLine 116: Insert ‘their’ before ‘subsequent’.\nLine 117: I don't really understand this sentence. Are you simply trying to say that before the\nJurassic (do you mean the Tr/J boundary?) there are broad similarities, and that afterwards the\ntrajectories are largely divergent?\nLine 118: I think the word ‘contradictory’ needs to replace with something else (very different?).\nLine 119: I suspect there is a better word/phrase to replace ‘resolutions’.\nLine 124. Replace ‘perturbated’ with ‘volatile’?\nLine 130: Replace ‘reinforced by’ with ‘found in’?\nLine 132: Replace ‘Besides’ with ‘Interestingly”, or just delete it without replacement?\nLine 138: Here is this fuzzy use of the Jurassic, which just seems a bit arbitrary, and a little odd\nsince itself in an interval, not a point in time.\nLine 138. Replace “Owing to’ with ‘Given the’.\nLine 154: I don't see what the relationship is between the brachiopod diversity drop at the G-L\nboundary and their drop at the PTME – yes, both are big drops, but the PTME drop may well have\nbeen just ad big even if there was no G-L drop.\nLine 175. Period missing.\nLines 187-202. This analysis is particularly nice.\nLines 203-229: This analysis is also particularly nice – it lovely to see the authors use their data to\ntease apart (and refute) past claims – now I feel I am really learning something.\nLine 221: Insert “They argue that ….” At the beginning of the new sentence on this line.\nLine 233: I don’t understand what you mean by ‘predictable’. Do you mean simply that they make\nsense? Or do you just want to back off a little and begin this section by simply stating that ‘For\nbrachiopods there is a heterogeneity of rates between the different lifestyles.”?\nLine 246: Delete ‘part’.\nLine 247: Begin with ‘of the’.\nLines 267-292: I am learning something again here, nice.\nLine 294: Insert ‘geographic’ in here, so it is crystal clear that this about geography?\nLine 312: What does MBD stand for?\nLine 313: What four smaller time windows? I don’t recall reading about these yet.\nLine 317: Replace ‘cannot’ with ‘may not’?\nLine 318: Then maybe delete ‘sufficiently’?\nLines 330-354. This is a little awkward (but it is OK). FYI, I came up with the same conclusions\nlooking at your data (Fig.5) – nice analysis, and nice to see you have not fallen into the trap of a\nsimplistic explanation.\nLines 355-362. Nice reasoning: this is what I wanted to know now.\nLine 363. Delete ‘Being’?\nLines 375-391. Well done; I find Liow’s arguments ‘awkward’, and you have made explicit the\n‘awkwardness’ very well, even though it is not quite resolved.\nLine 424: I don’t know what you mean by ‘rare’.\nLine 4242: The word ‘caused’ feels like the wrong word. Correlated?\nLine 424: Here you mention the long and short time window analysis, but I don’t feel I have been\ntold what this means. Oh, do you simple mean Phanerozoic (minus the Cambrian) versus Permian-\nto-Jurassic? Please make explicit.\nLine 433: Do you mean ‘(or maybe even stimulated)’?\nLine 444-445. Interesting.\nLine 457: Replace ‘could’ with ‘may’?\nLine 466: Replace ‘was not revealed’ with ‘not seen’ (whatever short-time window analyses\nmeans)?\nLine 510: Replace ‘contradictory’ with ‘different’?\nLine 511: It seems to me that a more accurate statement would be that before the Jurassic they\nshow broadly similar diversification dynamics.\nLine 520: And due to the loss of the PTe brachiopod group?", "reviewer_2": "Reviewer #2 (Remarks to the Author):\nThe authors have satisfactorily addressed all of the comments from my previous review. I have no\nfurther comments or suggestions."}, "author_response": "RESPONSE TO REVIEWERS’ COMMENTS\nReviewer #1 (Remarks to the Author):\nSummary: This is the second time I have seen this paper. I felt the first version was too\ndiffuse, did not make proper contact with the literature, and coupled with my feeling\nthat this issue was already settled, to warrant publication in this venue. However, up\nreading this revision, I am convinced that this is indeed should be published in Nature\nCommunications – it is a very good paper, with several sophisticated yet clearly\nexpressed arguments relevant to establishing causality of diversity change in the fossil\nrecord, not just for bivalves versus brachiopods, but in general. This is an important\npaper, well executed.\nBelow are a series of largely minor comments that I think will help make the paper\nmore readable (frankly, I don't know how anyone can write a sophisticated paper in\nanother language, but these authors have done a fine job). In the interests of full\ndisclosure, I did not scrutinize the methods section.\nREPLY. Many thanks for your valuable comments and corrections which make our\nmanuscript more readable.\nSuggestions for improvement\nLines 21-22: The authors’ say that “we find that unexpectedly the two clades displayed\nhighly comparable trends of diversification before the Jurassic”. It seems to me that the\ndetails are actually quite different (I don't see the ‘highly comparable trends’), but they\nDO share similar same large-scale trends, so perhaps this generality (‘similar large-\nscale trends’) is the phrase that belongs in the abstract (not ‘highly comparable trends’).\nREPLY. Thanks for your suggestion. We have revised it.\nLine 38: The reference to the ‘Big Five’ (Raup and Sepkoski 1982) is over 40 years old;\nmaybe the authors’ could cite the 2023 review of the Big Five in the new journal\n“Cambridge Prisms: Extinction”.\nREPLY. Thanks. We have added this reference.\nLine 79: ‘dependence’ mis-spelled.\nREPLY. Corrected.\nLine 92: Replace ‘any difference arising between’ with ‘differences between’.\nREPLY. Corrected.\nLine 102: I still can’t see the ‘comparable trends of diversification dynamics, especially\nbefore the Jurassic’ in the Figure (my eye is going to the fine detail of the wiggles). But\nif the next sentence explains what you mean, then I agree. So maybe just insert ‘viz.’\nbetween the sentences to indicate that you are going to tell is in what way they are\nsimilar – then I am OK with the statement.\nREPLY. Thank you for your suggestions. We have modified it.\nLine 116: Insert ‘their’ before ‘subsequent’.\nREPLY. Corrected.\nLine 117: I don't really understand this sentence. Are you simply trying to say that\nbefore the Jurassic (do you mean the Tr/J boundary?) there are broad similarities, and\nthat afterwards the trajectories are largely divergent?\nREPLY. Yes, we are simply saying that the large-scale trends are similar from the\nOrdovician to Triassic (before the Tr/J boundary), but the trajectories are largely\ndivergent from the Jurassic to recent (after the Tr/J boundary). We have added some\nsentences to make it clear.\nLine 118: I think the word ‘contradictory’ needs to replace with something else (very\ndifferent?).\nREPLY. Thanks. Corrected.\nLine 119: I suspect there is a better word/phrase to replace ‘resolutions’.\nREPLY. We have deleted this word. The word ‘volatility’ has already expressed our\nmeaning.\nLine 124. Replace ‘perturbated’ with ‘volatile’?\nREPLY. Corrected.\nLine 130: Replace ‘reinforced by’ with ‘found in’?\nREPLY Corrected.\nLine 132: Replace ‘Besides’ with ‘Interestingly”, or just delete it without replacement?\nREPLY. Corrected.\nLine 138: Here is this fuzzy use of the Jurassic, which just seems a bit arbitrary, and a\nlittle odd since itself in an interval, not a point in time.\nREPLY. Thanks for your suggestion. We have addressed the Tr/J boundary in the\nrevised text and re-worded this sentence.\nLine 138. Replace “Owing to’ with ‘Given the’.\nREPLY. Corrected.\nLine 154: I don't see what the relationship is between the brachiopod diversity drop at\nthe G-L boundary and their drop at the PTME – yes, both are big drops, but the PTME\ndrop may well have been just ad big even if there was no G-L drop.\nREPLY. Sorry for the misleading. The reason why brachiopod diversity declined more\nsharply is because the brachiopod extinction rate is much higher than that of bivalve in\nthe PTME. We have deleted the statement about the G-L extinction.\nLine 175. Period missing.\nREPLY. Corrected.\nLines 187-202. This analysis is particularly nice.\nREPLY. Thank you!\nLines 203-229: This analysis is also particularly nice – it lovely to see the authors use\ntheir data to tease apart (and refute) past claims – now I feel I am really learning\nsomething.\nREPLY. Thank you for your encouragement.\nLine 221: Insert “They argue that ….” At the beginning of the new sentence on this line.\nREPLY. Revised.\nLine 233: I don’t understand what you mean by ‘predictable’. Do you mean simply that\nthey make sense? Or do you just want to back off a little and begin this section by\nsimply stating that ‘For brachiopods there is a heterogeneity of rates between the\ndifferent lifestyles.”?\nREPLY. The previous statement means that the heterogeneity can be easily predicted\nfrom their diversity changes. For example, the cemented and reclining brachiopods\nshould have a higher extinction rate than the pedicle-attached ones in the PTME\nbecause many orders of them became extinct in the PTME. To make the statement clear\nand readable, we have deleted the word ‘predictable’ and simply stating the presence of\nrate heterogeneity.\nLine 246: Delete ‘part’.\nREPLY. Corrected.\nLine 247: Begin with ‘of the’.\nREPLY. Corrected.\nLines 267-292: I am learning something again here, nice.\nREPLY. Thank you.\nLine 294: Insert ‘geographic’ in here, so it is crystal clear that this about geography?\nREPLY. Revised.\nLine 312: What does MBD stand for?\nREPLY. We mentioned this word in the previous ‘regional analysis’ section and\ntherefore the abbreviation was used. Since this section is all about the result of MBD\nanalysis, we have re-written the complete words (multivariate birth-death) at the\nbeginning of this section to make it clear.\nLine 313: What four smaller time windows? I don’t recall reading about these yet.\nREPLY. The long-time window in our study refers to the entire Permian–Jurassic period.\nWe then subdivided this long period into four shorter time windows (Asselian–Wordian,\nCapitanian–Ladinian, Carnian–Toarcian, Aalenian–Tithonian) to investigate if the\ncorrelations between rates and factors are changeable through time. Both long- and\nshort-time window analyses are useful, as explained in the manuscript. The four short\ntime windows were described in the Method section. We have added their names in the\nmain text to make the discussion more readable.\nLine 317: Replace ‘cannot’ with ‘may not’?\nREPLY. Corrected.\nLine 318: Then maybe delete ‘sufficiently’?\nREPLY. Corrected.\nLines 330-354. This is a little awkward (but it is OK). FYI, I came up with the same\nconclusions looking at your data (Fig.5) – nice analysis, and nice to see you have not\nfallen into the trap of a simplistic explanation.\nREPLY. Yes, it is bit difficult to explain the result, and we have tried to make our\nexplanations clear and more understandable.\nLines 355-362. Nice reasoning: this is what I wanted to know now.\nREPLY. Thanks.\nLine 363. Delete ‘Being’?\nREPLY. Corrected.\nLines 375-391. Well done; I find Liow’s arguments ‘awkward’, and you have made\nexplicit the ‘awkwardness’ very well, even though it is not quite resolved.\nREPLY. Thanks. The Reitan and Liow’s multivariate analysis is also very complex, and\nmultiple relationships were analysed at the same time. More studies are needed to\nexplain their results.\nLine 424: I don’t know what you mean by ‘rare’.\nREPLY. ‘Rare’ means ‘almost no’. Most of the factors don’t show the selectivity of\nbivalves against brachiopods, but temperature is an exceptional factor.\nLine 4242: The word ‘caused’ feels like the wrong word. Correlated?\nREPLY. Corrected.\nLine 424: Here you mention the long and short time window analysis, but I don’t feel I\nhave been told what this means. Oh, do you simple mean Phanerozoic (minus the\nCambrian) versus Permian-to-Jurassic? Please make explicit.\nREPLY. Following your suggestion, we have mentioned them at the beginning of the\n‘The role of competition’ section.\nLine 433: Do you mean ‘(or maybe even stimulated)’?\nREPLY. Yes. Revised.\nLine 444-445. Interesting.\nREPLY: Thanks.\nLine 457: Replace ‘could’ with ‘may’?\nREPLY. Revised.\nLine 466: Replace ‘was not revealed’ with ‘not seen’ (whatever short-time window\nanalyses means)?\nREPLY. Revised. The short-time window has been explained.\nLine 510: Replace ‘contradictory’ with ‘different’?\nREPLY. Revised.\nLine 511: It seems to me that a more accurate statement would be that before the\nJurassic they show broadly similar diversification dynamics.\nREPLY. Thanks. Corrected.\nLine 520: And due to the loss of the PTe brachiopod group?\nREPLY. Revised.\nReviewer #2 (Remarks to the Author):\nThe authors have satisfactorily addressed all of the comments from my previous review.\nI have no further comments or suggestions.\nREPLY. Thank you very much for reviewing our manuscript again and positive\ncomments."}]

Decision letter