Deconstructing the latitudinal diversity gradient of North American mammals by nominal order

Abstract

North American mammals follow a well-established latitudinal diversity gradient in species richness. However, the degree to which species in different mammal clades follow the same latitudinal gradient—and to which each clade contributes to the pattern observed for all mammals remains unknown. Here, we separate the overall mammalian latitudinal diversity gradient by mammal orders and investigate the impact of climate and topography on the distribution of each major mammal clade. We joined an equal-area grid (100 × 100 km cells) of continental North America embedded with environmental variables (n = 10) with mammalian species ranges (n = 753). We used spatial regression models to quantify the relationship between species richness and latitude for all mammals, all mammals excluding select clades, and for each individual subordinate clade (n = 9). We used multiple linear regression and simultaneous autoregressive regression models to determine which environmental variables best explained patterns of species richness for each mammal order. Whereas North American mammals altogether exhibit a strong latitudinal diversity gradient in species richness, most orders deviate from the species richness pattern observed for all mammals and their gradients are weak or entirely absent. Bats (Chiroptera) exhibit the strongest latitudinal gradient—their removal from the pattern for all mammals substantially weakens the total mammalian gradient, more so than when rodents are removed. Environmental variables explain patterns of species richness well for some clades, but poorly for others. The gradient we observe for North American mammals today is likely a combined product of multiple diversification events, dispersals, and climatic and tectonic histories.

La riqueza de especies de mamíferos en Norte América sigue un gradiente latitudinal de diversidad bien establecido. Sin embargo, se desconoce si la riqueza de especies entre clados de mamíferos sigue el mismo patrón latitudinal. También se desconoce la contribución individuales de los clados al patrón observado de la riqueza de especies. En este artículo, se separa el gradiente de diversidad latitudinal general de mamíferos por órdenes de mamíferos y se investiga el impacto del clima y la topografía en la distribución de cada clado principal de mamíferos. La riqueza de especies de mamíferos se obtuvo mediante la sobreposición de sus áreas de distribución (n = 753) sobre una cuadricula (celdas de 100 × 100 km) integrada con variables ambientales (n = 10) de América del Norte continental. Se utilizan modelos de regresión espacial para cuantificar la relación entre la riqueza de especies y la latitud para todos los mamíferos. Este procedimiento se repitió para todos los mamíferos excluyendo clados seleccionados y para cada clado subordinado individual (n = 9). Se utilizaron modelos de regresión lineal múltiple y modelos de regresión autorregresivos para determinar que combinación de variables ambientales explicaban mejor los patrones de riqueza para cada orden de mamíferos. Mientras que los mamíferos de América del Norte en conjunto exhiben un fuerte gradiente de diversidad latitudinal en la riqueza de especies, la mayoría de los órdenes se desvían del patrón de riqueza de especies observado para todos los mamíferos y sus gradientes son débiles o están completamente ausentes. Nuestros resultados muestran que mientras los mamíferos de América del Norte en su conjunto exhiben un fuerte gradiente de diversidad latitudinal en la riqueza de especies, la mayoría de los órdenes se desvían del patrón de riqueza observado para todos los mamíferos y sus gradientes son débiles o están completamente ausentes. Los murciélagos (orden Chiroptera) exhiben el gradiente latitudinal más fuerte: su remoción del patron debilita substancialmente el gradiente latitudinal observado de riqueza de todos los mamíferos, es más acentuada que cuando se remueve el orden de roedores. Las variables ambientales explican bien los patrones de riqueza de especies para algunos clados, pero debilmente para otros. El gradiente observado en la actualidad para los mamíferos de América del Norte es probablemente el resultado combinado de múltiples eventos de diversificación, dispersiones e historias climáticas y tectónicas.

The latitudinal diversity gradient is the most pervasive biogeographic pattern on Earth today (Pianka 1966; Willig et al. 2003; Hillebrand 2004; Fine 2015). This pattern, which describes the increase in species richness with declining latitude, has been documented in numerous clades of vertebrates (e.g., Simpson 1964), invertebrates (e.g., Rex et al. 1993), plants (e.g., Currie and Paquin 1987), and bacteria (e.g., Fuhrman et al. 2008). However, the main driver(s) of this phenomenon—whether biotic, abiotic, or a combination of the two—remains unclear (Rohde 1992; Blackburn and Gaston 1996, Stevens et al. 2019). Hypotheses concerning the origin of the latitudinal diversity gradient generally relate species richness to latitudinal variation in three different categories: ecological limits; diversification rates; and time for species accumulation (Fine 2015; Pontarp et al. 2019). The geographic distribution of water–energy availability is a strong predictor of mammal richness (Currie 1991; Kerr and Packer 1997; Hawkins et al. 2003). Energy is thought to limit species richness at far northern latitudes, whereas water availability is thought to be limiting at southern latitudes (Hawkins et al. 2003; Barreto et al. 2021). These hypotheses propose that regions with greater productivity (e.g., the tropics) are able to support a higher number of individuals and therefore species (Currie et al. 2004; Pontarp et al. 2019). Buckley et al. (2010) reexamined this climate-richness relationship and instead proposed that the ecological limits of clades are evolutionarily conserved (see also Barreto et al. 2019). They suggested that there were more mammal clades of tropical origin during periods of warming during the Eocene when tropical climates were more pervasive, contributing to the large proportion of clades whose richness peaks in the tropics today. In contrast, the frequency of clades with temperate richness peaks increased during periods of cooling following the Early Eocene Climatic Optimum. This suggests that phylogenetic niche conservatism is a major driver of global mammal distributions (Buckley et al. 2010; Barreto et al. 2021). However, the temperature at which mammal families originated explains very little (<3%) variation among clades in the slope of richness gradients, casting doubt on the degree to which niche conservatism underpins latitudinal diversity gradients (Boucher-Lalonde et al. 2015). Others have argued that the higher temperatures seen in the tropics contribute to greater evolutionary rates that increase speciation through shorter generation times, higher mutation rates, and accelerated selection (Rohde 1992; Rolland et al. 2014). Furthermore, the time-for-speciation effect has been cited as a driver of tropical richness, reasoning that older communities have greater species richness because there has been more time for speciation to occur (Willig et al. 2003; Fine 2015). However, some research suggests evolutionary time is more influential in shaping functional diversity rather than species richness (Oliveira et al. 2016).

A common weakness of many of these hypotheses is that they tend to assume that species are equivalent. Yet, regions with equivalent species richness are likely to exhibit substantial ecological differences if, for example, one is comprised of species from the same taxonomic group whereas another is comprised of species from multiple taxonomic groups (Kaufman 1995; Blackburn and Gaston 1996; Marquet et al. 2004). Despite following the same general pattern of increasing species richness closer to the equator, taxonomic groups have been shown to vary widely in the strength of the slope of their gradients, even across similar climatic gradients (Hillebrand 2004). This variation suggests the existence of traits and/or processes specific to individual clades that impact the degree to which they follow latitudinal diversity gradients (Henriques-Silva et al. 2019; Pontarp et al. 2019).

To elucidate what these traits and/or processes might be, we can evaluate whether biogeographic patterns observed at higher taxonomic levels are maintained within that clade at lower taxonomic levels. Relatively few studies have done this, but those investigating species richness patterns in squamate reptiles (Whiting and Fox 2021), trees ((Weiser et al. 2018), rodents (Maestri and Patterson 2016), bats (Stevens 2004), bivalves (Crame 2000), and planktonic crustaceans (Turner 1981) have all shown that the species richness of just one or a small subset of taxa tends to drive the species richness pattern observed for the overall higher clade. Identifying and comparing differing patterns of species richness within the same clade of closely related species may provide critical insights as to which ecological and evolutionary mechanisms bring about certain species richness patterns (Kaufman 1995; Marquet et al. 2004; Stevens 2004). Therefore, deconstructing patterns of biodiversity into smaller taxonomic units can provide the basis for formulating and testing more targeted hypotheses surrounding the origins of the latitudinal diversity gradient and other macroecological phenomena.

North American mammals have repeatedly been shown to exhibit a strong latitudinal gradient in species richness (Simpson 1964; Kaufman and Willig 1998; Badgley and Fox 2000). Within North American mammals, bats exhibit an exceptionally strong latitudinal diversity gradient (Willig and Selcer 1989; Ramos Pereira and Palmeirim 2013; Alroy 2019). The pattern is so strong that many researchers have separated bats from nonvolant, quadrupedal mammals when studying the overall mammalian gradient, suggesting that the gradient for bats drives the pattern observed for mammals as a whole (Wilson 1974; Kaufman 1995; Kaufman and Willig 1998; Rodríguez and Arita 2004; Stevens 2004). Yet, beyond the separation of bats, the latitudinal diversity gradient for North American mammals has yet to be studied comprehensively at the ordinal level (Stevens 2004). Kaufman (1995) established that mammal ordinal richness follows a latitudinal gradient (i.e., the number of mammal orders increases toward the equator) and that richness generally peaks in the tropics within mammal orders. However, comparisons among orders were not explored quantitatively. Among-order comparisons deserve attention because species richness varies substantially among the nominal orders of North American mammals. Rodents are the most species-rich clade of mammals in North America, comprising ca. 51% of North American mammal species, whereas bats comprise ca. 24% of total North American mammalian species richness. Yet, it is unclear how, if at all, the disproportionate richness of rodents contributes to the latitudinal gradient for mammals as a whole. The remaining ca. 25% of North American mammal diversity is divided unevenly among eight other orders, but species richness patterns have yet to be compared across orders and with the overall mammal species richness pattern.

Our study advances prior work by: (1) analyzing the relationship between species richness and latitude for each of the 10 nominal orders of North American terrestrial mammals; and (2) evaluating the degree to which climatic and topographic gradients may have been responsible for driving these underlying patterns. For our study, we first assembled distributional data for mammals inhabiting continental North America and estimated the number of co-occurring species within equal-area grid cells of 100 × 100 km. We then applied spatial autoregressive models to evaluate the effects of individual climatic and topographic variables in driving the patterns of mammal species richness while accounting for the spatial structure in the data.

By limiting our analysis to continental North America, we contribute to a long history of studying patterns of North American mammal diversity (Simpson 1964; Wilson 1974; Badgley and Fox 2000; Qian et al. 2009), but we neglect to consider South American mammal species richness patterns. The evolutionary histories of North and South America are closely related as they have been contiguous land masses for several million years (Montes et al. 2015), permitting dispersal between the continents (Marshall et al. 1982; Woodburne 2010; Bacon et al. 2015). Several authors caution against studying North America in isolation, noting that North America does not extend far enough equatorially to completely capture patterns of low-latitude diversity (Willig and Selcer 1989; Kaufman 1995; Kaufman and Willig 1998). They argue that this reduces the ability of such analyses to generalize about latitudinal trends because they have not been corroborated by comparable data from southern latitudes. Yet, these studies demonstrated that the latitudinal gradient in species richness for mammals as a whole (Kaufman 1995; Kaufman and Willig 1998) and for bats (Willig and Selcer 1989; Alroy 2019) holds when each continent is studied in isolation, when they are studied together, and when continental area is controlled. Therefore, we are confident that our focus on North American mammals will still provide valuable insight into the processes that shape the distribution of mammal species.

Our study addresses the following key questions: (1) To what degree do Chiroptera and/or Rodentia drive the latitudinal diversity gradient for all North American mammals; (2) Does the geographic distribution of species richness for each mammal order mirror the pattern observed for all mammals; (3) How much variation in species richness for each mammal order can be explained by climate and/or topography? We predict that: (A) the North American latitudinal diversity gradient is strongly influenced by Chiroptera and Rodentia because they are the most species-rich orders; (B) different clades will exhibit different geographic distributions of species richness, attributable to differences in evolutionary histories, environmental tolerances, and/or tectonic history and landscape evolution; and (C) climate and topography explain most of the variation in species richness for individual North American mammalian orders. However, we also expect the strength and direction of the relationships—between these variables and species richness—to vary among mammal clades.

Materials and Methods

Mammal distributional data

We considered 753 extant North American mammal species spanning the 10 nominal orders present in North America (Table 1), excluding introduced species (e.g., Myocastor coypus) and island endemics (e.g., Urocyon littoralis). We compiled geographic range shapefiles for each species from Mammals of the Western Hemisphere Version 3.0, published by NatureServe (Patterson et al. 2007). For species recognized by the American Society of Mammalogists Mammal Diversity Database (MDD) but not included in the NatureServe file, we obtained additional geographic range shapefiles from the IUCN Red List of Threatened Species (IUCN 2019), when possible. The taxonomy in our study reflects the MDD v1.0 (Burgin et al. 2018). Since the time of data collection and initial analysis, newer versions of the MDD have been published. Whereas the newer MDD v1.9 recognizes an additional 96 extant species globally, of which approximately 27 have Neotropical and/or Nearctic ranges (Mammal Diversity Database 2022), these tend to be cryptic species based on gene sequence divergences among populations previously assigned to recognized species (e.g., Peromyscus amplus; Hernández-Canchola et al. 2022). Because reliable geographic ranges for many newer species are not available, we used the MDD v1.0 data set in which all species have published geographic ranges. In fact, the most up-to-date geographic range data set for mammals globally (Marsh et al. 2022) uses the taxonomy of the MDD v1.2. Moreover, given the number of new species in MDD v1.9 and the apparent small geographic ranges for most of these, the inclusion of the newer species would not substantially alter our results.

We used an Eckert IV 100 × 100 km equal-area grid projection (EPSG: 54012) of continental North America following Whiting and Fox (2021). To reduce bias resulting from the species-area effect, we excluded coastal grid cells with less than 50% land area (Flessa 1975; Badgley and Fox 2000). In total, we sampled 1,990 grid cells. Utilizing both QGIS 2.18.14 and ArcGIS 10.5.1, we joined the geographic range shapefile for each species to the North American grid and generated species counts for all mammals and for each individual mammal order within every grid cell. Here, we define species richness as the total number of unique geographic ranges that co-occur in each grid cell.

Environmental data

We assembled a series of high-resolution (30 arc-seconds) environmental variables (n = 10) from WorldClim (Hijmans et al. 2005), ENVIREM (Title and Bemmels 2018), and CGIAR-CSI (Trabucco and Zomer 2010) that summarize variation in temperature, moisture/water availability, and topography across continental North America (Table 2). These variables have been found to be particularly influential in explaining patterns of mammalian biodiversity in previous studies of mammal biogeography and macroecology (Currie 1991; Badgley and Fox 2000; Tognelli and Kelt 2004; Qian et al. 2009). Mean annual temperature (MAT), annual evapotranspiration (AET), and potential evapotranspiration (PET) are indicators of available, ambient environmental energy (Currie 1991; Kerr and Packer 1997; Qian et al. 2009)—whereas coldest month minimum temperature (CMMinT) and warmest month maximum temperature (WMMaxT) represent annual temperature extremes (Badgley and Fox 2000). Mean annual precipitation (MAP) is a measure of moisture and water availability (Qian et al. 2009). Precipitation seasonality (PS) and temperature seasonality (TS) sort mammals according to the breadth of their environmental tolerances, and are also important in determining the timing of resource availability (Qian et al. 2009; Alroy 2019). Elevation and topographic relief (i.e., the difference between maximum and minimum elevation in a grid cell) describe the geologic history and habitat heterogeneity of a region (Badgley and Fox 2000). Temperature variables were converted to Kelvin so that all temperature values were positive, and all environmental variables were natural log-transformed prior to our analyses. We reduced the dimensionality of the climate variables using principal component analysis (PCA) on the correlation matrix because the units of measurement differ among the variables (Badgley and Fox 2000). For this ordination, we considered all climatic variables (n = 8) from grid cells in which mammals were present (n = 1,990) and excluded the topographic variables (n = 2). From there, we identified and pruned variables that were highly correlated with each other (pairs with r > 0.75).

Statistical analysis

We evaluated the relationship between species richness and latitude for all mammals and for all individual mammal clades using ordinary least-square (OLS) regression models. For taxonomic groups in which the change in species richness with latitude was better fit by a logarithmic curve than a linear best-fit line, we performed the regression using the natural log of the latitude term. Perissodactyla was included in our analyses of all mammals, but was excluded from our individual, clade-level statistical analyses because there is only one endemic perissodactyl species present in North America (Tapirus bairdii; Table 1). We also assessed the relationship between species richness and latitude for all mammals excluding bats, all mammals excluding rodents, and all mammals excluding both rodents and bats to better understand how the addition and removal of species in these particularly species-rich clades affects the overall strength of the mammalian latitudinal gradient. For each OLS model, we only included grid cells in which richness was greater than or equal to one for each clade/group under consideration.

To ascertain how well climate and topography explain the species richness pattern of each mammal group, we constructed multiple linear regression (MLR) models using the pruned set of climate variables, as well as the two topographic variables. Following Badgley and Fox (2000), all species richness values were first natural log-transformed to ensure similar variances among the groups. We instituted two separate modes of model selection to choose the best combination of climatic and topographic variables for each model, following Whiting and Fox (2021). We implemented stepwise selection (n = 10,000 steps) using the stepAIC function from the ‘MASS’ package in R (Venables and Ripley 2002). Stepwise selection combines forward and backward selection, iteratively adding and then subtracting predictors based on the degree to which they improve the model fit, which is assessed through comparison of Akaike information criterion (AIC) values. We also implemented exhaustive best subsets regression from the ‘leaps’ package in R (Lumley 2017). Exhaustive best subsets regression tests all possible combinations of the predictor variables and then selects the best model by comparing the values of a chosen statistical criterion, in this case the coefficient of determination (R²). We compared the best models calculated by each model selection method. Often, the methods produced the same best model (i.e., the model with the lowest AIC value was also the model with the highest R² value). If the best models differed between the methods, we selected the exhaustive best subsets regression model with the fewest number of predictors whose R² value differed from the R² of the best stepwise selection model by less than 1%.

OLS regression underestimates standard errors and inflates type I errors when positive spatial autocorrelation is present in the data, especially at short distance classes (Diniz-Filho et al. 2003; Dale and Fortin 2014). Accordingly, we explored the spatial autocorrelation in our data by plotting Moran’s I correlograms using the residuals from each regression model using the “pgirmess” package (Giraudoux 2018). We then used simultaneous autoregressive (SAR) error models (Dormann et al. 2007; Kissling and Carl 2008; Ver Hoef et al. 2017) to account for any remaining structure in the residuals of each model for each clade/group. A comprehensive description of SAR models is beyond the scope of this paper but is available elsewhere (e.g., Wall 2004; Dormann et al. 2007; Kissling and Carl 2008). SAR error (SAR_err) models are extensions of OLS regression which assume any autoregressive processes occur in the error term, not in the predictor or response variables (Dormann et al. 2007). Therefore, they add an additional term (λWμ) to the typical OLS regression model $(Y=Xβ+ε)$⁠. This additional term (λWμ) characterizes the spatial structure (λW) of the spatially dependent error term (μ) through the use of a spatial weights matrix (W). In other words, the additional λWμ term represents a spatial weights matrix that weighs the neighboring locations or grid cells, where closer neighbors receive more weight (Dale and Fortin 2014). Following standard procedures (see Kissling and Carl 2008) we identified the neighborhood structure of our grid cells based on their latitude and longitude coordinates, quantifying their proximity using Euclidean distances. Then, we assigned spatial weights using row standardization. Finally, we fit the SAR error models using the errorsarlm function in the “spatialreg” package (Bivand et al. 2021).

See Supplementary Data SD1 for our full data availability statement, including links to all data sources. Our complete, reproducible code and species richness data are available in the Supplementary Data SD2 and SD3, respectively. For further discussion of our environmental variable selection criteria following the PCA, see Supplementary Data SD4.

Results

Mammal comparisons

North American mammal species richness is centered in the low-latitude tropics (Fig. 1A). This pattern persists even when rodents are excluded (Fig. 1B). However, removal of bats from the mammalian gradient reveals a different pattern. Nonbat mammal species richness peaks at mid-latitudes, especially in the Cascade Range, the Sierra Nevada, and the southern Rocky Mountains (Fig. 1C). Removal of both rodents and bats from the mammalian gradient illuminates two regions of elevated species richness: the Pacific Northwest and southernmost Central America (Fig. 1D).

The regression models reflect these observations (Table 3). Latitudinal gradients are characterized by a negative relationship between species richness and latitude, which becomes stronger as the slope of the best-fit line becomes steeper (further away from zero). A slope that is relatively flat (near zero) is indicative of a weak or nonexistent latitudinal gradient. We focus on the results of the SAR models because they generally exhibit little to no spatial autocorrelation in the residuals (Supplementary Data SD5) and tend to explain more variation in species richness than their OLS counterparts (Table 3). SAR models generally result in higher R² values relative to OLS models because they consider complex spatial relationships among neighborhoods of grid cells that OLS models are unable to detect and are therefore able to account for additional variation in species richness that is associated with spatial processes, such as dispersal. There are a few clades in our study for which the SAR models explained less variation in species richness (lower R² value) than the OLS models, which we discuss later in the text.

Considered as a whole, mammals exhibit a strong latitudinal gradient as the slope is steeply negative (b = −59.03) and latitude accounts for 86% of variation in species richness (Fig. 1E). In the absence of rodents, the strength of the mammalian latitudinal gradient persists. The slope is still strongly negative (but slightly flatter than for all mammals) and latitude accounts for slightly more variation in species richness (b = −41.12, R² = 0.89; Fig. 1F). In the absence of bats, the strength of the mammal latitudinal gradient declines sharply and is best described by a linear best-fit line rather than a logarithmic model (b = −0.77; R² = 0.85; Fig. 1G). Mammals in the absence of both bats and rodents follow an even weaker latitudinal gradient (b = −0.31, R² = 0.87; Fig. 1H). All four of these SAR models and their respective slope coefficients are statistically significant (P < 0.001; Table 3).

Among-order comparisons

Bats (Chiroptera) exhibit the strongest latitudinal diversity gradient among the North American mammal orders (b = −29.17, R² = 0.90; Fig. 2A). Their species richness is heavily concentrated in the Neotropics and steadily declines with increasing latitude until becoming absent in the northernmost reaches of the continent (Fig. 3A), paralleling the distribution of Mammalia.

Didelphimorphia, Xenarthra, and Primates all have primarily Neotropical distributions, with species richness peaks at the lowest latitudes (Fig. 3B–D) and modest total species richness (Table 1). We recovered a moderate latitudinal diversity gradient for Didelphimorphia (b = −1.30, R² = 0.82; Fig. 2B). In Xenarthra, the SAR model explains a vast majority of variation in species richness (R² = 0.92; Fig. 2C), but the slope is not significantly different from zero, indicating a lack of discernable latitudinal diversity gradient (Table 3). For Primates, the SAR model explains little variation in species richness (R² = 0.36) and the slope becomes closer to zero (b = −0.14), indicating a weaker gradient. Primates were the only mammal order for which the SAR model failed to eliminate residual spatial autocorrelation (Supplementary Data SD5).

Rodentia displays a mid-latitude species richness peak in the American Southwest and the Sierra Nevada (Fig. 3E), rather than a low-latitude, tropical peak as seen in Mammalia, Chiroptera, Didelphimorphia, Xenarthra, and Primates. As a result, the latitudinal diversity gradient for rodents is weak to moderate in strength (b = −0.46, R² = 0.84; Fig. 2E).

Carnivores are broadly and relatively evenly distributed across the continent, lacking a strong species richness peak at any latitude (Fig. 3F). Lagomorpha species richness peaks at mid-to-high latitudes, near the Pacific Northwest (Fig. 3G). Eulipotyphla exhibits a longitudinal, bimodal distribution of species richness, with peaks in the Pacific Northwest and Appalachia (Fig. 3H). Artiodactyla is fairly evenly distributed latitudinally, but is mostly concentrated in the western half of the continent (Fig. 3I). It follows that Carnivora, Lagomorpha, Eulipotyphla, and Artiodactyla all exhibit near zero, but still statistically significant (P < 0.001; Table 3), slopes. Therefore, these clades exhibit exceedingly weak latitudinal diversity gradients that deviate from the pattern observed in all mammals, even though the SAR models were able to explain a large portion of variation in the species richness of each clade (R² ≥ 0.82; Table 3).

Correlation with climate and topography

The results of our initial PCA, which considered all eight climatic variables, indicate that the first two axes account for 84.6% of the variance in the climatic data (Supplementary Data SD4). MAT, CMMinT, PET, WMMaxT, and TS load most highly on PC1 (64.2% of explained variance), whereas PS, MAP, and AET load most highly on PC2 (20.4% of explained variance). We only consider these first two axes, as the remaining principal components explain much less variance. After pruning highly correlated (r > |0.75|) variables, the subsequent PCA on the smaller climate set explains a comparable amount of variance to the PCA on all eight climatic variables; PET and TS load most highly on PC1 (69.3% of explained variance) and MAP loads most highly on PC2 (20.7% of explained variance; Supplementary Data SD4). Therefore, we reduce the total suite of environmental variables for use in our regression models to elevation, topographic relief, MAP, PET, and TS.

The entire reduced suite of environmental variables (elevation, topographic relief, PET, MAP, and TS) best explains the species richness of all mammals, all mammals excluding rodents, and all mammals excluding bats (Table 4). Overall, our SAR models explain a majority of variation in species richness in each of these clades/groups (R² ≥ 0.79; Table 4). When both rodents and bats are excluded, species richness was best explained by elevation, topographic relief, PET, and MAP (but not TS; Table 4). The SAR model explains 83% of variation in species richness of mammals excluding rodents and bats. All coefficients are statistically significant (P < 0.001) with the exception of TS, which is nonsignificant (Table 4). In all four models, PET contributes the most to explaining species richness.

Climate and topography explain the most variation in species richness for Xenarthra (R² = 0.89). In this clade, the best model includes just three variables (elevation, PET, and TS), of which only TS is statistically significant in the SAR model (P < 0.05; Table 4). Most variation in the species richness of Carnivora is explained by the SAR model incorporating all five environmental variables (R² = 0.81), much more than when evaluated by the equivalent nonspatial regression (R² = 0.32; Table 4). For Chiroptera, elevation, topographic relief, PET, MAP, and TS explain 77% of variation in bat species richness when evaluated using a SAR model (Table 4). PET is the strongest contributor to the model out of all of the environmental variables included. We also found strong correlations between climate, topography, and species richness for Artiodactyla (R² = 0.73), Didelphimorphia (R² = 0.71), Rodentia (R² = 0.68), Lagomorpha (R² = 0.67), and Eulipotyphla (R² = 0.67), although the combination of variables that best explains species richness is unique to each clade (Table 4). In Artiodactyla, we find elevation, PET, and MAP are all significant (P < 0.001) predictors of species richness. In Didelphimorphia, the best model includes PET, MAP, and TS, but no measures of topography (i.e., elevation, topographic relief). Of these three variables, only PET is statistically significant (P < 0.001; Table 4). For Rodentia, Lagomorpha, and Eulipotyphla, each best model includes all five environmental variables (Table 4). The SAR models explain less variation in species richness for rodents and didelphimorph marsupials than the nonspatial regressions, whereas the spatial regressions performed better for lagomorphs and euliptoyphlans (Table 4). In Lagomorpha, only elevation, topographic relief, and PET are statistically significant (P < 0.001) in the SAR model, with PET contributing the most. In Eulipotyphla, all variables are statistically significant, but TS contributes the most to explaining species richness, followed by MAP and PET. The SAR model performed poorly when explaining the species richness of Primates (R² = 0.33; Table 4). The best model includes topographic relief, MAP, and TS, although none of these environmental variables were statistically significant and spatial autocorrelation remains in the residuals of the SAR model at the largest distance classes (Supplementary Data SD5).

Discussion

The latitudinal diversity gradient in species richness for Mammalia in North America is strong, consistent with previous research on North American mammals (Simpson 1964; Kaufman 1995; Kaufman and Willig 1998; Badgley and Fox 2000), as well as many other clades of organisms globally (Currie and Paquin 1987; Rex et al. 1993; Fuhrman et al. 2008). Decomposition of the latitudinal diversity gradient for all mammals into nominal orders reveals that the species richness of each clade has its own unique biogeographic pattern. Chiroptera exhibits the strongest gradient, with species richness increasing dramatically with decreasing latitude. When bats are removed, the latitudinal gradient of the remaining mammals weakens substantially, and the pattern broadly resembles the distribution of Rodentia, the most species-rich clade of mammals in North America. Despite their high species richness (ca. 51% of all North American mammals), when rodents alone are excluded, the latitudinal diversity gradient for the remaining mammals is still quite strong, driven mostly by bats. When both rodents and bats are excluded, the subsequent mammalian latitudinal gradient is weak and mostly reflects the Neotropical species richness of Didelphimorphia, Primates, and Xenarthra. Altogether, our findings show that the mammalian latitudinal diversity gradient appears to be primarily driven by the species richness of bats (Chiroptera), which comprise almost a quarter (ca. 24%) of North American mammalian species richness, and does not reflect the multiple, unique species richness patterns exhibited by other subordinate mammal clades.

Out of the nine North American mammal orders analyzed here, only Chiroptera and—to a lesser degree—Didelphimorphia follow the canonical latitudinal diversity gradient displayed by all mammals. The majority of North American mammal orders deviate from this trend, exhibiting a diverse array of species richness patterns, challenging the perception of the latitudinal diversity gradient as a ‘pervasive’ feature of modern biodiversity. Rodentia is one such clade that deviates from the overall mammal trend. Despite following a moderate latitudinal diversity gradient, rodent species richness peaks at mid-latitudes in the American Southwest and the Sierra Nevada. Additionally, the SAR model exploring the relationship between rodent species richness and environmental variables explains less variation in richness than the nonspatial regression alone. When accounting for spatial autocorrelation, spatially structured variables tend to exhibit weaker effects (Lichstein et al. 2002; Tognelli and Kelt 2004). Therefore, this result may suggest that there is a substantial proportion of variation in North American rodent species richness that cannot be explained by spatial processes. Rather, geohistorical and macroevolutionary processes are likely important drivers of observed rodent species richness. The American Southwest is located at the intersection of several major biotic regions (e.g., the Great Plains, Great Basin, Rocky Mountains, Madrean Sky Islands, and the Chihuahuan and Sonoran deserts) and is underlain by the Southern Basin and Range physiographic province. This region is topographically complex and has been tectonically active throughout the Cenozoic (Dickinson 2002, 2006; Bahadori and Holt 2019), which has had a profound effect on the biodiversity of the region (Coblentz and Ritters 2004; Eronen et al. 2015; Badgley et al. 2017). Our results are thus in alignment with both paleontological (Finarelli and Badgley 2010; Badgley and Finarelli 2013; Badgley et al. 2014) and modern biogeographic (Smiley et al. 2020) analyses suggesting that the disproportionately high rodent species richness found in the American West and Southwest stems from intervals of intensified diversification during periods of increased tectonic activity and global warming over the last 25 million years.

The possibility that biodiversity continues to be underestimated in the tropics, as evidenced by the increasing recognition of cryptic species through molecular methods (Isaac et al. 2004; Adams et al. 2014; Gill et al. 2016), may complicate our interpretation of mammalian latitudinal diversity gradients. Newly delineated taxa tend to be small-bodied (Patterson 2000; Parsons et al. 2022), so we might expect clades like Rodentia to be particularly impacted. Recognition of this “hidden” biodiversity could reveal that rodents follow a stronger latitudinal diversity gradient than is currently recognized. However, we believe that cryptic speciation would not have a significant impact on the overall rodent trend for several reasons. The recognition of cryptic species often reveals that newly delineated species occupy smaller, distinct subdivisions of the original species range, rather than coexisting over a broad geographic area (Arbogast et al. 2017; Eme et al. 2018; Robuchon et al. 2019). Whereas this change could increase local species richness within a given grid cell, it is unlikely that it would increase species richness in every grid cell across a wider region, thereby preserving the existing pattern. Additionally, cryptic species are common in temperate latitudes (Bickford et al. 2007), especially for North American rodents (Neiswenter and Riddle 2010, 2011; Riddle et al. 2014). Even if cryptic tropical species richness did exist that has not been fully accounted for, it is plausible that increases in species richness in the south could be balanced out by the presence of cryptic species richness in the north. Furthermore, it would take an inordinate amount of unrecognized cryptic species richness in the tropics to shift the existing rodent pattern to resemble the pattern for bats.

We found that the latitudinal gradient in species richness for Carnivora is especially weak in North America, contrary to a global survey of mammalian species richness that concluded that carnivoran species richness peaks in the tropics (Rolland et al. 2015). Our results and maps reveal that this pattern is actually much subtler and, rather than species richness solely peaking in the tropics, carnivores appear to be more evenly distributed across the continent. It is possible that this pattern is the result of phylogenetic niche conservatism exhibited by two major subclades within Carnivora—Feliformia and Caniformia. Buckley et al. (2010) demonstrated that, globally, feliforms—which originated and diversified in the tropics—exhibit a strong, positive climate-richness relationship (e.g., more species in the tropics); whereas caniforms—which diversified in temperate environments—exhibit a shallow, negative climate-richness relationship (e.g., fewer species in tropics). Therefore, the geographic distribution of these subclades appears to be limited by the environments in which they diversified (Rolland et al. 2014). If this is the case, we would expect to see high caniform species richness at northern temperate latitudes balanced out by high feliform species richness at southern tropical latitudes in North America, thus producing the relatively flattened latitudinal gradient we observe for Carnivora as a whole. We did not plot species richness at the level of suborder, but a global survey of species richness patterns for Felifomia and Caniformia at a much broader spatial scale (482 × 482 km quadrats) provides preliminary support for this prediction (Pedersen et al. 2014).

Lagomorpha species richness is primarily concentrated at temperate latitudes in the American West and Pacific coast, producing a weak latitudinal diversity gradient. Globally, lagomorph richness peaks at northern temperate latitudes (Rolland et al. 2014). Eulipotyphla displays little latitudinal variation in species richness. Rather, their species richness is concentrated in the Pacific Northwest and Appalachia, resulting in a longitudinally bimodal distribution. This is consistent with a previously conducted spatial analysis of shrew species richness in North America, which found that they follow a parabolic latitudinal trend wherein species richness peaks near 48°N and is most strongly correlated with environmental moisture, not temperature (Berman et al. 2007). The latitudinal diversity gradient of Artiodactyla is also weak. This was briefly acknowledged by Simpson (1964), who noted that artiodactyl richness peaks outside the tropics in North America and actually declines toward the equator. Our SAR models were able to account for additional variation in artiodactyl species richness that the nonspatial models were unable to detect. This variation could have originated from biotic processes, which can generate spatial patterns in species distribution data (Gaspard et al. 2019). Particularly for artiodactyls, we suggest that processes such as dispersal (Jacques and Jenks 2010; Killeen et al. 2014; Lutz et al. 2015) and predation by large mammalian carnivores (Ripple and Beschta 2012; Letnic and Ripple 2017) may be influential in determining their biogeographic patterns.

We found that neither latitude nor environmental variables were good predictors of Primate species richness in North America. Even after accounting for spatial autocorrelation using SAR models, substantial spatial structure remained in the residuals. Residual spatial autocorrelation that remains even after spatial regression is performed is sometimes the result of variables that affect species richness at the local level, which only explain variation at short distances and are more difficult to account for at macroscales (Diniz-Filho et al. 2003). Indeed, these results are consistent with a previous study that suggested forest structure—specifically forest canopy height—is the strongest predictor of primate richness (Gouveia et al. 2014). These results further highlight the importance of utilizing spatial regression models when analyzing species distribution data, as the nonspatial models (e.g., OLS models) overestimated the effects of latitude and environment on primate richness.

Similarly, across the extent of its geographic range, Xenarthra is the only clade studied that does not follow a significant latitudinal gradient in species richness. Whereas the SAR model relating species richness and latitude explains most variation in xenarthran species richness (R² = 0.92), latitude itself was not a significant predictor of their distribution in North America. Our understanding of the drivers of xenarthran distribution may be hindered by the loss of ancient diversity in the clade. Xenarthra has the second lowest species richness of the orders studied here (following Perissodactyla), but in the late Pleistocene there would have been at least seven additional genera of Xenarthra inhabiting North America, including giant armadillos and giant ground sloths, some of whose geographic ranges extended as far north as Alaska (Grayson 1991). Had some of these more temperate species not gone extinct, it is possible that we would observe a different biogeographic pattern overall for Xenarthra, which points to the possible role of historical contingency (i.e., extinction patterns) in driving some biogeographic patterns at low taxonomic levels. Still, there is reason to believe that distributions of both modern and extinct xenarthrans would have been impacted by climate. The SAR model relating environmental variables and species richness explains a great deal of variation (R² = 0.89) in xenarthran species richness; TS is the only significant coefficient (P < 0.05; Table 4). Cold temperature extremes have been postulated to limit the distributions of modern armadillos, the xenarthran species with the northernmost range (Taulman and Robbins 2014), as well as to control the distributions of extinct giant ground sloths in North America (McDonald and Jefferson 2008). The ancestors of xenarthrans lost a key protein involved in thermoregulation early during mammalian evolution, contributing to their poor ability to regulate body temperature and overall reduction in metabolic intensity (Gaudry et al. 2017). Consequently, it is possible that the tendency of xenarthrans to follow a latitudinal gradient depends on metabolic costs associated with living in certain habitats, the distribution of which may be spatially structured (Kilpatrick et al. 2006).

An understanding of the evolutionary diversification of each mammal clade alongside past climatic and tectonic events in North America provides the necessary context with which to evaluate their roles in generating the modern mammal latitudinal diversity gradient. Within Chiroptera, the frugivorous and species-rich Phyllostomidae (leaf-nosed bats) are primarily responsible for driving the latitudinal diversity gradient for all New World bats (Stevens 2004; Villalobos and Arita 2010; Alroy 2019)—consistent with the work of Buckley et al. (2010), who found that New World bats are strong contributors to the global climate-richness gradient and exhibit significant tropical niche conservatism (see also Villalobos and Arita 2010). Phyllostomid bats likely originated in South America (Lim 2009) and proceeded to diversify in both North and South America (i.e., the ‘dual centers hypothesis’; Arita et al. 2014; Rojas et al. 2016; López-Aguirre et al. 2018). Prior to the Quaternary, phyllostomid bats experienced intensified diversification, owing to rapid speciation in the subfamily Stenodermatinae during the Miocene ~18 million years ago (Ma; Dumont et al. 2012; Shi and Rabosky 2015; Rojas et al. 2016). Additional diversification events continued into the Pliocene within the genus Sturnira (Velazco and Patterson 2013). Phyllostomids as a whole likely entered North America prior to the closure of the Isthmus of Panama (Arita et al. 2014; López-Aguirre et al. 2018), but the subsequent radiation of Sturnira in Central America during the Quaternary is hypothesized to have been facilitated by the Great American Biotic Interchange (GABI; Marshall et al. 1982; Woodburne 2010; Velazco and Patterson 2013). Stevens (2006) invoked the time-for-speciation hypothesis to explain high species richness in Phyllostomidae given their tropical origin. The subsequent restriction of phyllostomids to low latitudes has been linked to their shift to frugivorous and nectivorous diets, which would have: (1) limited their geographic range to tropical regions where fruiting plants are plentiful; and (2) increased their metabolic rate, requiring the adaptation of precise thermoregulatory strategies that may not easily translate to temperate regions (Stevens 2006; Czenze and Dunbar 2020). The South American clades Didelphimorphia, Primates, and Xenarthra also immigrated to North America during the GABI (Alfaro et al. 2015; Moraes-Barros and Arteaga 2015; Dias and Perini 2018), though there is evidence that some platyrrhine primates had reached North America by the early Miocene (Bloch et al. 2016). These three clades continued to diversify in North America postdispersal, though to a lesser degree in Xenarthra (Jansa et al. 2014; Aristide et al. 2015; Moraes-Barros and Arteaga 2015). Therefore, the contributions of Neotropical clades to the modern mammalian gradient in North America appears to be relatively recent in origin, geologically speaking, with bursts of speciation in Chiroptera beginning in the Neogene (Rojas et al. 2016) and migrations across the Isthmus of Panama associated with the GABI beginning approximately 2.6 Ma, at the start of the Quaternary (Marshall et al. 1982; Woodburne 2010). Following the GABI, North American mammal communities were profoundly restructured by Late Pleistocene extinction events (Graham et al. 1996; Tóth et al. 2019). Unlike other mass extinctions, these Late Pleistocene extinctions disproportionately impacted large-bodied terrestrial mammals, including North American carnivorans (e.g., saber-toothed cats, American cheetah, dire wolf, short-faced bear), artiodactyls (e.g., camels, stag-moose, long-horned bison), xenarthrans (e.g., giant grounds sloths, giant armadillos), and proboscideans (e.g., mammoths, mastodons; Grayson 1991; Stuart 1991). Many of these species are best known from more temperate, northern latitudes (Grayson 1991, 2006). Given that mammalian geographic range size increases with both body size (Brown and Nicoletto 1991; Pagel et al. 1991) and latitude (i.e., Rapoport’s rule; Arita et al. 2005), the extinction of these temperate megafauna would have impacted patterns of species richness across large swaths of North America and therefore contributed in particular to the flattened gradients observed for modern carnivorans and artiodactyls today.

These histories are consistent with previous analyses of latitudinal diversity gradients in deep time, which found that the mammalian latitudinal diversity gradient was absent in western North America during the early Paleocene Epoch (63–58 Ma; Rose et al. 2011) and likely was not fully established at the continental scale until at least 4 Ma (Marcot et al. 2016; but see Jones et al. 2021). Marcot et al. (2016) suggested that the geologically recent appearance of the latitudinal diversity gradient for North American mammals may have been associated with the onset of Northern Hemisphere glaciation in the Late Neogene. Our results do not contradict or falsify this hypothesis, but do suggest an important role for diversification histories within some lower-level clades for establishing continental richness gradients. These diversification histories conceivably may also relate to climatic histories, either directly or indirectly as mediated by the diversification of myriad plant foods, such as suitable fruiting plants in the Neotropics (Rojas et al. 2012; Chen et al. 2017; Sánchez and Giannini 2018).

Whereas the strength of the latitudinal diversity gradient of North American mammals appears to hold at higher taxonomic levels, our analyses reveal that the species richness of many mammal clades is not distributed along a monotonic latitudinal gradient when studied at lower taxonomic levels. We find that bats, which are very species-rich and follow an extremely strong latitudinal diversity gradient, mask multiple ordinal level patterns when they are included in analysis of all North American mammals. The evolutionary radiation of the bat clade largely responsible for driving the North American latitudinal diversity gradient (Phyllostomidae) is geologically recent, suggesting that the latitudinal diversity gradient may not have always been a pervasive feature throughout Earth history (Mannion et al. 2014). Rather, the latitudinal diversity gradient currently exhibited by North American mammals appears to be an epiphenomenon, driven by the diversification histories of the two most species-rich clades (bats and rodents), and amplified and/or modulated by climatic and/or tectonic histories. It is thus a historically contingent phenomenon without a single, simple driver—such as temperature, precipitation, or energy. Future studies on latitudinal diversity gradients should therefore consider multiple environmental variables, regional and/or continental geologic histories, and the evolutionary histories of individual clades, rather than searching for singular drivers of these complex macroecological patterns.

Supplementary Data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—Data availability statement.

Supplementary Data SD2.—Complete R script for all numerical operations, statistical analyses, and resulting plots.

Supplementary Data SD3.—Mammal species richness counts and environmental variable values for each grid cell.

Supplementary Data SD4.—Ordination of climate variables using principal component analysis across continental North America.

Supplementary Data SD5.—Moran’s I correlograms testing for spatial autocorrelation in the residuals of our regressions.

Supplementary Data SD6.—Comparisons of ordinary least-square (OLS) regression and simultaneous autoregressive (SAR) models for all mammals (Mammalia) and for all mammal orders/subgroups relating species richness and latitude.

Supplementary Data SD7.—Comparisons of multiple linear regression (MLR) and simultaneous autoregressive (SAR) models for all mammals (Mammalia) and for all mammal orders/subgroups relating species richness and latitude.

Supplementary Data SD8.—Perissodactyla species richness grid map.

Acknowledgments

The authors thank U-Spatial at the University of Minnesota for providing student access to ArcGIS software. We also thank K. McNulty, J. Cavender-Bares, and S. Broadfoot for their support and helpful discussions and insights. Finally, we thank the anonymous reviewers whose constructive feedback helped improve the manuscript during revisions.

Literature Cited