Climate plays a crucial role in shaping species distribution and evolution over time. Dr. Elisabeth Vrba’s Resource-Use hypothesis posited that zones at the extremes of temperature and precipitation conditions should host a greater number of climate specialist species than other zones because of higher historical fragmentation. Here, we tested this hypothesis by examining climate-induced fragmentation over the past 5 million years. Our findings revealed that, as stated by Vrba, the number of climate specialist species increases with historical regional climate fragmentation, whereas climate generalist species richness decreases. This relationship is approximately 40% stronger than the correlation between current climate and species richness for climate specialist species and 77% stronger for generalist species. These evidences suggest that the effect of climate historical fragmentation is more significant than that of current climate conditions in explaining mammal biogeography. These results provide empirical support for the role of historical climate fragmentation and physiography in shaping the distribution and evolution of life on Earth.

Climate Data and Classification

In this study, we employed the Köppen-Geiger climate classification to categorize climate zones. This system relies on climatic parameters, specifically monthly mean temperature (ºC) and total precipitation (mm), to define climate types (Beck et al., 2018; Köppen, 1884). Given the close correlation between climate and vegetation, these climate zones tend to align closely with global biome patterns (Belda et al., 2014), providing a proxy for examining how climate shapes biome distributions (Mucina, 2019). The Köppen-Geiger climate classification recognises 23 distinct climate regimes, grouped into five major zones: Tropical, Arid, Temperate, Cold, and Polar (Figure 1A). These zones served as the basis for our analysis of the impact of climate change on environmental fragmentation.

Climate data for the last 5 million years were obtained from the high-resolution paleoclimate emulator, PALEO-PGEM (Holden et al., 2019). This dataset offers monthly climate information at a spatial resolution of 0.5º and temporal resolution of 1,000 years, beginning from the pre-industrial era (ca. 1760). We reclassified the climate data into the five major climate zones (tropical, arid, temperate, cold, and polar) for each 1,000-year interval following the methodology outlined by Beck et al. (2018). To facilitate computational operations, we introduced a "-99" value for missing data and made specific adjustments to the function 'KoppenGeiger.m' (Beck et al., 2018), as communicated by H. Beck (personal communication, December 18, 2021), to align with defined precipitation thresholds: "Pthreshold = 2×MAT if >70% of precipitation falls in winter, Pthreshold = 2×MAT+28 if >70% of precipitation falls in summer, otherwise Pthreshold = 2×MAT+14 (Galván et al., 2023). This change was made to rectify a previous code typo that prevented the accurate assignment of some pixels to their climate zone. "Pthreshold" refers to the precipitation threshold for determining the aridity of a climate zone. Meanwhile, "MAT" corresponds to Mean Annual Temperature.

Geographical Framework

This study was conducted on a global scale, to assess whether similar climate zones behaved consistently across different continents. To facilitate these comparisons, we divided the world into three distinct landmasses, hereafter referred to as Americas, Africa, and Eurasia+Oceania (EurOc). The rationale behind this division was to partition our planet into distinct landmasses, each of which would encompass a tropical zone. Upon delimiting the three main landmasses, the different islands were assigned to the nearest landmass in a straight line. This, in turn, corresponds with other biogeographical criteria based on the similarity of flora and fauna. Thus, the three studied landmasses were established as follows:

·       Americas: This category encompasses continental North, Central, and South America, as well as the Caribbean Islands. In the North (Bering Strait), we have included Aleutian St. Matthew, St. Paul, St. George, and Nunivak Islands. St. Lawrence Island is excluded due to its proximity to Europe. In the West, we encompass the Islands off the Mexican West Coast, the Galápagos Islands, and Easter Island. To the South, the Malvinas Islands are included. In the East, we consider Fernando de Noronha, Atol das Rocas Biological Reserve, and Boi Islands. Greenland is part of this category, while Iceland is excluded.

·       Africa: This category covers continental Africa and Madagascar. In the West, it includes the Canary, Madeira, and Savage Islands, the Cabo Verde archipelago, St. Helena, Tristan da Cunha, and Ascension Islands. In the East, Socotra, Seychelles, and the Mayotte archipelago, Comoros, and Mascarene Islands are encompassed.

·       Eurasia + Oceania: This category comprises continental Eurasia, the Arabian Peninsula, Iceland, St. Lawrence Island, Japan, Philippines, Indonesia Australia, New Zealand, and Papua-New Guinea Islands. All the islands of the Pacific Ocean, including the Hawaii Archipelago, are also included. In the Indian Ocean, we consider the Laccadive, Maldives, and Chagos Islands in the West and Ceylon and the Andaman Islands in the East.

The French Austral and Antarctic Lands islands that are closer to the Antarctic region were excluded from the study.

Measuring Fragmentation

To assess climate zone fragmentation, we used the R package landscapemetrics v1.5.4 (Hesselbarth et al., 2019), employing the equal-area Mollweide projection. We applied the “lsm_p_area” function to calculate the number of fragments within each climate zone for each time interval, classifying them based on their area into four size categories:

-        Small fragments (S): Those with an area of up to 3,000 km2, approximately equivalent to the area of a single pixel under our 0.5º resolution.

-        Medium fragments (M): Those with an area between 3,000 and 30,000 km2.

-        Large fragments (L): Those with an area between 30,000 and 600,000 km2.

-        Extra-large fragments (XL): Those with an area exceeding 600,000 km2.

Upon confirming that the number of fragments in the different climate zones followed a normal distribution but did not meet the assumption of variance homogeneity, we conducted the corresponding Welch One-Way ANOVA tests to determine the significance of the results. Given that we were comparing five climate zones, we applied Bonferroni correction to post-hoc results significance. Statistical analyses were conducted using the R library jmv (v2.3.4;53)

In addition to quantifying the number of fragments within climate zones at each time in our series, we computed several additional measures to assess fragmentation:

-Fragmentation Events: The count of instances when the number of fragments increased compared to the previous point in time.

-Fragmentation Strength: The median number of fragments generated in each fragmentation event.

-Maximum Fragmentation: The highest number of fragments produced in a single fragmentation event.

Higher levels of climatic fragmentation are operationally defined as a prevalence of small (S) and medium (M) fragments, while lower levels of fragmentation are characterised by a greater abundance of larger patches (L and XL).

Fragmentation vs. Richness

To explore the relationship between climate fragmentation and specialist mammal richness we sourced mammal range maps from IUCN polygons (IUCN, 2022). Terrestrial mammal data was downloaded on 24th January 2022, while freshwater mammal data was obtained on 21th September 2022. We imported these range maps in shapefile format into R using the ‘rgdal’ package version 1.5-28 (Bivand et al., 2021). We excluded polygons associated with certain families such as Delphinidae, Iniidae, Phocidae, Phocoenidae, Platanistidae, Trichechidae, and the possibly extinct Lipotidae, due to their predominantly aquatic habits. We further excluded species range polygons with presence values of 3 (“possibly extant”) and 6 (“presence uncertain”), as well as range values of 3 (“introduced”) and 4 (“vagrant”) to retain only reliable natural range data (Miraldo et al., 2016).

Range data for each species were converted into a 0.5º raster using the ‘terra’ R package version 1.5-21 (Hijmans, 2022). Mammal species were classified according to their range into specialists, those species that are restricted to a single climate zone, and generalists, which are found in more than one climate zone. To this end, we considered the current distribution of terrestrial mammal species as a reliable representation of their climatic specificities. We then quantified the richness of specialist and generalist mammal species within each climate zone on every continent. We considered various factors of climate fragmentation, including the total number of fragments categorized by size (S, M, L, and XL), the frequency of fragmentation events, as well as the fragmentation strength and maximum fragmentation within each fragment size, climate zone, and continent. In addition, we calculated the mean annual temperature and mean annual precipitation for each climate zone on each continent.

To explore the relationship between these variables and specialist mammal richness, we employed a generalized linear model (GLM). To refine our model and identify the most influential predictors, we employed a bidirectional stepwise regression. This method systematically evaluates interaction terms, ensuring the final model contains only strong predictors or those involved in substantial interactions (Gelman & Hill, 2006). The stepwise regression process continues until no further terms can enhance the model. The selected variables were subsequently evaluated through significance tests, residual analysis, and sensitivity assessments.