Effects of range-through genera.

Our γ diversity estimation within each stratigraphic stage and land mammal age implicitly assumes that time intervals are independent of each other. This assumption increases the likelihood of underestimating diversity especially in short time intervals with few locations, because these often miss genera that were recorded before and after the interval but not within it. Here, we define range-through genera as those that have an occurrence gap during one or more time intervals but are present in intervals either side of the gap (note that this definition differs from that of ref. 59 and others, in which such genera are termed “Lazarus taxa,” denoted “N_r−”). Given that occurrence gaps may reflect either true fluctuations in diversity or data biases, we counted the number of range-through genera for gaps of either one or two intervals in length (note that one-interval range-through genera as defined here are equivalent to “part-timers” in ref. 59, denoted “N_pt”). Then, as a simple test for underestimation of genus richness by the diversity estimation methods, we investigated whether the total genus number (the sum of observed and range-through genera) was above our γ diversity estimates in each of the time intervals.

We found that the number of range-through genera for gaps of either one or two intervals in length were usually low relative to the number of observed genera (see Fig. S1 for range-through genera with gaps of a single interval). For gaps of a single stratigraphic stage, the median number of range-through genera was 1.5 (maximum = 13) in continents and 0.5 (maximum = 11) in regions. The median number of range-through genera with single-interval gaps in continent-specific land mammal ages was 3 (maximum = 20) for continents and was 1 (maximum = 23) for regions. For gaps lasting two time intervals, the numbers decreased markedly: the median number of range-through genera was 0 (maximum = 2) in stratigraphic stages for the continental data, and 0 (maximum = 3) for the regional data. The median number of range-through genera with gaps lasting two land mammal ages was 3 (maximum = 7) for continents and was 1 (maximum = 5) for regions. Virtually all γ diversity estimates with both Jackknife and Chao1 estimators within a given time interval were above its total genus number for both continental and regional datasets (Fig. S2 K and L). The only exceptions were the stratigraphic stage Langhian and the MN05 zone in Eastern Europe (Fig. S2 K and L), where the number of locations was low and the proportion of range-through was genera high (10 locations with 19 genera observed, and 10 range-through genera in the Langhian and 8 in MN05) (Fig. S1G). However, because the numbers of range-through genera were low, especially in the stratigraphic stages, and because our diversity estimates were greater than the total genus number for almost all time intervals and regions, we argue that our diversity estimation effectively took the relatively weak effects of range-through genera into account.

Effects of different spatial location definitions in North American and European datasets.

Spatial locations in the Eurasian source dataset were the unique spatial fossil sites, but locations in North America were based on fossil sites that were grouped if they were close in space and belonged to the same stratigraphic formation (see ref. 30 for details). To assess the impact of this difference in location definition between our source datasets, we ran a sensitivity analysis across the best-sampled part of Europe (following refs. 36 and 43) and determined the extent to which diversity estimates with first-order Jackknife varied when locations are spatially grouped. We sampled all locations in Western and Eastern Europe into the 1° grid (Fig. 1). This grouping of locations is probably at a much larger spatial scale than the grouping inherent in the North American dataset, so it should represent a robust test case for the effects of spatial grouping on diversity estimation. Up to 31 locations (within a given time interval) were grouped into one grid cell (Fig. S3A), and 8,875 occurrences of 411 genera in 1,222 locations across Western and Eastern Europe in the original location dataset were transformed into 6,208 occurrences in 257 grid cells (of the same genera).

There was a strong and significant correlation between Jackknife γ diversity estimates calculated using the original and grouped locations (Pearson’s r = 0.996, t = 71.9, df = 40, P < 0.001) (Fig. S3B). Temporal diversity patterns for the two European regions (Fig. S3 C and D) were highly robust to spatial grouping (median difference between estimates across stratigraphic stages and MN units: 2 genera; maximum: 12), although the confidence intervals were slightly larger with spatial grouping (mean error: 8 genera with locations, 10 genera with grid cells). Given these results, we conclude that our comparisons of North American and European diversity estimates are not affected by differences in spatial grouping of locations between the source datasets.

Performance assessment with simulations based on present-day data.

We conducted a sensitivity analysis of Jackknife estimator performance with the present-day mammalian occurrence data (see Materials and Methods in the main text and discussion below for data description). To maintain the nonrandom spatial sampling structure of our fossil data in the present-day simulation dataset, we extracted present-day species lists by overlaying the species’ range maps with the actual locations from the Neogene fossil dataset (298 locations in North America and 1,225 in Europe which could be assigned to a single continent-specific land mammal age; compare with Fig. 1). With this dataset, we assessed the effects of varying (i) the number of locations and (ii) the preservation of genera within locations on the estimation of γ diversity. To simulate the variation in number of locations observed in our fossil record, we randomly picked 100 sets of locations within each continent consisting of 10, 50, 100, or 200 locations each. These location numbers reflect typical values from the Neogene fossil data: Stratigraphic stages contained 166 locations on average (minimum: 38 locations; maximum: 485 locations) in Europe, and 47 locations on average in North America (minimum: 21 locations; maximum: 93 locations), whereas the more finely resolved continent-specific time intervals contained fewer locations (European MN units contained 1–190 locations and an average of 86, and North American NALMAs contained 11–48 locations and 25 on average). We also explored preservation of genera, because fossil locations preserve only a subset of the genera that actually lived there. For this exploration, we randomly deleted 20, 40, or 60% of the genus occurrences in each location from our randomized subsets of 100 locations. We then calculated γ diversity with the first-order Jackknife algorithm across the 100 randomizations in each location and preservation set and compared the results with the true number of genera from the present-day range maps.

Our simulations showed a robust performance of γ diversity estimation (Fig. S4). With only 50 locations, a median of 93% of the genus richness in North America and 78% in Europe was estimated (Fig. S4A). Diversity in North America was estimated well with as few as 10 locations (median 81%); in Europe performance with 10 locations was lower (median 58%) but increased to median values of 85% and 89% with 100 and 200 locations, respectively. Random deletion of genus occurrences within each location had comparatively little effect (Fig. S4B). In North America, reducing the genus occurrences by as much as 60% led to virtually the same median estimate as that of the full set (25.98 vs. 25.99 genera). In Europe, a decrease in performance was observed with reduction of genus occurrences, but it was weak (median estimate 74% of genera estimated with 60% of occurrences deleted, vs. 85% performance with all occurrences) (Fig. S4B). Our simulations also showed that even though not all genera known to be present in the focal regions were sampled in the complete location set, we were able to approximate the correct number of genera with a reasonable number of locations (compare the ranges of estimates and the horizontal dotted lines in Fig. S4 C and D).

These results indicate that γ diversity can be well estimated, even when data are much more incomplete than our mammalian fossil data are likely to be. Our analyses of the fossil diversity–productivity relationship focused on the stratigraphic stages, which on average contained more locations than necessary to achieve a performance of above 85% in our simulations. Further, provided enough locations were selected, the effect of undetected or nonpreserved genera on diversity estimation was low, assuming random loss of genera and unbiased preservation rates. Recent evidence shows that these assumptions are reasonable for large datasets such as our mammalian fossil record (60). Although for practical and philosophical reasons we can never be sure that we are estimating the true diversity (39, 40), our evaluation of the diversity estimation clearly indicated that estimating mammalian γ diversity from the fossil occurrence data with a first-order Jackknife approach constituted a robust and appropriate way to quantify regional and continental genus richness.

Accounting for temporal uncertainty in allocation of stratigraphic stages.

We obtained paleo-climatic data for temperature and precipitation from several public datasets, all of which were compilations across many original source studies that derived paleo-climatic estimates from fossil plant communities (18, 45⇓–47). We identified and removed duplicate climatic locations (i.e., those included in several of these compilation datasets) based on the original source literature citations and on spatial and temporal proximity and overlap. All climatic records for the Pliocene were allocated to the stratigraphic stage of the Early Pliocene following previous work (46, 47). Age estimates for the other climatic records were usually given as absolute start and end times of stratigraphic stages (18, 45). Because absolute age definitions for any given stratigraphic time interval can change substantially depending on the stratigraphic reference scheme, we had to assign records from these datasets to our stratigraphic stages in a consistent fashion rather than following the original, disparate stratigraphic schemes in the different compilations, which were often not consistent with each other. We therefore allocated the non-Pliocene climatic records to our stratigraphic stage scheme (Table S1) following two different approaches to account for the uncertainty of date assignment. In both approaches, all records allocated to more than two stratigraphic stages were removed to reduce the temporal uncertainty of paleo-climatic estimates.

Our first approach (hereafter “automatic allocation”) was to assign stages automatically based on the minimum and maximum absolute age estimates given in the source data using a temporal accuracy of 0.3 Mya. To do so, we identified the uppermost and lowermost stage assigned to a climatic record as the most recent stages with start dates smaller or equal to the maximum and minimum estimated age ±0.3 Mya, respectively. The upper- and lowermost stage and all stratigraphic stages between these then were assigned to the climatic record. This first approach assigned many stages to most climatic records, so when we excluded records that were allocated to more than two stages, the automatic allocation dataset was severely reduced. As an alternative second approach (hereafter “stratigraphic allocation”), we assigned the nearest equivalent stage manually based on the stratigraphic reference scheme (or more detailed stratigraphic information if given) in the compilation datasets (18, 45). A table showing the allocated stages for each age range of paleo-climatic records is available online; a link to the files is given in Materials and Methods in the main text. In both approaches we excluded climatic records that had no age information or that were allocated to stages older than Aquitanian or younger than Early Pliocene.

After deriving the NPP from the paleo-climatic estimates (see Materials and Methods in the main text and details below), we evaluated the influence of temporal uncertainty in the assignment of stratigraphic stages on the final NPP estimates for continents and regions. We found that the temporal allocation of stratigraphic stages had little effect on the final NPP estimates, which were very similar when the automatic and stratigraphic allocation methods were compared (correlation across the stratigraphic stages: Pearson’s r = 0.96, t = 10.5, df = 10, P < 0.001) (see Fig. 2 I and J for continental estimates). The mean absolute difference in continental NPP between the two datasets was only 48.1 g dry matter⋅m⁻²⋅y⁻¹ (maximum 191.8 g dry matter⋅m⁻²⋅y⁻¹ in the North American Messinian, where estimates were based on two grid cells) (Fig. 2I). We therefore used the dataset based on manual stratigraphic allocation in the diversity–productivity analyses, because it was based on more original paleo-climatic data points.

Accounting for climatic uncertainty in deriving paleo-climatic estimates for each grid cell.

Averaging paleo-climatic records across whole continents or regions directly does not account for spatial variation in sampling or uncertainty of climatic values. Instead, we averaged paleo-climatic estimates across locations within 1° grid cells taking climatic uncertainty into account, calculated the NPP and its uncertainty within these grid cells, and then averaged NPP estimates across grid cells weighted by the uncertainty (the latter two steps are described in the next section). Almost all paleo-climatic records provided ranges between a minimum and a maximum climatic estimate for a given fossil plant community, and the distribution of possible climatic values usually is assumed to be uniform between these endpoints (18). We used this distribution of climatic estimates, or uncertainty range, to calculate mean climatic estimates within each grid cell (Fig. S5). We generated a distribution of climatic values for a given grid cell by allocating each climatic range between a given minimum and maximum value into bins (bin size for temperature: 0.01 °C; bin size for precipitation: 1 mm/y). The mean climatic value assigned to the cell was the highest-density point of the resulting frequency distribution of binned values (Fig. S5B). To incorporate the variation within the grid cell and the uncertainty associated with the climatic estimates, we also extracted the upper and lower limits of the credibility interval containing 50% of the binned values (Fig. S5B).

In the two datasets, i.e., automatic and stratigraphic stage allocation, between 103 and 129 grid cell-by-stratigraphic stage combinations (in 64–82 grid cells) contained only one record for either temperature or precipitation. In these cases, the mean for the grid cell was calculated as the average of the given minimum and maximum climatic value. To derive a credibility interval around this mean, the limits of the 50%-interval of climatic values were approximated as the mean ± 1 SD (this interval should contain 68.2% of values if data are normally distributed), with the SD calculated as the total range (maximum minus minimum) divided by 4 (because in a normal distribution 95% of values fall into the interval of 4 SD units around the mean). Further, a few climatic records did not supply a range of values but instead supplied point estimates: Between 13 and 16 grid cell-by-stratigraphic stage combinations (in 12–14 grid cells) contained only a single record with a climatic point estimate. In these cases, the SD was set to 1 °C or 120 mm/y to derive credibility intervals as above, based on the approximate median values across all climatic estimates given as ranges in the two datasets.

Main models of the fossil relationship.

We ran separate analyses for continents and regions (see also Materials and Methods in the main text). These spatial units were not simple replicates at different spatial aggregation levels but instead represented the best-available current knowledge at a given spatial scale. Because of the lower data coverage in some regions, such as Caucasus and Eastern North America, we excluded these regions in regional analyses. We used GLMMs with Poisson-distributed errors because the response was a count variable (γ diversity estimates were rounded to integers to meet the model assumptions). We tested different model specifications (following ref. 56) to account for the temporal and spatial data structure (Table S2) and determined the best model based on comparison of the AIC as a random slope and intercept GLMM. This kind of model assumes that random effects act on the variance of fixed-effect estimates but not on their mean values (see Table S2 for model comparisons and formula).

We fitted random slope effects of the midpoint of the stratigraphic stage (age in Mya) within each continent or region and random intercept effects for continent or region identity. Preliminary analyses with generalized least-square models containing autoregressive terms indicated that temporal autocorrelation had no significant effect of on variance heterogeneity and that the fitted random effects in the GLMM should account for the temporal trend and the spatial structure in the data (56). The variable of interest (NPP, log-transformed), as well as the two covariates area (in square kilometers, log-transformed) and duration of the stratigraphic stage (in My, log-transformed), were fitted as fixed effects. To approximate the area for which γ diversity was estimated in each stratigraphic stage, area was measured with a convex hull around all the mammalian locations occurring in each stage (using a Mollweide projection). We accounted for area and stage duration because these parameters are known to influence diversity in fossil samples significantly even after controlling for spatial and temporal sampling biases. Correlation coefficients between the variables fitted as fixed effects were low (all absolute r < 0.34 for continents and r < 0.24 for regions), indicating that collinearity was unlikely to affect the models.

Supplemental models of the fossil relationship.

We fitted the same GLMMs to γ diversity estimates averaged across land mammal ages within each stratigraphic bin to test the effect of estimating fossil mammalian diversity in the relatively long stratigraphic stages as opposed to the finer temporal resolution of the land mammal ages (supplemental models in Table S2). Because of the strong effects of area in our models, we also ran additional models to check that the tectonic extension of Western North America during the Miocene did not drive our results. Between 15 and 5 Mya, Western North America experienced a major extension in area that is estimated to have doubled its east–west extent over that timeframe (61). The area estimates for Western North America were rescaled to approximate the extension process with the following multiplication factors: Aquitanian and Burdigalian 0.5×, Langhian 0.6×, Serravallian 0.7×, Tortonian 0.8×, Messinian 0.9×. We then fitted the GLMMs across stratigraphic stages as before and found that results of the continental analysis did not change qualitatively when area values were adjusted for the tectonic extension of North America (supplemental models in Table S2). The NPP remained as a significant predictor of mammalian γ diversity in the continental model but was nonsignificant in the regional model that accounted for tectonic extension (Table S2). We conclude that the diversity–productivity relationship found in the continental models was robust to the influence of North American tectonic changes on the covariate area, but the effect of area extension on the regional level warrants further study.

Models of the present-day spatial relationship.

The substantial differences between fossil and present-day data, especially in their sampling and underlying timescale, precluded a joint statistical analysis of the diversity–productivity relationship (see Results and Discussion in the main text). Therefore, we compared the present-day estimates with the model predictions from the fossil record only visually (Fig. 3). To evaluate the present-day diversity–productivity relationship across the globe, we followed a purely spatial approach, as commonly applied in macroecology, fitting generalized linear models with a Poisson error distribution across the 1° grid cells (Fig. S7). To check the influence of spatial autocorrelation, we refitted these models including a spatial autocovariate as additional predictor variable, which was calculated from the distance-weighted spatial correlation matrix of the residuals from the initial models that included only the NPP as predictor variable (62). The models controlling for spatial autocorrelation were virtually identical to the models presented in Fig. S7: The present-day diversity–productivity relationship across all grid cells globally was significant and strongly positive (n = 12,730, residual deviance = 6,455, slope estimate = 0.0005, z = 149.6, P < 0.001), as was the relationship across all cells falling within the fossil NPP gradient (n = 4,531, residual deviance = 3,553.9, slope estimate = 0.0006, z = 75.2, P < 0.001). The relationship across only the grid cells within our focal regions was significant but much weaker (n = 1,460, residual deviance = 112, slope estimate < 0.0001, z = 3.2, P < 0.01).