The correlation of subsampled diversity with within-region geographic sample spread could be explained in one of two possible ways, as discussed above for “global” diversity. The first is that variation in paleogeographic spread among regions and intervals systematically biases our subsampled diversity estimates, artifactually enhancing the trend towards increasing subsampled diversity through time, and increasing the scatter of subsampled diversities. The second is that the correlation between paleogeographic spread and subsampled regional diversity results from a common-cause process, in which a third variable such as continental flooding drives both genuine biodiversity and changes in regional sample spreads, as is well-documented in shallow marine biotas (e.g., [ 55 ]). We cannot presently distinguish between these two hypotheses, so we discuss results for the relationship between time and subsampled diversity estimates under both models.

(A) Regional subsampled species diversity versus paleogeographic spread. (B) Regional subsampled species diversity versus paleolatitudinal centroid (positivised value). (C) Regional paleogeographic spreads versus geological age (Ma). (D) Regional paleolatitudinal centroids versus geological age (Ma). Paleogeographic spreads are the minimum spanning tree lengths in km. Subsampled values were obtained using a quorum of 0.4. Correlation coefficients and p-values from Pearson’s correlation tests are reported in the top-left panels A and B. Abbreviated interval names are given in full in S1 Table . The data displayed in this figure can be accessed at http://doi.org/10.5061/dryad.9fr76 [ 48 ].

No individual region is represented sufficiently well to provide a continuous time series of subsampled diversity. However, variation in paleogeographic sample spread and paleolatitude among regions and through time could bias a pooled regional analysis. To investigate this, we compared regional subsampled diversities to regional minimum spanning tree lengths and paleolatitudinal centroids (median paleolatitudes of collections) across the Mesozoic, using general linear models with a Gaussian error distribution and ln() link function ( Fig 4 ; Table 1 ). In univariate comparisons, subsampled diversity has significant positive relationships with geographic spread, but not with absolute paleolatitude, and a significant negative relationship with geological time ( Table 4 ), indicating higher subsampled diversity estimates in younger time bins. Comparisons of the AICc-weights [ 56 ] of regression models including combinations of time, geographic spread and paleolatitude as explanatory variables indicate that a univariate model explaining regional subsampled diversities using only geographic spread (AICc-weight = 0.50) is more likely than one using geological time (= 0.12), and substantially more likely than one using paleolatitude (= 0.02). The second best model (AICc-weight = 0.18), is one in which time is included as an explanatory variable together with geographic spread in a multiple regression. In this model, the slope of geological time is reduced to approximately half of its value in a univariate model ( Table 4 : from -0.00335 to -0.00148) and is non-significant, whereas the slope of geographic spread is reduced by less, and becomes marginally non-significant ( Table 4 : from 0.0000489 to 0.0000356). The slope of paleolatitude is non-significant in all models ( Table 4 ).

The results described above suggest that both directly counted and subsampled diversity measures are at least partly explained by the number of distinct global regions sampled. This indicates that sampling bias determines global richness patterns in the tetrapod fossil record. In fact, these “global” patterns result from a heterogeneous assemblage of regional patterns, and the presence of sample spread bias seriously undermines the use of global fossil occurrence data as an adequate summary of standing diversity through time. Therefore, we report subsampling results for contiguous continental regions defined in S2 Table ( S1 Appendix ), and not for the entire planet taken as a whole.

All three of our geographic distribution measures are strongly and significantly correlated with face-value counts of both genera and species over the full range of thresholds examined ( Table 3 ). Counts of long branches retain strong and significant relationships with taxon counts when conditioned on counts of short branches, but they have a non-significant relationship with taxon counts when conditioned on the summed lengths of short branches at thresholds of 500 and 1,000 km. Counts of long branches have strong and significant or near significant (p = 0.054; subsampled species | threshold = 1,000 km) correlations with subsampled species and genus diversity estimates ( Table 3 ), and retain substantial power to explain subsampled diversity estimates (R 2 > 0.45 [species]; R 2 > 0.35 [genera]) when conditioned on both of our of short branch measures. The variance explained by partial correlation for long branch counts is greater than that for short branch measures when conditioned on counts of long branches in almost all cases, with exceptions only at a threshold of 1,000 km ( Table 3 ).

We tested between these alternatives by examining the correlations and partial correlations between taxonomic richness measures and three measures of the geographic distribution of localities: (1) counts of geographically “long” branches of our minimum spanning trees, representing the number of distinct global regions sampled; (2) counts of geographically “short” branches of our minimum spanning trees, representing the addition of fossil localities to regions that have already been sampled; and (3) the summed lengths of geographically “short” branches of our minimum spanning trees, representing the geographic spread of localities within local regions. There is no clear distinction between “short” and “long” branches in our minimum spanning trees. However, histograms of the frequency distributions of branch lengths in the intervals with the greatest total geographic spread suggest that a frequency transition occurs somewhere between 100 and 1,000 km ( S1 Fig ), so we performed analyses at several “threshold” branch lengths between these values, and up to 2,000 km. In general, counts of long branches are strongly and significantly correlated with both of our “short branch” measures ( Table 3 ), indicating that intervals with more sampled regions also have greater total local sampling.

The minimum spanning trees for global fossil localities circumscribe planetary spatial scales from 9,500 km to 68,800 km, with a range that is comparable to the circumference of the Earth (~40,075 km). At this scale, the correlation between paleogeographic sample spread and diversity could result from either or both of two possibilities: (1) A direct bias model, in which the correlation is due to variation in the number of distinct global regions for which data are available in each time bin (assuming that each region existed even during intervals in which it is unsampled). (2) A “common cause” model, in which processes such as sea level changes, orogeny, and rifting determine the absolute sizes of individual regions, thereby determining both paleogeographic sample spread and species richness via species-area effects (e.g., [ 55 ]). Both models assume that a species-area relationship exists (i.e., that available land area constrains species diversity). They differ in that the direct bias model assumes that differences in the area sampled for each interval result from bias, whereas the common cause model assumes that they result from actual changes in the ancient Earth system.

Using (A) genus counts, (B) subsampled genera, (C) species counts, and (D) subsampled species. Paleogeographic spreads are the minimum spanning tree lengths in km. Subsampled values were obtained using a quorum of 0.4. Correlation coefficients and p-values from Pearson’s correlation tests are reported in the top-left of each panel. Abbreviated interval names are given in full in S1 Table . The data displayed in this figure can be accessed at http://doi.org/10.5061/dryad.9fr76 [ 48 ].

Uneven fossil record sampling may substantially bias directly counted diversity patterns. To address this, we applied equal-coverage or shareholder quorum subsampling (SQS) to standardise sampling among time bins. Most subsampling approaches, including SQS, require that geographic spread is held approximately constant to allow meaningful comparison of regional gamma diversities (e.g., [ 52 ]). However, both directly counted and subsampled patterns of “global” richness are biased by differences in paleogeographic sample spread, estimated as the length of the minimum spanning tree uniting the paleocoordinates of non-marine tetrapod-bearing localities for each interval ( Fig 3 ). This bias is evidenced by strong, statistically significant, positive correlations, a pattern that is similar to geographic bias in the marine invertebrate [ 53 ] and Miocene North American mammal [ 54 ] fossil records and which explains 72% of the variance in global subsampled species counts ( Fig 3D ) and 62% of the variance in directly observed species counts ( Fig 3C ).

Pooled regional face-value genus and species counts ( Fig 2 ) also show no evidence for a Mesozoic trend of increase. General linear models predicting these counts using geological age across the Mesozoic have non-significant slopes ( Table 2 ). Taxon counts for the first time bin (Tr1) in Asia and Africa, and for the last two time bins (K7, K8; S1 Table ) in North America are identified as influential data points with high leverage. Significant negative slopes are obtained only when the influential Tr1 data points are excluded from the analysis on their own, and not when all influential data points are excluded together ( Table 2 ). The absence of a trend of increasing regional taxon counts through the Mesozoic is consistent with the observation that most continental regions lack any Late Cretaceous increase ( Fig 2 ), with the exception of North America, where Campanian and Maastrichtian deposits are disproportionately well sampled, with approximately six to eight times as many collections as are known from the most intensively sampled time intervals in regions outside of North America, or two to three times as many collections as the most highly sampled other North American intervals ( Fig 2C ).

General linear models assuming a negative binomial error distribution and ln() link function were used to predict global family counts from geological age across the entire Mesozoic. We found a statistically significant, negative slope that is robust to the exclusion of influential data points ( Table 1 ), indicating a long-term trend of increasing family counts through time. By contrast, statistically significant trends in genus and species counts are only supported if the first Triassic time bin (Tr1; S1 Table ), an influential data point with high leverage, is excluded ( Table 1 ). Furthermore, the significance of this increase is largely due to the occurrence of high taxon counts in the final two time bins of the Cretaceous (K7 and K8), and the slope becomes marginally non-significant when these time bins are also excluded ( Table 1 ). Notably, Late Triassic and latest Jurassic taxon counts also exceed those of most Cretaceous time bins. These observations indicate either that substantial species and genus diversification occurred only in the latest Mesozoic or that oversampling of latest Mesozoic terrestrial faunas has inflated face-value diversity counts for the Campanian and Maastrichtian. The latter possibility is more consistent with our further analytical results, described below. The observation that counts of lower-level taxa do not show a robust trend of Mesozoic increase ( Fig 1B and 1C ; Table 1 ) refutes the proposition that species counts should reveal the hypothesised underlying exponential nature of diversification more prominently than do studies at higher taxonomic levels [ 30 , 49 , 50 ].

“Face-value” observed counts of genera and species occurring globally within time bins approximating 9 million years (Myr) ( S1 Table ; S1 Appendix ) provide little support for exponential diversity increases during the Mesozoic. These counts resemble previously reported global tetrapod family counts [ 12 , 38 ] in several details, including the occurrence of Paleogene diversity levels that are several times greater than those of most Mesozoic intervals ( Fig 1 ). Furthermore, within the Mesozoic, direct counts of families, genera, and species are all highest in the final two stages, the Campanian and Maastrichtian ( Fig 1 ). Nevertheless, counts of genera and species show different long-term patterns than counts of families.

One possibility is that a taxonomic restructuring of terrestrial ecosystems at the K/Pg boundary rapidly established a new dinosaurian/mammalian diversity equilibrium that substantially exceeded the Mesozoic baseline. Such rapid equilibration could only be possible under strong diversity-dependence of diversification rates (e.g., [ 28 , 46 ]). However, mammals, which have increased proportional representation in Cenozoic ecosystems, possess complex teeth, allowing more precise taxonomic identifications from highly fragmentary material than can be diagnosed in fossils of the other highly diverse extant clades (lissamphibians, squamates, and birds), and potentially also Mesozoic dinosaurs. The relatively greater ability to diagnose mammalian species based on fragmentary fossil remains compared to non-mammalian tetrapods is evident in our results: the ratio of mammalian species to species of lissamphibians plus squamates on Earth today is about 1:3, but our subsampled diversity estimates from fossil data yield a ratio of >5:1 in the Paleocene. This suggests that an increase in the proportion of mammalian species within the total terrestrial tetrapod fauna should result in an increase in apparent species diversity in the fossil record, even in the absence of any change in actual tetrapod diversity. It is therefore possible that this apparent Paleocene increase in diversity at least partly reflects a change in the nature of terrestrial vertebrate taxonomy, and is not necessarily a genuine evolutionary phenomenon.

An abrupt and substantial increase in regional subsampled diversity is apparent in the earliest Paleogene, following the end-Cretaceous mass extinction 66 Mya ( Fig 5A ). This increment cannot be explained by bias due to paleogeographic sample spread, which does not change over the boundary ( Fig 4C ). It results entirely from an increase in mammalian species diversity ( Fig 5B ) [ 62 , 63 ], which is disproportionately large compared to the loss of dinosaur diversity ( Fig 5D ). The diversity of non-dinosaurian, non-mammalian tetrapods (“herps”; Fig 5C ) does not change substantially over the Cretaceous/Paleogene (K/Pg) boundary on the temporal resolution of our study, although a major, short-term K/Pg turnover certainly occurred among herps, including squamates [ 64 ]. Nevertheless, our subsampling results tentatively suggest a doubling of herp diversity around the Jurassic/Cretaceous boundary ( Fig 5 ; S1 Fig ; S1 Appendix ), consistent with patterns of subsampled fossil turtle diversity [ 65 ].

Furthermore, equation A25 of reference [ 61 ] (see Materials and Methods : Raup Equation) gives the expected diversity of a clade after a specified time under specified birth and death rates, conditioned on the observation that the clade survived until that time had elapsed. We assumed that the tetrapod crown group originated 100 Myr earlier in the Late Carboniferous, and then specified a net diversification rate of 0.00335 ln(species)/Myr (conservatively assuming the higher diversification rate implied by a direct bias model), and per-lineage death rates of 0.10, 0.15, 0.20, 0.25, and 0.30 ln(extinctions)/Myr (centred on values estimated for Cenozoic North American mammals [ 22 ]). This gives expected Early Triassic regional diversities of 13.3, 19.2, 25.2, 31.1, and 37.0 species. We do not know the actual (rather than observed or subsampled) regional diversities of any studied intervals. However, these expected values of Late Palaeozoic regional diversity obtained under the estimated Mesozoic net diversification rate are lower than the face-value regional species counts of Tr1 in Asia (121 species), Africa (96 species), and Europe (41 species), and of Tr2 in South America (34 species). The diversity counts for these relatively well-sampled regions are not corrected for the possible existence of multiple chronofaunas that could cause counts to exceed the standing diversity at any single time horizon, and they immediately follow the Permian/Triassic extinction event rather than representing Late Permian diversity.

Our estimated long-term diversification rate of 0.00335 ln(species)/Myr is particularly striking in context of the increase in tetrapod diversity that must have occurred during the c.130 million years prior to our study interval, from the Late Devonian origin of tetrapods to the early Mesozoic, which entailed substantially more than a doubling of diversity (e.g., [ 38 , 58 , 59 ]). This can be demonstrated by approximation, assuming that Late Permian diversity was comparable to Early Triassic diversity, which is estimated from the general linear model of subsampled diversity on time as 2.62 ln(subsampled species). The transition from ln(1) to ln(2.62) over 130 million years implies a Paleozoic long-term net diversification rate of 0.0202 ln(species)/Myr, which would generate more than a 40-fold increase in diversity over 190 Myr of Mesozoic time. This estimate is conservative: it could only increase if Late Permian diversity was higher than that of the Early Triassic, as is possible due to the occurrence of the Permian/Triassic boundary mass extinction event (e.g., [ 60 ]). The overall pattern therefore seems to be one of substantial reductions in the long-term net diversification rate of tetrapods during the Paleozoic and Mesozoic, representing the first 87% of their evolutionary history.

The general linear model using geological time to explain subsampled diversity, pooled across geographic regions for the entire Mesozoic, demonstrates only a very weak, but significant slope ( Fig 5A ; Table 4 ; p = 0.02), indicating a trivial net diversification rate of 0.00335 ln(species)/Myr (±2 standard errors yields 0.00089–0.00581 ln(species)/Myr). This implies an expected increase in species richness of 0.637 ln(species), or 89% over c.190 Myr (±2 standard errors yields 18%–202%; and when within-region geographic spread is considered to be a bias the net diversification rate is reduced to 0.00148 ln(species)/Myr, predicting a diversity increase of 32%; Table 4 ). This expected value is equivalent to less than one net speciation per lineage, and comparable to three standard deviations of the regression residuals (s.d. = 0.28). Therefore, short-term diversity fluctuations and statistical counting error have a similar magnitude to our estimate of the long-term expansion of diversity through the Mesozoic. The failure of short-term and inter-regional diversity variations to sum to a greater long-term change would be direct evidence of diversity dependence if we could demonstrate that the proportion of these short-term variations attributable to counting error was low [ 22 , 57 ].

(A) Subsampled species diversity within continental regions for a quorum of 0.4. The dashed line in A is the general linear model predicting subsampled regional diversity from geological age for the entire Mesozoic, modelling taxon counts as a Gaussian distribution and using a ln() link function (coefficients in Table 4 ). (B–D) Subsampled diversities for mammals (B), herps (C; non-mammalian, non-dinosaurian tetrapods), and dinosaurs (D). (E–G) Subsampling curves for (E) the Triassic–Early Jurassic of North America, Asia, and South Africa, (F) the Jurassic–Cretaceous of North America and Europe, and (G) the Cretaceous–Palaeogene of North America. The vertical, dashed grey lines in E–G indicate the target quorum of 0.4. An asterisk is placed in the same location of plots E–G to aid comparison. The data displayed in this figure can be accessed at http://doi.org/10.5061/dryad.9fr76 [ 48 ].

Regional subsampling results indicate a protracted interval of only limited increases in standing diversity spanning the entire Mesozoic ( Fig 5A ). Similar subsampled diversity estimates were obtained for widely separated time intervals such as the Maastrichtian (72.1–66 mega-annum [Ma]) and Kimmeridgian–Tithonian (157–145 Ma) of Europe, and the Kimmeridgian–Tithonian and Norian (228–208 Ma) of North America ( Fig 5A and 5E–5G ). The sporadic availability of data that is rich enough for rigorous diversity estimates makes it difficult to infer short-term patterns of change in standing diversity, although Cretaceous values seem generally higher than those of the Triassic and Jurassic. Nevertheless, we are able to estimate the resultant long-term net diversification rate using general linear models ( Table 4 ), acknowledging that this represents a simplification of potentially more complex short-term patterns.

The results of our alpha diversity analyses should be treated as “first pass” estimates that should be investigated in more detail by future studies, because (1) we did not apply subsampling methods, (2) we did not consider potential environmental or paleogeographical differences between these localities that might affect diversity counts, (3) we did not study the specimens known from these localities to determine the minimum taxon count based on unreported or undiagnosed material, and (4) we did not quantify biases resulting from the likely increased ability of taxonomists to identify fragmentary specimens belonging to extant clades—a bias that could cause relative underestimation of alpha diversity in the Triassic, preceding the origins of most tetrapod crown groups.

It is notable that the maximal within-locality counts generally occur within those intervals containing the greatest numbers of localities such as the Norian (Triassic 4), Kimmeridgian–Tithonian (Jurassic 6), Campanian and Maastrichtian (Cretaceous 7 and 8), suggesting that the intensity of fossil collection activities plays a role in determining the apparent diversity of local communities sampled in the fossil record. Notably, almost all the localities exhibiting high species richness have been intensively bulk sampled for microvertebrate remains—the highest maximal species richness occurs in the latest Cretaceous (Campanian and Maastrichtian) North American localities, which have been intensively bulk sampled ( Fig 6D ). At present, within-locality taxon counts do not suggest any increase in diversity during the early Cenozoic.

(A) Alpha diversity excluding records that are indeterminate at the species level. (B) Alpha diversity including records that are indeterminate at the species level. (C) Per-interval global locality counts (black). (D) Per-interval global bulk sampled locality counts. In all panels, localities that have not been bulk sampled for microvertebrate remains are shown in grey and localities that have been bulk sampled are shown in red. The data displayed in this figure can be accessed at http://doi.org/10.5061/dryad.9fr76 [ 48 ].

Patterns of tetrapod alpha diversity, measured as counts of taxa found within individual fossil localities, are consistent with slow Mesozoic diversification among non-flying tetrapods. Specimens that are taxonomically determinate at the species level are present in 4,357 Mesozoic localities. Of these, just a handful of localities yield substantially greater species counts than most others, including the Late Triassic Placerias Quarry of North America (e.g., [ 66 ]), Late Jurassic Como Bluff Quarry 9 of North America (e.g., [ 67 ]) and Guimarota Mine of Portugal [ 68 ], and the Late Cretaceous Lull 2 Quarry and Bushy Tailed Blowout [ 69 ] of North America ( Fig 6A ). The rare and sporadic occurrence of maximally-diverse fossil localities presents a challenge concerning our ability to resolve patterns of local diversity in the fossil record. Nevertheless, the diversities of these maximally diverse localities increases approximately 3-fold through the Mesozoic, or 2-fold if specifically indeterminate occurrences, which can represent the occurrences of distinct clades and are therefore relevant to diversity counts, are included ( Fig 6B ). These values are comparable to the diversity increase estimated from regional subsampled diversities.

Near-Stasis in Mesozoic Tetrapod Diversification

The Mesozoic–early Cenozoic has previously been regarded as an episode of unbounded diversification, culminating in substantially increased tetrapod diversity on land [2,10–13]. In contrast to this paradigm, our analyses indicate less than a doubling of tetrapod diversity through the Mesozoic, and imply a near-zero long-term net diversification rate. Substantial increases in regional tetrapod diversity were absent from the entire Mesozoic, both for directly counted and subsampled fossil species (Figs 1C and 5). Furthermore, a possible dramatic expansion of tetrapod diversity occurred in the early Paleogene. Our conclusions are strongest if differences in paleogeographic sample spread among regions and intervals are considered to be a bias, in which case Mesozoic regional tetrapod diversity is estimated as being almost static on long time scales (Table 4). Furthermore, first-pass maximal alpha diversity estimates also indicate slow diversification rates (Fig 6), demonstrating that similar patterns occur at local and regional geographic scales. This is consistent with, though not conclusive for, the hypothesis that ecological constraints within local communities could slow diversity increases and thereby regulate diversity at larger scales [28,70]. Our results do not exclude the possibility that an increase in the number of distinct biogeographic regions due to continental fragmentation during the Cretaceous resulted in a greater increase in global diversity than that seen in regional diversities.

Our estimated long-term net diversification rates of 0.00335 or 0.00148 ln(species)/Myr are 1–2 orders of magnitude less than those reported from studies of extant tetrapods (e.g., [1,6]). This difference is unlikely to be explained by underestimation of absolute biodiversity resulting from the incompleteness of the fossil record: estimated diversification rates rely only on accurate inference of relative, not absolute, changes in diversity through time; although we cannot altogether rule out any contribution of fossil record biases (e.g., the possibility that preservational biases could mask an increase in the diversity of small-bodied taxa). Regardless of fossil record biases, a difference between net diversification rates estimated from fossils and those from extant taxa might be expected, because even the richest phylogenies of living taxa lack information on the contributions of entirely extinct clades to diversity dynamics [23,71,72]. Specifically, the contributions of extinct clades and stem groups to total extinction rates cannot easily be estimated from extant-only datasets, and the proportion of entirely extinct clades is likely to increase systematically further back in time from the present. This should cause over-estimation of net diversification rates within inclusive and ancient clades such as Tetrapoda based on the study of living taxa alone. The discrepancy between Mesozoic diversification rates inferred from fossils and diversification rates inferred from living tetrapod phylogenies could also be explained if Cenozoic diversification rates (which are the primary object of inference from living tetrapod phylogenies) substantially exceeded those of the Mesozoic.

Another explanation is plausible if tetrapod subclades show waxing/waning dynamics, as documented in invertebrate genera and mammalian families [23,73,74]. If the dynamics of subclades were asynchronous, whether this were due to diversity dependent interactions [25,75,76], variable environmental tolerances [77], or stochasticity, then the large diversity increases resulting from the waxing phases leading to speciose modern groups could be balanced on long timescales by the waning dynamics of groups that are extinct or depauperate today. This must have occurred in Cenozoic mammals, which show static and diversity-dependent diversification on large scales [22], which apparently results from a zero-sum game among smaller clades that individually exhibit waxing/waning dynamics [23,25].

Near-stasis in Mesozoic tetrapod diversification could be explained by any of three prominent hypotheses: (1) diversity-dependence of diversification rates, or “equilibrial” models, under which speciation and extinction rates become balanced at equilibrial diversity [15,26,27]; (2) the possibility of relatively stable long-term environments during the Mesozoic, which could lead to nearly static diversity under Vrba’s Turnover Pulse hypothesis [77]; or (3) a “damped exponential” model, in which unconstrained diversification is held in check by frequent, stochastic downwards perturbations [2,37,50]. Determining which of these alternatives, if any, provides a good explanation of the pattern requires further work, and we discuss each of them below.

Diversity dependence. The concept of diversity-dependence has been influential in the development of paleontological studies of diversification [15–17] and received significant recent attention from evolutionary biologists studying extant groups [21,28,32]. Recently, diversity dependence has been statistically demonstrated among those vertebrate and non-vertebrate groups that have sufficiently rich fossil records [15–25]. We note that “equilibrial” diversity is attained when the balance of speciation and extinction rates results in an approximately zero net diversification rate, as seen in our analyses. This need not imply that ecosystems are absolutely “saturated,” with all their niches filled. Slow increases in diversity equilibria are possible under a “damped increase” or similar diversity-dependent models in which ecological constraints impose limits to diversity that are increased by the evolutionary discovery of new niche spaces [70]. In the absence of a complete interval-to-interval data series allowing reliable estimation of short-term diversity changes, we cannot directly demonstrate diversity-dependence of tetrapod diversification rates on land using correlation tests (e.g., [16,17,20]) or other methods [22,25]. Nevertheless, rapid recovery of regional diversities following the Cretaceous/Palaeogene extinction event is predicted by diversity-dependence, and observed in our results. Furthermore, strong correlations with physically limiting variables such as palaeogeographic area can be seen as evidence of diversity-dependence because they demonstrate that standing diversity equilibriates rapidly to the availability of environmental resources [28,46]. Indeed, land area, one example of a potentially limiting environmental resource, is a key variable in MacArthur and Wilson’s equilibrial model of island biogeography [26], which is the foundation of diversity-dependent models in paleobiology and evolution [15,27]. The correlations documented here (Table 4) demonstrate scaling of diversity with geographic area, whether the geographic area spanned by fossil localities results from geographic sampling bias or from actual changes in land mass area. We also document that a substantial decrease in long-term net diversification rates occurred through the Paleozoic–Mesozoic. Under diversity-dependence, this could be explained by the low initial diversity of Paleozoic tetrapods (presumably a single species), which would result in high net diversification rates compared to those of Mesozoic tetrapods. Decelerating diversification rates can also occur under alternative models [78]. Of these, environment-driven bursts of diversification postulated under Vrba’s Turnover Pulse hypothesis [77] are one alternative that is relevant to fossil record studies (i.e., it is not an analytical artefact of analysing phylogenies containing only extant taxa), and long time scales (i.e., it does not invoke short-term limits to speciation such as delayed post-speciation range expansion due to physiological conservatism or reproductive interference [78]).

Environment-driven bursts of speciation. The Turnover Pulse Hypothesis, originally formulated to describe diversification among Neogene mammals in Africa during glacial/interglacial cycles [77], proposes that global climatic forcing influences patterns of diversification. Specifically, the appearance and removal of environmental barriers to species distributions via climate change is a prerequisite of turnover: lineages generally exhibit phenotypic stasis in the absence of turnover (i.e., a punctuated equilibrial mode occurs), and environmental oscillations past a threshold amplitude are therefore necessary to drive evolutionary innovation and diversification [77]. This model could explain the variation in long-term net diversification rates documented here in the absence of diversity-dependent interactions among clades, if Mesozoic climates and environments were relatively stable compared to those of the Paleozoic and K/Pg boundary. Although Mesozoic climates are often inferred to have been relatively stable, the Mesozoic was not free from climatic variation, and witnessed apparently extreme climatic events such as the early Turonian thermal maximum around 93 Mya, as well as subsequent global cooling towards the end of the Cretaceous (~75–66 Mya) [79]. A key question, however, is how the timescale and amplitude of Mesozoic climate oscillations [80,81] compares to those occurring during the glacial/interglacial cycles of icehouse regimes, such as that of the late Palaeozoic (mid Carboniferous–early Late Permian; e.g., [82]), and during abrupt environmental deterioration at the end of the Cretaceous (e.g., [83]). The association between climatic variation and net diversification requires thorough investigation to address the question of whether Mesozoic climatic stability is a plausible explanation of slow net diversification rates.

“Damped exponential” model. The damped exponential model was quantitatively formulated to describe clades whose diversification rates depended not only on within-clade standing diversity but also on the diversities of ecologically similar clades [75,76]. This represents a form of diversity dependence similar to that modelled recently in caniform mammals by reference [25]. Nevertheless, the most frequent and recent model referred to as “damped exponential” is one in which fundamentally expansionist diversification, lacking diversity-dependence, is held in check by stochastic downward perturbations caused by frequent mass extinctions [2,37,50]. The predictions of this model, and its fit to real data, are not well constrained because it has not been subjected to significant quantitative examination. Nevertheless, it could result in small or negative long-term diversification rates across large clades, even if some individual subclades show high positive diversification rates on shorter timescales. The near-zero long-term net diversification rates recovered here during the Mesozoic indicate highly balanced speciation and extinction rates. Unlike the diversity-dependent and environment-driven models, the “damped exponential” model does not invoke any terms that specifically act to regulate diversification rates around zero. Therefore, under the “damped exponential” model, near-zero net diversification rates would be coincidental rather than expected, unless the model was modified in such a way that the timings or magnitudes of downward perturbations were diversity dependent.