Introduction
Latitudinal richness gradients (LRGs; i.e., loss of richness from low latitudes to high) are nearly ubiquitous among modern terrestrial organisms, having been observed in many groups, including angiosperms, birds, mammals, insects, and other invertebrates (McCoy and Connor Reference McCoy and Connor1980; Stevens Reference Stevens1989; Currie and Fritz Reference Currie and Fritz1993; Roy et al. Reference Roy, Jablonski, Valentine and Rosenberg1998; Currie et al. Reference Currie, Francis and Kerr1999; Engle and Summers Reference Engle and Summers1999; Condit et al. Reference Condit, Pitman, Chave, Terborgh, Foster, Aguilar, Valencia, Villa, Muller-Landau, Losos and Hubbell2002; Hawkins et al. Reference Hawkins, Field, Cornell, Currie, Guegan, Kaufman, Kerr, Mittelbach, Oberdorff, O’Brien, Porter and Turner2003; Qian et al. Reference Qian, Badgley and Fox2009; Baselga et al. Reference Baselga, Lobo, Svenning, Aragón and Araújo2012; Condamine et al. Reference Condamine, Sperling, Wahlberg, Rasplus and Kergoat2012). Terrestrial LRGs have been studied among both extant and extinct organisms (Jablonski et al. Reference Jablonski, Roy and Valentine2006, Reference Jablonski, Belanger, Berke, Huang, Krug, Roy, Tomašových and Valentine2013; Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014; Mannion et al. Reference Mannion, Upchurch, Benson and Goswami2014). The study of LRGs in deep time (millions of years) is motivated by interest in the origins and evolutionary history of marine and terrestrial macroecological phenomena (Roy et al. Reference Roy, Jablonski, Valentine and Rosenberg1998; Jablonski et al. Reference Jablonski, Roy and Valentine2006, Reference Jablonski, Belanger, Berke, Huang, Krug, Roy, Tomašových and Valentine2013; Rose et al. Reference Rose, Fox, Marcot and Badgley2011) and the influence of climate, tectonics, and macroevolution (Rose et al. Reference Rose, Fox, Marcot and Badgley2011; Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014; Mannion et al. Reference Mannion, Upchurch, Benson and Goswami2014). The fossil record is the result of a natural experiment that potentially records the evolutionary history of LRGs. However, the link between biodiversity change and its driving mechanisms is filtered through the processes of fossilization, which introduce considerable bias (Newell Reference Newell1959; Raup 1972, Reference Raup1975, 1976, Reference Raup1979; Kidwell and Flessa Reference Kidwell and Flessa1996; Alroy et al. Reference Alroy, Marshall, Bambach, Bezusko, Foote, Fürsich, Hansen, Holland, Ivany, Jablonski, Jacobs, Jones, Kosnik, Lidgard, Low, Miller, Novack-Gottshall, Olszewski, Patzkowsky, Raup, Roy, Sepkoski, Sommers, Wagner and Webber2001; Kidwell and Holland Reference Kidwell and Holland2002; Tomašových and Kidwell 2009, Reference Tomašových and Kidwell2010; Alroy Reference Alroy2010b; Benton et al. Reference Benton, Dunhill, Lloyd and Marx2011; Miller et al. Reference Miller, Behrensmeyer, Du, Lyons, Patterson, Tóth, Villaseñor, Kanga and Reed2014). Our understanding of fossil LRGs therefore hinges on the assumption that we can apply methods of sampling correction to reduce the influence of bias, thus interpreting remaining patterns as biologically relevant.
Characteristic terrestrial fossil-record biases include, but are not limited to, bias against species without mineralized skeletons, body and range-size bias, geographic bias in sampling effort, depositional environmental bias, and rock record–volume bias (Newell Reference Newell1959; Raup 1972, Reference Raup1976; Smith Reference Smith2001; Crampton et al. Reference Crampton, Beu, Cooper, Jones, Marshall and Maxwell2003; Benton et al. Reference Benton, Dunhill, Lloyd and Marx2011; Mannion et al. Reference Mannion, Upchurch, Carrano and Barrett2011). At individual fossil localities, biases against species of small body and geographic-range sizes are also well documented (Behrensmeyer et al. Reference Behrensmeyer, Western and Boaz1979, Reference Behrensmeyer, Kidwell and Gastaldo2000; Donovan and Paul Reference Donovan and Paul1998; Kidwell and Flessa Reference Kidwell and Flessa1996; Kidwell and Holland Reference Kidwell and Holland2002; Cooper et al. Reference Cooper, Maxwell, Crampton, Beu, Jones and Marshall2006). Between 5 and 50% of species in small size classes can go missing from the marine fossil record (Cooper et al. Reference Cooper, Maxwell, Crampton, Beu, Jones and Marshall2006). In terrestrial ecosystems, depending on their trophic group, 30 to 80% of small species are lost from death assemblages (Behrensmeyer and Dechant Boaz Reference Behrensmeyer and Dechant Boaz1980; Kidwell and Flessa Reference Kidwell and Flessa1996). Less abundant species (e.g., species higher in the food chain) with small geographic ranges are similarly less likely to enter the fossil record, because they might not occur in appropriate depositional environments or have large enough populations to ensure preservation (Newell Reference Newell1959; Raup 1972, Reference Raup1976; Smith Reference Smith2001; Crampton et al. Reference Crampton, Beu, Cooper, Jones, Marshall and Maxwell2003; Benton et al. Reference Benton, Dunhill, Lloyd and Marx2011; Mannion et al. Reference Mannion, Upchurch, Carrano and Barrett2011).
In North America, the numbers and geographic coverage of terrestrial fossil localities increase toward the present (e.g., from the Miocene to the Pleistocene, the number of high-latitude fossil localities increases dramatically), which is problematic given the positive correlation between species diversity and area (the so called species-area effect) (Brown Reference Brown1995; Rosenzweig Reference Rosenzweig1995; Barnosky et al. Reference Barnosky, Carrasco and Davis2005; Alroy Reference Alroy2010a). Furthermore, regional sampling changes through time (e.g., the early Cenozoic North American mammal record is mostly limited to the Great Plains and Central Lowlands of North America, while the late Cenozoic record is much more widespread), which is problematic due to likely topographic differences and thus biodiversity differences among regions (Davis Reference Davis2005; Badgley Reference Badgley2010). Our interpretation of biodiversity time series is made more difficult by rock-record bias; fossil-bearing outcrop areas are often correlated with terrestrial vertebrate taxonomic diversity (Raup 1972; Smith Reference Smith2001; Crampton et al. Reference Crampton, Beu, Cooper, Jones, Marshall and Maxwell2003; Benton et al. Reference Benton, Dunhill, Lloyd and Marx2011; Mannion et al. Reference Mannion, Upchurch, Carrano and Barrett2011). It therefore behooves us to understand how the many biases typical of the terrestrial fossil record affect our ability to detect changes in macroecological phenomena through time.
In recognition of fossil-record bias, paleoecologists use a variety of compensatory approaches, including rarefaction, subsampling, and model selection (Raup Reference Raup1975; Colwell and Coddington Reference Colwell and Coddington1994; Alroy Reference Alroy1996; Miller and Foote Reference Miller and Foote1996; Alroy Reference Alroy2010a; Rose et al. Reference Rose, Fox, Marcot and Badgley2011; Benson and Mannion Reference Benson and Mannion2012; Valentine et al. Reference Valentine, Jablonski, Krug and Berke2013; Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014, Reference Fraser, Gorelick and Rybczynski2015; Mannion et al. Reference Mannion, Upchurch, Benson and Goswami2014), the last of which does not attempt to remove bias but allows partitioning of variance among biological and bias-related model terms. Although many compensatory methods reduce sample size differences among localities and time periods, the extent to which true biodiversity patterns are retained after fossilization and whether applying rarefaction or modeling approaches allows apparent diversity patterns to be interpreted as true biological signals remain unknown.
Herein, I divide “latitudinal diversity gradient” into two components: LRGs (i.e., change in the number of species) and latitudinal turnover gradients (LTGs; i.e., change in the species composition of communities). An LTG, as defined here and elsewhere (Bowman Reference Bowman1996; Koleff et al. Reference Koleff, Lennon and Gaston2003), refers to the rate of change in community composition or similarity/dissimilarity among sites along a latitudinal axis. The distinction is important, because species composition can turn over entirely with no change in richness and show more rapid change than richness under environmental perturbation (Dornelas et al. Reference Dornelas, Gotelli, McGill, Shimadzu, Moyes, Sievers and Magurran2014). Estimates of community turnover are also more robust to the formation of death assemblages (Tomašových and Kidwell 2009). Although both richness and community composition vary along environmental gradients (Qian and Ricklefs Reference Qian and Ricklefs2007, 2012; Qian et al. Reference Qian, Badgley and Fox2009; Qian, 2012), for the reasons outlined in the two preceding sentences, they are not necessarily correlated.
For the purposes of this paper, I use extant North American mammals as a test case, because they show a negative LRG (Fig. 1A) that is well explained by latitudinal changes in climate (e.g., mean annual precipitation, temperature, and potential evapotranspiration) (McCoy and Connor Reference McCoy and Connor1980; Currie Reference Currie1991; Kaufman Reference Kaufman1995; Badgley and Fox Reference Badgley and Fox2000; Hawkins et al. Reference Hawkins, Field, Cornell, Currie, Guegan, Kaufman, Kerr, Mittelbach, Oberdorff, O’Brien, Porter and Turner2003). North American mammals also show considerable latitudinal change in community composition (Fig. 1B) that is often attributed to environmental filtering (Qian et al. Reference Qian, Badgley and Fox2009; Kent et al. Reference Kent, Bar-Massada and Carmel2011), a process whereby community composition is partly determined by the environmental tolerances of constituent species (Kraft et al. Reference Kraft, Cornwell, Webb and Ackerly2007; Qian et al. Reference Qian, Badgley and Fox2009; Soininen Reference Soininen2010; Peres-Neto et al. Reference Peres-Neto, Leibold and Dray2012). Furthermore, the North American Cenozoic (65 Ma–present) terrestrial mammal fossil record has been the focus of many paleoecological studies (e.g., Alroy Reference Alroy, Koch and Zachos2000; Figueirido Reference Figueirido, Janis, Pérez-Claros, Renzi and Palmqvist2012; Raia Reference Raia, Carotenuto, Passaro, Piras, Fulgione, Werdelin, Saarinen and Fortelius2012a,Reference Raia, Passaro, Fulgione and Carotenutob; Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014). Although I focus on mammals, the method I develop here allows for the rigorous comparison of modern and fossil communities for any terrestrial group with a fossil and modern record.
Figure 1 Modern North American mammal diversity gradients: A, latitudinal richness gradient; B, longitudinal richness gradient; and C, latitudinal and longitudinal gradients in community turnover.
To address whether latitudinal diversity gradients measured using the fossil record might reflect the true gradients from which they are drawn, I developed a simulated fossilization approach. Using the modern mammal record (i.e., spatially referenced geographic-range data), I created simulated fossil localities based on the geographic distribution of real North American mammal fossil localities from the late Cenozoic and introduce body- and geographic-range size bias as well as variations in sampling effort (0, 25, 50, and 75% species loss). I generated fossil localities iteratively and recalculated the slope of the LRG for each sampling effort and the magnitude of the latitudinal turnover gradient (LTG). It is, however, beyond the scope of this paper to address the causal mechanisms for LRGs.
Methods
I downloaded spatially referenced geographic-range data for modern North American mammals from NatureServe Canada (Patterson et al. Reference Patterson, Ceballos, Sechrest, Tognelli, Brooks, Luna, Ortega, Salazar and Young2007). The entire data set for the Western Hemisphere included 1300 species of non-volant Western Hemisphere mammals after the exclusion of a small number of unreadable or corrupted files. Limiting the study to the continental United States and southern Canada (30–50°N) reduced the data set to 370 species.
Measuring the Modern Latitudinal Diversity Gradient
To quantify the slope of the modern North American mammal LRG prior to the creation of simulated fossil localities, I sampled the ranges of extant mammals using a 1° grid (100×100 km at the equator) under a cylindrical equal-area projection. I calculated total richness in each grid cell and regressed richness against latitude to estimate the slope of the modern mammal LRG in the continental United States and southern Canada (30–50°N). I have chosen this latitudinal range due to the high concentration of North American late Cenozoic mammal fossil localities between 30°N and 50°N. When estimating the slope of the LRG, I did not use a method such as generalized least squares regression, which corrects for spatial autocorrelation, because ordinary least squares is an unbiased estimator of regression coefficients (e.g., slope) even if there is a correlation between error terms such as would be introduced by spatial autocorrelation (Dormann et al. Reference Dormann, McPherson, Araújo, Bivand, Bolliger, Carl, Davies, Hirzel, Jetz, Daniel Kissling, Kühn, Ohlemüller, Peres-Neto, Reineking, Schröder, Schurr and Wilson2007). Furthermore, failing to control for spatial autocorrelation does not introduce systematic bias (Diniz-Filho et al. Reference Diniz-Filho, Bini and Hawkins2003).
As a measure of faunal turnover with latitude (herein termed the latitudinal turnover gradient or LTG) of modern North American mammals, I used detrended correspondence analysis (DCA) (Kent et al. Reference Kent, Bar-Massada and Carmel2011; Oksanen et al. Reference Oksanen, Blanchet, Roeland Kindt, Minchin, O’Hara, Simpson, Solymos, Stevens and Wagner2012) and nonmetric multidimensional scaling (NMDS), which were applied to the grid cell by species-occurrence matrix.
DCA is a form of indirect gradient analysis that is often used in ecological studies to quantify the degree of community change along environmental gradients (Ejrnæs Reference Ejrnæs2000) but has been infrequently applied to the fossil record (Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014). DCA is an ordination method; each fossil site with its complement of species occupies a particular position along the first and second axes. Often, the position of each sampled site along the first DCA axis is regressed against the environmental or geographic variable of interest (in this case, latitude) (Bowman Reference Bowman1996). However, I have used the envfit function from the ‘vegan’ R package (Oksanen et al. Reference Oksanen, Blanchet, Roeland Kindt, Minchin, O’Hara, Simpson, Solymos, Stevens and Wagner2012), which fits environmental vectors to the ordination, effectively testing for a linear relationship between the position of each site in ordination space and the environmental or geographic variable. For the envfit method, the dependent variable is a matrix of ordination scores for sites and the independent variable is the environmental variable measured at each site (e.g., mean annual temperature, latitude, longitude). The envfit function then uses a permutation test to assess the significance of environmental variables. The result is a vector or series of vectors that have both direction (the gradient) and length (magnitude of the gradient) that can be plotted in ordination space (e.g., Fig. 1B). The function also calculates how much of the variation among sites is explained by the environmental variable (R 2; herein also referred to as the magnitude or strength of the LTG, which are both used as synonyms for effect size). High values of R 2 and long vectors indicate strong LTGs (Tuomisto and Ruokolainen Reference Tuomisto and Ruokolainen2006). In this context, R 2 values are precisely correlated with the vector lengths, which are scaled according to the correlation between the environmental vectors and the ordination points (Oksanen et al. Reference Oksanen, Blanchet, Roeland Kindt, Minchin, O’Hara, Simpson, Solymos, Stevens and Wagner2012). Slopes are not associated with the fitted environmental vectors as they are with a standard linear regression, and so I do not interpret the slope of the LTG as I do with the LRG. I use the DCA approach as opposed to other methods of measuring community turnover (e.g., calculating faunal similarity or dissimilarity among grid cells or latitudinal bands, distance decay of similarity) (Soininen et al. Reference Soininen, McDonald and Hillebrand2007; Qian et al. Reference Qian, Badgley and Fox2009), because DCA can be applied to species-by-grid cell occurrence matrices (i.e., equally spaced, even coverage) of extant species and species-by-fossil locality occurrence matrices (i.e., unequally spaced, uneven coverage).
NMDS differs from DCA and other ordination methods, such as principle coordinates analysis, in that it uses the rank orders of distances among communities rather than the raw Euclidean distances (Hammer and Harper Reference Hammer and Harper2006). Unlike DCA, NMDS requires a matrix of dissimilarities among sites. Here, I use the Bray-Curtis dissimilarity index. NMDS first positions the dissimilar data according to Euclidean distances and transforms the distances into rank distances. NMDS then creates a new ordination plot and attempts to place the data in a way that preserves the rank distances from the original plot, which is done iteratively. The configuration of the data on the new plot is compared with the original plot based on the rank order of the Euclidean distances in an attempt to reduce the disagreement between them (referred to as stress). As a result, NMDS is often more robust than approaches such as DCA, because it deals more effectively with missing dissimilarities and can accommodate any dissimilarity measure. As with DCA, I employed the envfit function to test for latitudinal turnover in the dissimilarity among sites. R 2 values are interpreted the same as above.
I then tested for the effects of typical fossil-record biases by developing a simulated fossilization resampling approach in R (R Development Core Team 2016) as described below.
Creating Simulated Fossil Localities
To create “fossil localities” from the modern mammal record, I used an iterative, point-sampling approach (R code in Appendix A in the Supplementary Material). I used the latitude and longitude of late Miocene (~7 Ma; early late Hemphillian North American Land Mammal Age) fossil localities as a model. I downloaded the latitudes and longitudes of the 46 early late Hemphillian sites from the MioMAP database (Carrasco et al. Reference Carrasco, Kraatz, Davis and Barnosky2005). My choice to use the late Miocene is derived from my interest in mammal community changes from the mid- and late Miocene of North America, but my simulated fossilization code allows for input of fossil locality data from any time period. Furthermore, late Miocene mammal fossil localities are fairly well distributed across the continental United States. Due to limits on computing time, it is beyond the scope of this paper to simulate site distributions from other Cenozoic epochs. My approach can be easily modified for any continent and any group for which sufficient extant distribution data exist.
To ensure that I generate “fossil localities” with comparable spatial distributions to the North American fossil record for mammals, I fit frequency distributions (normal, γ, or β) to the latitudes and longitudes of the real fossil localities. The best-fit distributions were chosen using the Akaike information criterion. For the late Miocene, the best-fit distributions were β and normal for the latitudes and longitudes, respectively. I then generated random sites using the fitted probability density functions for both the latitude and longitude of the late Miocene sites. My approach of generating random sites has the benefit of allowing me to evaluate the effect of small changes in the position of simulated fossil localities and to capture different taxa across iterations (i.e., a species with small range may never show up in the simulated fossil record if a single set of sites is generated based on the exact location of the late Miocene sites as performed by Marcot et al. [Reference Marcot, Fox and Niebuhr2016]). I then created locality-by-species occurrence matrices using the over function in the ‘sp’ R package to return the species present at each “fossil locality.” I repeated the procedure 10,000 times (see R code in Appendix A in the Supplementary Material).
Simulating Fossil-Record Bias
The point-sampling approach I use herein is inherently biased against species with small geographic ranges (i.e., species with small ranges are less likely to overlap with the point sample or simulated fossil locality). Therefore, with the sampling method I employ here, no special correction is needed to ensure that species with small ranges are sampled less frequently than those with large ranges.
I weighted the probability of North American mammal species appearing in the simulated fossil record based on body size and trophic categories. I used body-mass categories rather than body-mass estimates because I did not want to make the assumption that preservation potential varies continuously with body mass. Populations from different trophic levels exist at different relative abundances (e.g., carnivores are typically less abundant than herbivores in terrestrial African ecosystems) and thus differ in their preservation potential (Peters Reference Peters1983). Herbivores (plus omnivores) and carnivores larger than 15 kg were given a 95% chance of appearing in the “fossil record” (Behrensmeyer and Dechant Boaz Reference Behrensmeyer and Dechant Boaz1980; Kidwell and Flessa Reference Kidwell and Flessa1996; Kidwell and Holland Reference Kidwell and Holland2002). Small herbivores (<15 kg) were assigned a 60% chance of appearing in the fossil record, while small carnivores (<15 kg) were assigned a 21% chance (Behrensmeyer and Dechant Boaz Reference Behrensmeyer and Dechant Boaz1980; Kidwell and Flessa Reference Kidwell and Flessa1996; Kidwell and Holland Reference Kidwell and Holland2002). Species were placed into body-mass categories based on their estimated species average body mass from PanTHERIA (Jones et al. Reference Jones, Bielby, Cardillo, Fritz, O’Dell, Orme, Safi, Sechrest, Boakes, Carbone, Connolly, Cutts, Foster, Grenyer, Habib, Plaster, Price, Rigby, Rist, Teacher, Bininda-Emonds, Gittleman, Mace and Purvis2009) and Elton Traits 1.0 (Wilman et al. Reference Wilman, Belmaker, Simpson, de la Rosa, Rivadeneira and Jetz2014). The probabilities and body-size categories (Appendix B in the Supplementary Material) I have used here are based on field studies of live–dead assemblages in African savannas (Behrensmeyer and Dechant Boaz Reference Behrensmeyer and Dechant Boaz1980).
I then sampled the extant-mammal simulated fossil record, assuming 100, 75, 50, and 25% rates of species preservation. Because the true sampling rate of the fossil record is unknowable, I use sampling rates that fall within the range of preservation rates suggested by Behrensmeyer and Dechant Boaz (Reference Behrensmeyer and Dechant Boaz1980), Kidwell and Flessa (Reference Kidwell and Flessa1996), and Kidwell and Holland (Reference Kidwell and Holland2002). For these analyses, the probability of species being removed or “fossilized” was based on their body-size categories. If a species is not selected for “fossilization,” it is removed from the analysis globally (removed from all sites at which it may have occurred). I repeated the procedure 10,000 times each for a total of 40,000 iterations (see R code in Appendix A in the Supplementary Material). Globally removing species from the analyses rather than removing species on a site-by-site basis has no effect on the outcome of the analyses described herein (40,000 iterations; Figs. S1, S2).
I calculated total richness at each simulated fossil locality both before and after application of nonparametric richness estimation, which was only applied when using a species sampling rate of less than 100%. The comparison of raw and estimated richness was made to determine whether the application of methods designed to estimate richness when sampling is incomplete can increase the accuracy with which the slope of the richness gradient is estimated. Due to computational load, I only applied first-order jackknife estimation from the ‘fossil’ R package (Vavrek Reference Vavrek2012) and recalculated the LRG slope. I did not use shareholder quorum subsampling (Alroy Reference Alroy2010a,Reference Alroyb), because quorums (combined relative abundances of sampled species) of well below 0.4 were required, which is not recommended for the SQS method (Alroy Reference Alroy2010a). Additionally, I used a latitudinal band method similar to Marcot et al. (Reference Marcot, Fox and Niebuhr2016), in which I divided North America into bands at 2° intervals, calculated richness within each band, both with and without first-order jackknife richness estimation, and calculated the slope of the LRG.
I also applied DCA to the fossil locality by species occurrence matrices to estimate the magnitude of the LTG both with and without simulated fossil-record bias. Rarefaction and richness estimation, which are methods developed to improve the comparability of richness estimates among sites, were not applied in this case, because DCA relies on taxonomic compositions rather than differences in richness. I did not use NMDS, because resulting stress levels were typically above 0.3, which means the ordination is a random representation of the data, thus indicating poor performance of the NMDS (Hammer and Harper Reference Hammer and Harper2006).
Because climate is thought to be one of the primary drivers of latitudinal diversity gradients, I also extracted mean annual temperature and mean annual precipitation data from Climate Wizard for the period of 1951–2006 (Mitchell et al. Reference Mitchell and Jones2005; Girvetz et al. Reference Girvetz, Zganjar, Raber, Maurer, Kareiva and Lawler2009). I calculated the slopes of latitudinal temperature and precipitation gradients (slope of the relationship between climate and richness estimated using ordinary least squares) and variance in faunal turnover explained by precipitation and temperature (R 2 values) for each iteration. That is, temperature and precipitation values were extracted at the latitudinal and longitudinal position of each new simulated fossil locality. The slope of the latitudinal temperature and precipitation richness gradients and the variance explained by these gradients (R 2 values) were then re-estimated for each iteration. As a result, there is variation in the estimated relationship between richness and taxonomic turnover with climate among iterations. The estimated slope and R 2 values were then included as independent variables in the statistical analyses described below. Including the slopes of the climate-richness and climate-turnover gradients allowed me to test whether apparent changes in these relationships explained apparent changes in the LRG among iterations (i.e., do observed changes in the latitudinal-richness gradient relate to changes in the position of the fossil sites and resultant change in the estimated diversity–climate gradients?).
Statistical Analysis of Results
I used an information theoretic approach (automated selection of the best-fit generalized linear regression model using dredge from the ‘MuMIn’ R package) (Bartoń Reference Bartoń2013) to model which combination of factors accounts for changes in LRG slope and LTG magnitude among iterations under simulated fossil-record bias. I included the mean latitude and longitude of all simulated localities, slopes of the latitudinal temperature and precipitation gradients, slopes of the temperature and precipitation richness gradients (or magnitude of the temperature and precipitation turnover gradients), and preservation rate as independent variables in all of the full models. Best-fit linear models were selected using AICc. I also partitioned the explained variance using adjusted R 2 values (Oksanen et al. Reference Oksanen, Blanchet, Roeland Kindt, Minchin, O’Hara, Simpson, Solymos, Stevens and Wagner2012), calculating model averaged regression coefficients as the average slope values across all candidate models with ΔAICc of less than 10 and variable importance for all model terms as the sum of Akaike weights over all the candidate models with ΔAICc of less than 10 (Burnham and Anderson Reference Burnham and Anderson2002; Bartoń Reference Bartoń2013). I only investigated main effects to reduce computation time.
Results
I estimated the slope of the modern North American mammal LRG to be −0.68±0.10 (standard error) species per degree of latitude (R 2=0.07, p<0.0001; Fig. 1A) and the longitudinal richness gradient to be −0.96±0.03 species per degree longitude (R 2=0.54, p<0.0001; Fig. 1B) between 30°N and 50°N. Using DCA, I estimated the magnitude of the LTG in the same region as 0.67 (p<0.001) and the longitudinal turnover gradient as 0.86 (p=0.001), indicating strong gradients in community change for modern mammals (Fig. 1C). Using NMDS, I estimated the magnitude of the LTG in the same region as 0.45 (p<0.001) and the longitudinal turnover gradient as 0.85 (p<0.001).
Among iterations of simulated fossilization, there is considerable variation in estimates of LRG slope and LTG. I also found that the variance in LRG slope estimates decreases at low rates of species sampling, likely due to the loss of all or most of the small-bodied, narrow-range species. Most of the variation in slope estimates among iterations is explained by change in the correlation between diversity and temperature (Table 1, Fig. 2A–D). In other words, iterations in which climate is strongly associated with richness and species turnover also show steep LDGs and vice versa. Only a small percentage of the variation in the LRG and LTG was explained by changes in the mean latitudes and longitudes of the simulated fossil localities, indicating that iteratively generating new simulated fossil localities did not impact the findings of this study (Table 1). Additional unexplained variation probably reflects both the fact that latitude is not a perfect surrogate for the climate gradient in North America and the fact that the steepness of climate gradients varies longitudinally in North America (i.e., chance variation related to minor changes in the placement of simulated fossil localities).
Figure 2 Estimated mammal LRG slopes against (A) slope of the temperature–richness gradient and (B) slope of the latitudinal temperature gradient. Correlation of estimated LTGs with estimates of the (C) temperature–turnover gradient and (D) precipitation–turnover gradient for modern North American mammals. A temperature or precipitation-turnover gradient refers to faunal turnover along a spatial gradient in temperature or precipitation, respectively. All estimates are derived from iterations of simulated fossilization with complete species sampling.
Table 1 Best-fit generalized linear models relating estimates of the slope of the latitudinal richness gradient (LRG) and magnitude of the latitudinal turnover gradient (LTG) to metrics for the geographic coverage of simulated fossil localities, the magnitude of climate–diversity gradients, and climate changes across the landscape. Note: The terms “temperature–richness gradient” and “precipitation–richness gradient” refer to the slope of the climate–richness relationship. The terms “temperature–turnover gradient” and “precipitation–turnover gradient” refer to the correlation between species turnover and climate.
When applying a bias against species with small body size, reducing the species sampling rate of the simulated mammal fossil record from 100% to 25% increased the estimated slope of the LRG (Fig. 3A) and increased the proportion of positive slope estimates (i.e., richness increasing with latitude). Species sampling rate is therefore one of the top explanatory variables in the fitted models (Table 2A). However, the slopes of the climate-richness gradients remain the top explanatory variables (Table 2A).
Figure 3 Box plots of estimated North American mammal LRG slope when species sampling is reduced from 100% to 25% with and without the application of first-order jackknife richness estimation. LRG slope estimates (A) without and (B) with first-order jackknife richness estimation. LRG slope estimates when richness is estimated within latitudinal bands (C) without and (D) with first-order jackknife richness estimation. The dashed line represents the true value of the modern North American mammal LRG between 30°N and 50°N.
Table 2 Best-fit generalized linear models relating estimates of the slope of the latitudinal richness gradient (LRG) and magnitude of the latitudinal turnover gradient (LTG) (A) without nonparametric richness estimation and (B) with nonparametric richness estimation to metrics for the species sampling rate of the simulated fossil record, geographic coverage of simulated fossil localities, the magnitude of climate–diversity gradients, and climate changes across the landscape. Note: The terms “temperature–richness gradient” and “precipitation–richness gradient” refer to the slope of the climate–richness relationship. The terms “temperature–turnover gradient” and “precipitation–turnover gradient” refer to the correlation between species turnover and climate.
Using the method of dividing simulated fossil localities among latitudinal bands produced a very high range of slope estimates for the North American mammal LRG. In fact, the majority of estimates fall on the positive size of a flat gradient (Fig. 3C,D). Application of first-order jackknife richness estimation reduced the explanatory power of species sampling rate (Fig. 2B,D, Table 2B), but not enough to ameliorate the apparent flattening of the LRG slope at the lowest rates of species sampling (Fig. 2B).
LTG estimates also appear weaker as species sampling rate is reduced (Fig. 4) but converge on the true value measured for extant mammals before the application of simulated fossilization. Species sampling rate explained a similar amount of model variance for LTGs as LRGs (Table 2A).
Figure 4 Box plots of estimated North American mammal LTG magnitude when species sampling is reduced from 100% to 25%. The dashed line represents the true magnitude of the modern North American mammal LTG between 30°N and 50°N.
Discussion
Interest in the origins of macroecological phenomena such as LRGs in deep time is burgeoning (Mannion et al. Reference Mannion, Upchurch, Benson and Goswami2014). However, the fossil record is characterized by considerable taxonomic, temporal, and spatial bias (Newell Reference Newell1959; Raup 1972, Reference Raup1976; Smith Reference Smith2001; Crampton et al. Reference Crampton, Beu, Cooper, Jones, Marshall and Maxwell2003; Benton et al. Reference Benton, Dunhill, Lloyd and Marx2011; Mannion et al. Reference Mannion, Upchurch, Carrano and Barrett2011). Herein, I used a simulated fossilization approach to test whether estimates of the LRG and LTG of modern North American mammals are robust to loss of small-bodied, small-range species. In general, I found that LRG estimates are highly impacted by loss of diversity during simulated fossilization, particularly at really low sampling rates. That is, the probability of calculating the true slope of the LRG declines with species sampling rate. In many cases, simulated fossil-record bias completely obliterated evidence of the negative LRG of modern North American mammals, a phenomenon that was not ameliorated by the application of nonparametric richness estimation. The LTG showed a similar loss of signal when species sampling rate fell. However, at medium (50%) to high (75%) rates of species sampling, the magnitude of the latitudinal diversity gradient can be reliably estimated. Furthermore, the relationship of richness and turnover with climate remained strong, even when the species sampling rate was low, suggesting that these relationships might be retained in a fossil record that is biased against small-bodied and small-range species.
Modern North American mammals show a significant decline in richness between 30°N and 50°N (Fig. 1A). Peaks in mammalian diversity occur in the high-altitude regions of Arizona and New Mexico (~35°N, 110°W) and at the most southern parts of the Appalachian region (~32°N, 85°W; Georgia and Alabama) (Fig. 1B). Modern mammal richness gradients are also associated with significant latitudinal and longitudinal faunal turnover gradients (Fig. 1C). That is, there is considerable change in richness and community composition on both the south–north and west–east axes. The best-supported published models suggest that climate and productivity are the primary drivers of LRGs (Currie et al. Reference Currie, Francis and Kerr1999; Currie and Fritz Reference Currie and Fritz1993; Engle and Summers Reference Engle and Summers1999; Condit et al. Reference Condit, Pitman, Chave, Terborgh, Foster, Aguilar, Valencia, Villa, Muller-Landau, Losos and Hubbell2002; Hawkins et al. Reference Hawkins, Field, Cornell, Currie, Guegan, Kaufman, Kerr, Mittelbach, Oberdorff, O’Brien, Porter and Turner2003; Qian et al. Reference Qian, Badgley and Fox2009; Baselga et al. Reference Baselga, Lobo, Svenning, Aragón and Araújo2012; Condamine et al. Reference Condamine, Sperling, Wahlberg, Rasplus and Kergoat2012). The longitudinal gradient in North America reflects the east–west elevation gradient and tendency for topographically complex regions to show steeper richness and β-diversity gradients (Badgley and Fox Reference Badgley and Fox2000; Badgley Reference Badgley2010).
On average, simulated fossilization with a species sampling rate of 100% yielded negative estimates of the LRG slope (Fig. 3A,B) and LTG magnitudes indicative of considerable species turnover (Figs. 4, 5). Only small amounts of variation were explained by changes in the mean latitude and longitude of the simulated fossil localities, suggesting that estimation of the LRG slope is not impacted by small changes in the location of fossil localities (Table 1). The majority of variation in estimated LRG slopes and LTG magnitudes was explained by change in the correlation between diversity and temperature across the landscape (Table 1); within iterations where the relationship between diversity and temperature was tightest, I found steeper LRGs and LTGs.
Figure 5 Relationship of species geographic-range size measured as species occupancy or the number of grid cells occupied with the median latitude of each species’ geographic range.
As I reduced the species sampling rate, I retrieved a higher proportion of positive estimates for the LRG slope (Fig. 3A,B). Positive LRG slopes are particularly interesting, because early Cenozoic North American mammals may have shown flat or positive LRG, which is consistent with contemporaneous climate reconstructions (Rose et al. Reference Rose, Fox, Marcot and Badgley2011). Species sampling rate explained a relatively large proportion of the model variance (Table 2A,B), suggesting that fossil-record bias could be one of the processes leading to apparently positive LRGs. However, the majority of model variance was still explained by the tightness of the relationship between richness and climate (Table 2A,B), also suggesting that the diversity–climate relationship is fairly robust to bias against species with small body and geographic-range sizes. Additionally, at medium (50%) to high (75%) rates of fossilization, the true LRG slope is found well within the upper and lower quartiles of the estimated LRG slopes among iterations (Fig. 3A,B). On average, the “latitudinal band” method produced a very high range of slope estimates for the North American mammal LRG, the majority of which were positive, but modeling results were qualitatively similar (Fig. 3C,D).
Flattening of the LRG as species sampling rate falls is likely explained by both Bergmann’s rule, that is, the increase in mean body size with latitude (McNab Reference McNab1979; Blackburn et al. Reference Blackburn, Gaston and Loder1999; Freckleton et al. Reference Freckleton, Harvey and Pagel2003), and Rapoport’s rule, that is, the increase in mean geographic-range size with latitude (Stevens Reference Stevens1989; Kaufman Reference Kaufman1995; Whittaker et al. Reference Whittaker, Willis and Field2001). Between 30°N and 50°N, there is an apparent increase in the mean and variation of mammal geographic-range sizes (measured here as species occupancy or the number of grid cells occupied by a species; Fig. 5), which is most often attributed to loss of productivity toward the poles (Blackburn et al. Reference Blackburn, Gaston and Loder1999; Stevens Reference Stevens1989). Selective removal of small-bodied, narrow-range species from the simulated fossil record therefore means disproportionate loss of low latitude species and a flattening of the estimated LRG slope.
To what extent might the loss of power to detect the true LRG slope of mammals at low rates of preservation be a product of modern mammal distributions? Modern and Quaternary mammal communities may show steeper gradients in body size (Carotenuto et al. Reference Carotenuto, Diniz-Filho and Raia2015) and range size (Veter et al. Reference Veter, DeSantis, Yann, Donohue, Haupt, Corapi, Fathel, Gootee, Loffredo, Romer and Velkovsky2013) than communities in the past; there were likely times in the past when gradients in body and range size were shallower, such as during the early and mid-Cenozoic (except perhaps when sampling is restricted to a particularly small region of North America, as during Paleocene and early Eocene). Thus, estimates of the richness and turnover gradients for terrestrial mammals may have been less biased That is, if richness gradients were less correlated with body and range size, flattening of the richness and turnover gradients may not have occurred to the same degree as I record herein for modern mammals.
Mean estimates of the LTG magnitude also decrease as preservation rates decreases (Fig. 5A,B). For medium (50%) to high (75–100%) rates of fossilization the true LTG magnitude is found well within the upper and lower quartiles of the estimated LTG magnitudes among iterations (Fig. 4). Some estimates of the LTG, excepting a fossilization rate of 25% of species, were slightly higher than estimated for modern mammals prior to the application of simulated fossilization. Slight overestimation of LTG magnitude results from the discontinuity introduced by the creation of simulated fossil localities. As with the LRG, only a small amount of variation in LTG estimates was explained by the mean latitude and longitude of the simulated fossil localities (Table 2A,B). As for the LRG, the majority of model variance was explained by the tightness of the relationship between species turnover and climate (Table 2A,B), suggesting that measuring the LTG in the fossil record might reflect the underlying diversity–climate relationship even when bias against small-bodied, small-range species is fairly high.
I restricted my study to the region of North America between 30°N and 50°N, as have studies of ancient latitudinal diversity gradients (e.g., Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014; Marcot et al. Reference Marcot, Fox and Niebuhr2016). I chose this latitudinal range because: (1) it is one of the best-sampled regions for Cenozoic fossil mammals, (2) very few mammal fossil localities occur above 50°N in North America due to Quaternary glaciation, and (3) far fewer well-sampled localities occur below 30°N due to differences in preservation potential among the tropical and temperate zones. Restricting my analysis of modern mammals to the same latitudinal range as the fossil record therefore increases the comparability of the modern and fossil richness gradients. It is likely that sampling such a narrow range of latitudes, between which the mammalian LRG is comparatively shallow, has reduced the power to detect a negative LRG. However, I show that, regardless of restricting the latitudinal range of the study, on average, I am able to detect a negative LRG for terrestrial mammals, particularly if species sampling is somewhere between 50% and 100% (Fig. 3A,B). I predict, however, that performing the same analyses across a much wider range of latitudes (e.g., 0°N to 70°N) would increase the power to detect the true richness gradient of terrestrial mammals even at low rates of species preservation (e.g., 25%).
There have been several recent attempts to characterize the stability or change in latitudinal diversity gradients through time using the fossil record of Cenozoic North American mammals (e.g., Rose et al. Reference Rose, Fox, Marcot and Badgley2011; Fraser et al. Reference Fraser, Hassall, Gorelick and Rybczynski2014; Marcot et al. Reference Marcot, Fox and Niebuhr2016). Using different methods for characterizing the latitudinal diversity gradient of extinct Cenozoic mammals of North America, Fraser et al. (Reference Fraser, Hassall, Gorelick and Rybczynski2014) and Marcot et al. (Reference Marcot, Fox and Niebuhr2016) both concluded that shallow and steep latitudinal diversity gradients characterize the early and late Cenozoic, respectively. Both also suggest that that global cooling and drying as well as increased polar glaciation are likely drivers of the observed changes in gradient steepness. Given my finding that the true latitudinal diversity gradient for modern North American mammals is detectable even when the preservation rate is fairly low (50% of the fauna), the conclusions Fraser et al. (Reference Fraser, Hassall, Gorelick and Rybczynski2014) and Marcot et al (Reference Marcot, Fox and Niebuhr2016) are not in conflict with my current study. Furthermore, Marcot et al. (Reference Marcot, Fox and Niebuhr2016) replicate the findings of Fraser et al. (Reference Fraser, Hassall, Gorelick and Rybczynski2014), thus suggesting that steep latitudinal diversity gradients are characteristic of late Cenozoic and modern mammalian assemblages and not their early Cenozoic counterparts.
Conclusions
Studying the deep-time origins of macroecological phenomena can help us understand their long-term climatic and macroevolutionary drivers. Given the considerable taxonomic and spatiotemporal bias of the fossil record, we must understand how the inference of macroecological patterns from the fossil record are impacted by typical biases. I used a simulated fossilization approach to show that estimates of the LRG slope and magnitude of the LTG of modern North American mammals are sensitive to the loss of small-bodied, narrow-range species from the fossil record. However, changes in the slope and magnitude estimates were largely explained by differences in the diversity–climate relationship among iterations, suggesting that these relationships may be measurable in the fossil record using the LRG and LTG as surrogates, with the caveat that at least some of the change among time periods might relate to changes in fossil-record quality. Although I have focused solely on North American mammals, I am confident that the effects reported herein are generalizable to other groups of terrestrial vertebrates that show a decline in richness from low latitudes to high as well as numerical abundance of small-bodied, narrow-range species at low latitudes.
The simulated fossilization approach presented herein allows for the direct comparison of modern diversity gradients with gradients measured in the past. My approach allows users to estimate the modern North American mammal diversity gradient using a similar distribution of localities to the fossil record, when removing different percentages of species in the community, and when adjusting the weighted probabilities of different species entering the simulated fossil record. My approach therefore allows users to directly estimate whether a measured fossil diversity gradient might be reflective of biological reality versus fossil-record bias. Direct comparisons between modern and fossil mammal diversity gradients can be made by “fossilizing” modern mammals using the exact distribution of fossil localities or by fitting frequency distributions as done here and by Fraser et al. (Reference Fraser, Hassall, Gorelick and Rybczynski2014) from each time period of interest (removing varying percentages of species from the simulated fossil record) and comparing the simulation results to the estimated fossil gradient using methods such at the single sample t-test. Although the simulated fossil approach presented herein does not account for all possible sources of bias in the fossil record, it is a significant step toward increasing the rigor with which diversity of fossil taxa is being studied.
Acknowledgments
I thank David Polly, John Alroy, Catherine Badgley, and an anonymous reviewer for insightful review of this article. I also thank Christopher Hassall for considerable help with R coding for this project. I thank S. Kathleen Lyons and the Evolution of Terrestrial Ecosystems working group at the Smithsonian National Museum of Natural History for suggestions that led to the improvement of this paper. I would also like to thank Natalia Rybczynski and Root Gorelick for reading and copyediting earlier iterations of this paper. I thank David Currie and Lenore Fahrig for feedback related to this project and for reading an earlier version of this paper. The work was supported by a Natural Science and Engineering Research Council of Canada (NSERC) postgraduate fellowship, an Ontario Graduate Scholarship, and Sir James Lougheed Award of Distinction from the Government of Alberta.
Supplementary materials
Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.8p0s7
