Frontiers of Biogeography Conservation biogeography of the terrestrial mammals in Iran: diversity patterns, and vulnerability to climate change and extinction

Under the framework of a “conservation biogeography” approach, initially, I reviewed and updated the taxonomy and distribution of the rich but understudied mammalian diversity of Iran. This data then formed the basis for the biogeographical regionalization of this complex transitional area using hierarchical clustering and infomap network methods. I used linear models to explore the correlates of extinction risk for this threatened mammalian fauna. Functional grouping of target species was used to assess their vulnerability to the magnitude and velocity of climate change impacts. Both clustering and network methods successfully illuminated the intricate biogeographic patterns, while the network detected many more small bioregions, including two transition zones. The extinction risk analyses revealed that human activities, such as hunting and persecuting (direct impacts) played a major role in the decline of these taxa, as opposed to the minor effect of indirect and intrinsic and extrinsic factors. The magnitude and velocity of climate change impacts varied significantly between functional groups, with the highest risk of exposure to extreme climates in large and threatened species occurring in lowlands. This study provides a foundation for future biogeographic, systematics and ecological studies of Iranian mammals while simultaneously adding to the limited available information on the bioregionalization at regional scales. It also highlights the importance of incorporating threats in extinction risk models and functional trait information in climate change impact assessments.

models to explore the correlates of extinction risk for this threatened mammalian fauna.
Functional grouping of target species was used to assess their vulnerability to the magnitude and velocity of climate change impacts. Both clustering and network methods successfully illuminated the intricate biogeographic patterns, while the network detected many more small bioregions, including two transition zones. The extinction risk analyses revealed that human activities, such as hunting and persecuting (direct impacts) played a major role in the decline of these taxa, as opposed to the minor effect of indirect and intrinsic and extrinsic factors. The magnitude and velocity of climate change impacts varied significantly between functional groups, with the highest risk of exposure to extreme climates in large and threatened species occurring in lowlands. This study provides a foundation for future biogeographic, systematics and ecological studies of Iranian mammals while simultaneously adding to the limited available information on the bioregionalization at regional scales. It also highlights the importance of incorporating threats in extinction risk models and functional trait information in climate change impact assessments.

Introduction
Biogeographical transitional areas are prime targets for conservation biogeography research, both because of their conservation value as they often can serve as biodiversity hotspots and because of their high faunal complexity and diversity, which makes them ideal for biogeographical studies (Ferro andMorrone 2014, Morrone 2018). As a country located just at the meeting point of three main biogeographic regions-Palaearctic, Saharo-Arabian, and Oriental-Iran is a clear example of a biogeographical transition area. Even at the regional scale, Iran acts as a proverbial bridge between the Arabian Peninsula, Mediterranean, Central Asia, and Indian region (Fig. 1). Iran's unique biogeographical position further complicated by extreme landscape heterogeneity and climate diversity makes Iran one of the richest and most complex regions in western Asia from a biodiversity viewpoint, of which mammals are one of the most prominent faunal elements.
The Iranian land mammal fauna offers an ideal setting in which to analyse biogeographical patterns and to investigate the effects of increasing human activities and extreme climates for several reasons. First, the country is rich in terms of mammal diversity as it contains 192 terrestrial species from 34 families, perhaps one of the most diverse countries in the entire Palaearctic realm. Second, its mammal fauna constitutes a complex assemblage, with a mixture of species with different biogeographical affinities and life histories adapted to a wide range of environmental conditions from extremely high and cold mountains to extremely warm deserts (Firouz 2005). Third, the Iranian Plateau acts as a cradle of diversification and center of genetic diversity for some species and a phylogeographic hotspot for others, and is a vast contact zone where different species and subspecies meet (Hardouin et al. 2015, Schai-Braun and Hacklander 2016, Mahmoudi et al. 2017. Fourth, a wide range of high conservation interest species of carnivores and ungulates occur there, and the country serves as the last stronghold for some endangered species in the region; simultaneously, a substantial portion of west Asia's distribution ranges of large mammals are currently confined to its border (Firouz 2005). Fifth, the country's biodiversity is historically threatened by local human activities (Firouz 2005) and now face the mounting threat of climatic change, since the entire Middle East is at high risk of exposure to future extreme climates (Pal and Eltahir 2015). The country therefore deserves to be examined more extensively with current methodological approaches. Nevertheless, the lack of both reliable and adequately detailed taxonomic knowledge and distributional data hampers further progress in analysing biodiversity patterns and their vulnerability to extinction and climate change impacts.
Terrestrial mammals are a major component of global biodiversity, as well as excellent indicators of the global impact of humans (Ceballos et al. 2017), thus making them perfectly suitable for conservation biogeography research. This diverse group displays a considerable variation in biological, ecological, morphological, and physiological characteristics, which allow them to occupy a wide range of climatic and environmental conditions across different biogeographical regions (Davidson et al. 2017). Such variation in their characteristics especially in life history traits are known to influence their susceptibility to extinction and affect their response to environmental changes (Davidson et al. 2017). For instance, it has been recognized that species with certain traits such as low fecundity are more vulnerable to extinction (González-Suárez and Revilla 2013), and species with specialized diets are more likely than others to be affected by climate change (Pacifici et al. 2017). Mammals also respond differently to the direct threat posed by human activities (hunting and persecuting) and indirect threats (habitat destruction and fragmentation) .
Most importantly, this group is of great conservation concern because mammals are one of the most threatened taxonomic groups of animals, with almost a quarter of species-for which there are adequate data available-at risk of global extinction (IUCN 2019). This situation is exacerbated in the case of biodiversity hotspots, where diversity is highest and where most threatened mammals occur (Visconti et al. 2011). Fortunately, increasing availability of biological/ecological traits data-better known in mammals than other taxa-and high-resolution, long-term climatic data enables quantification of their vulnerability to extinction risk and to current and future climate change (Davidson et al. 2017).
Under the framework of a "conservation biogeography" approach, which aims to understand the spatial patterns of diversity and to assess its vulnerability to human activities (Richardson and Whittaker 2010), the purpose of my thesis research (Yusefi 2019a) was to analyse biogeographic patterns in species diversity of terrestrial mammals in Iran and assess the vulnerability of this rich, but understudied, fauna to extinction and climate change. This research, however, has required a sound knowledge of species diversity (i.e., accurate taxonomy of species) and distribution (i.e., accurate description of species distributions; Richardson and Whittaker 201o), which was not available to me. Thus, initially I incorporated all the available data, reviewed, and updated the taxonomy and distribution of all the known species to fill the knowledge gaps in this regard, also known as the Linnean and Wallacean shortfalls (Richardson and Whittaker 2010). This dataset was then used to analyse the biogeographical structure of this complex country, which, interestingly, was the subject of one of the first studies in regional biogeographical analysis (Blanford 1876a; Supplementary   Fig. S1), a study that happened at the same time that Blanford had a heated debate with Wallace over the inclusion of India in the Oriental realm (Blanford 1876b). Mapping biogeographical regions of the country was later attempted by a few others using different taxonomic groups (e.g., Misonne 1959, Zohary 1963, but poor-quality data and lack of quantitative analyses have left biogeographical regions ambiguously defined. In this study I also took advantage of this large dataset on regional mammal distribution from a complex transitional area to assess the applicability of the recently developed network-based infomap algorithm in bioregionalization analysis at regional spatial scales, and I compared the results with the commonly used hierarchical clustering method, which unlike the global/continental scales (e.g., Vilhena and Antonelli 2015, Edler et al. 2017, Droissart et al. 2018) remains unevaluated.
The database was then used to identify the main factors responsible for extinction and population decline. In this extinction risk modelling, I emphasized the potential role of threats, other than traits, because previous studies on extinction risk have mostly focused on species' intrinsic and extrinsic traits, and the role of human threats, especially direct ones, in explaining species extinctions remains largely unknown (Murray et al. 2014, Verde Arregoitia 2016. This subject remains particularly understudied at regional spatial scales (besides being biased to certain geographical areas) even though the regional modelling has advantages, such as revealing patterns that would be masked in large-scale comparisons (Keil et al. 2018) and producing more focused outcomes which could have more influence on conservation practices (Cardillo and Meijaard 2011). Iranian mammals constitute in particular an ideal group to find the effects of human threats because the fauna of Iran still retains almost a complete assemblage of large mammals. Whereas in many other regions/countries most of exploitable mammal species are long extinct, making it hard to link extinction risk with hunting, which could thus lead to overestimation of other predictors (Tingley et al. 2013).
Lastly, I used spatial distribution data of 188 species to evaluate how the predicted magnitude and velocity of climate change in the highly vulnerable region of the southwest Asia (or Middle East) might affect functional groups in mammals. Although climate change is predicted to harshly impact southwest Asia (Pal and Eltahir 2015), the potential impacts of extreme climates on the region's biodiversity remains largely unknown yet. As outlined above, I used functional groups rather than species because evidence confirms that species´ response to climate change, besides the level of exposure, is largely associated to life-history traits (Pacifici et al. 2017). Nevertheless, previous studies of climate change impacts focused mostly on changes in potential climate suitability at species level by constructing individual species distribution models (Pacifici et al. 2015).

The study area
Iran is a sizable country (1.65 million km 2 ) located in southwest Asia, between the latitudes of 25°30' and 40° north and the longitudes of 44° and 63°30' east, with elevation gradients ranging between -60 to 5671 meters. The country is known as the land of extremes, with ground surface temperatures range from −29 °C to 72 °C, and the annual rainfall varies from <100 mm in Central Basin up to >2,000 mm in the Caspian Sea coasts (Iran Meteorological Organization 1 ). The country holds two large mountain ranges, Alborz in the north and Zagros in the west, and a vast Central Basin, which has high physiographic complexity with several scattered large mountains and numerous other small mountains as well as large hyper-arid deserts (Firouz, 2005). The country hosts two (Caucasus and Irano-Anatolian) of the seven global biodiversity hotspots recognized in Asia ) and contains seven global terrestrial biomes and 19 ecoregions (Dinerstein et al. 2017; Fig. 1).

Systematics and distribution: reducing the Linnean and Wallacean shortfalls
Addressing both the diversity and distribution of Iran's mammal fauna has been the objective of researchers for long time, as its unique features have attracted the interest of many researchers since the early naturalists and explorers (Firouz 2005). Therefore, a growing body of taxonomic and distributional data has been accumulated, especially in recent decades owing to intensive surveys and increasing availability of molecular tools. The results, however, are widely scattered in numerous publications and unpublished obscure sources, most of which are not easily accessible. Additionally, uncertainty in species diversity estimates were apparent across different taxa, due to a paucity of systematic surveys and to taxonomic discrepancies (see discussion in Benda et al. 2012). To overcome this drawback and to provide a comprehensive picture of the country's mammalian diversity, I (at best) compiled and reviewed all the current available information on Iranian mammalian diversity and distribution since mid-18 century.
This resulted in the first comprehensive updated species account (Yusefi et al. 2019b and supplementary material), with complete taxonomic review (at the subspecies level) of all known species. Here, I cited over 850 publications, such as journal articles, books, dissertations, documents, and abstracts that had relevant information. The taxonomy was based on Wilson and Reeder (2005) and was updated whenever appropriate to include published findings, especially recent molecular studies. I followed "Handbook of the Mammals of the World" (Volumes: 1,2, 6-9; Wilson and Mittermeier 2009, 2018, 2019Wilson et al. 2016Wilson et al. , 2017 for taxonomic treatment of the species. This review has changed the species composition of Iranian mammal fauna significantly: the scientific names of 45 taxa (mostly rodents and bats) have changed, 13 of them are new species or new records, and 32 had changes in classification or nomenclature when compared to the last review by Karami et al. (2008). Additionally, the distribution maps of these species were determined with detailed accounts of the species occurrences. To do this, I compiled, standardized, and organized in a georeferenced database a total of 14,965 species occurrence records, mainly from published data, online databases (GBIF 2 ), Iran Department of Environment unpublished data, media news (internet websites), unpublished reports, personal communications, and my own field observations (see details in Yusefi et al. 2019b and supplementary data). Most of the reviewed literature and collected records came from sources and databases written in Persian and are therefore quite difficult to access.
The review concludes that knowledge about Iranian mammals, both in diversity and distribution, is far from complete. In particular, shrews (Soricidae), moles (Talpidae), some bats (e.g., Rhinolophidae and Vespertilionidae families), some rodents (e.g., Calomyscidae, Cricetidae, Dipodidae, and Muridae families) and meso-carnivores (particularly, Herpestidae and Mustelidae families) are poorly known. Within these taxa, which represent priority groups for future work, the discovery of new species is highly expected given their rate of recent discoveries and finding of hidden cryptic diversity, as noted previously (e.g., Benda et al. 2012 Urva due to lack of information and possible misidentification of species in each genus.

Biogeography: biodiversity patterns in a complex area
To classify bioregions, the conventional distance-based hierarchical clustering method and a network-based method built on a community-detection algorithm known as Infomap were used by incorporating distributional data from 188 native species (Yusefi et al. 2019c). The latter method was applied using the interactive web application "Infomap Bioregions" and the former was based on Simpson's dissimilarity index (βsim) and Unweighted Pair-Group Method, using Arithmetic averages (UPGMA) (Kreft and Jetz 2010). Infomap Bioregion adaptively organizes species records into discrete geographical grid cells such that the data density determines the spatial resolution. Next, it generates a bipartite network between species and grid cells, and the bipartite network were then clustered by applying the Infomap clustering algorithm (see more details in Edler et al. 2017).
The network-based simulations resulted into the final list of seven bioregions and two transition zones as the most likely number of clusters (Fig. 3). These bioregions were identified based on specimen records described in the previous section and default clustering parameters (maximum cell size = 1°, cell capacity = 10-300, for a total of 160 grid cells). Of these, two bioregions largely matched the Alborz-Zagros mountains and Central Basin lowlands. The two transition zones were located in boundary areas, one where the Iranian Plateau and Arabian Peninsula meet and the other between the two largest identified bioregions (Fig. 3). Infomap also identified the most indicative and common species of each bioregion. The distance-based method suggested five clusters as the most likely number of bioregions, including one large bioregion that covers the entire Central Basin and Alborz-Zagros mountains, but not the Caucasian region (Fig. 3).
This study documented the efficacy of Infomap for identifying bioregions at regional scales, and it was apparently more sensitive (by detecting nine biogeographical units) than the clustering approach which suggested five. The detection of transition zones by the Infomap, in particular, was an advantage (as found at larger scales by Bloomfield et al. 2018). The results also stressed the fact that the biogeography of Iran is highly complex, particularly in the south-eastern part where we can find contact areas between several bioregions.

Extinction risk: threats versus traits
By using an integrative model incorporating multiple predictor variables as well as phylogenetic information, I examined the relationships between response variable (IUCN Regional Red List status) and 11 predictors of extinction risk (Supplementary Table S1), including direct (hunting/persecution) and indirect (human influence index, HII) human impacts as well as intrinsic biological/morphological traits (e.g., fecundity, growth, size) and extrinsic ecological/environmental factors (e.g., diet, habitat breadth, primary productivity), using linear regressions in R (R Development Core Team 2018; G.H. Yusefi, J.C. Brito, and K. Safi unpublished). As a source for mammal phylogeny, I used an updated species-level supertree (Faurby and Svening 2015) to derive the pairwise (co-phenetic) phylogenetic distances using 'ape' in R (Paradis and Schliep 2018). The ecological, environmental, lifehistory and morphological variables were collated from the database PanTHERIA (Jones et al. 2009), using only those variables for which data were available for >50% of the cases in order to reduce bias that may arise from the missing data (González-Suárez et al. 2012). The species-specific mean HII was calculated by overlaying the HII of Iran (which was extracted from the global dataset; WCS, CIESIN 2005) with the extent of occurrence of each species (Table S2)  up to high penalties for violations involving the killing of special concern or charismatic species (for details see Supplementary Tables S1, S2). After accounting for autocorrelation (using 'car' R package;Fox et al. 2019), relationships between the response and 11 predictor variables were evaluated by linear regressions, and a stepwise heuristic procedure was used to identify the significant predictors and find the best-fitting models based on the Akaike Information Criterion (AIC). The adjusted R-squared (R 2 ) was then used to quantify the amount of explained variation by the fitted models (following Safi and Pettorelli 2010). I first analysed all species together with adult body mass as a continuous explanatory variable. I then ran analyses separately for small, and medium/ large species (with a cut-off of 1 kg) due to the pronounced differences in factors explaining extinction risk in these groups. Small mammals (n = 117) included all species belonging to orders Eulipotyphla, Chiroptera,  Table S2).

Rodentia (except
The evidence of phylogenetic non-independence in the model residuals was tested using generalized linear mixed models (GLMM), with glmmPQL function in 'mass' and "Moran's I" in 'ape' R packages (Paradis and Schliep 2018). Lastly, I estimated goodness-of-fit (using Gtest) and did a posterior heterogeneity analysis of the taxonomic distribution of the risk of decline (i.e., declining species or threatened ones: those categorised as Vulnerable, Endangered, or Critically Endangered) to see if declining risk is randomly distributed among orders (if so, the expected number of declining species in each order should be proportional to the number of declining species in the total species assemblage).

Comparisons of AIC values showed that the best performing model included Hunting
Vulnerability measured via DoE fines (HV1), Gestation Length (GL), Weaning Age (WA), and Actual Evapo-Transpiration (AET), and the fit of this model slightly decreased when excluding AET and WA (ΔAIC = 1.98). The best models explained a large amount (71-73%) of the variation in extinction risk (Table 1). In every model, extinction risk showed a strong positive association with HV1 (p ≤ 0.001). Predictors of extinction risk in small (< 1 kg) and medium/large (> 1 kg) mammals exhibited different patterns. In small mammals HII was the only significant predictor of extinction, but in medium/large mammals the best model included HV1, GL, WA, and AET, explaining a large amount of the variation (68-81%) in which the variable most strongly associated with species vulnerability was HV1 (p ≤ 0.01; Table 1).
No phylogenetic signal was detected for any set of variables tested (for example, in HV1; r = 1.37, SE = 0.12, df = 154, p = 0.000), and the Moran I coefficients also revealed no phylogenetic autocorrelation (Moran's I test all: p > 0.05). The observed distribution of threatened species across orders differed significantly with that expected by chance and were highly heterogeneous (G = 18.72, df = 6, p < 0.001). Carnivores (order Carnivora) and Ungulates (Artiodactyla and Perissodactyla orders) exhibited more declining species than expected, while the opposite trend was observed in other orders.
These results explicitly stressed the strong influence that human direct impacts (i.e., hunting/persecuting), rather than indirect impacts of human activities or species' intrinsic and extrinsic traits, have on the vulnerability of these taxa threatened to extinction, thus emphasizing the need to protect these species which are more impacted by humans, especially medium and large ones. The results further support the conclusion that the threats, which are currently in early stages of development in comparative extinction risk analyses (as noted by Murray et al. 2014), are important and should be seriously considered in such studies. Similar to global findings (Ripple et al. 2016), this study shows that hunting is an important predictor of decline for mammals at a regional scale. Moreover, understanding how human threat differentially affects larger and smaller mammal species may have important implications for conservation actions.
Climate change: functional diversity vulnerabilities in a high-risk area 16 Following previous studies (Newbold et al. 2013,2020, García-Llamas et al. 2019) and based on the species' life-history and ecological traits (body size, diet, climate breadth, habitat use, and locomotion types), the species assemblage was classified into 12 functional groups: small, large, carnivore, herbivore, insectivore, climate generalist, climate specialist, desert species, forest species, mountain species, non-volant and volant. Two additional groups of threatened and total species richness were also considered. The spatial distribution of these groups was determined, and the geographic location of their hotspots were mapped. Then, the level of exposure of each group's hotspots were quantified by intersecting their range with climate risk areas to evaluate how the predicted magnitude of climate change in the region might affect them. Climate risk areas or areas with high risk of exposure to extreme climate change were mapped using current and future (2070) mean annual precipitation (mm) and mean annual temperature (°C) data 3 . These risk areas were defined as areas above the upper critical temperature (UCT) threshold (after calculating average and standard deviation of UCT from the mammal species occurring in Iran based on the largest sets of endotherms UCT data available; Araújo et al. 2013, Bennett et al. 2018) and below dryland threshold (200 mm/year precipitation; Ward 2016). The future climate was projected using 19 different global climate models (GCMs) and two representative concentration pathway (RCP) scenarios (RCP4.5 and RCP 8.5). All geographic and statistical analyses were performed in ArcGIS ver. 10.4.1. (ESRI, 2016). The velocity of climate change (following Loarie et al. 2009) was calculated for the same climatic variables (precipitation and temperature) covering the extent of the study area as an average value of the available combinations of 19 GCMs, two RCPs (4.5 and 8.5) and two time periods (current time and 2070). The velocity of climate change was defined as the ratio of the temporal gradient of change to the spatial gradient (following Loarie et al. 2009). The histograms of the velocity of climate change for each functional group area of occupancy were plotted and compared with the mean velocity of climate change of each group. The workflow was implemented in R using 'raster' package (Hijman et al. 2018).
The results showed that vulnerability to the velocity and magnitude of climate change varied between the functional groups and the threatened and total species (G.H. Yusefi, K. Safi, P. Tarroso and J.C. Brito unpublished). Functional groups were more exposed to precipitation decrease than to temperature rise. In the present time no single functional groups' hotspot is exposed to temperatures above the UCT threshold, but in 2070, eight groups (except insectivore, small, specialist and volant), as well as threatened and total species, will likely be exposed to extreme temperatures, with the highest values found in the desert group (16.7% of hotspot areas; Fig. 4, Supplementary Fig. S2). Currently, all functional groups (except insectivore and volant), as well as threatened and total species, are already exposed to annual precipitation under 200 mm, but in 2070, the number of exposed cells is predicted to increase and affect all groups, with the highest values found in the carnivore, desert, and large-sized groups (43.8%, 53.6%, 49.1% of hotspot areas, respectively), along with the threatened group (47.7%; Fig. 4, Supplementary Fig. S2). These findings remained relatively unchanged when we only consider functional group ranges covered by existing protected areas given the large portions (above 50%) of currently protected desert hotspots; carnivore and large-sized functional groups and threatened species will be exposed to lower precipitation levels. Climate change velocities estimated for 12 functional groups, threatened and total species, were invariably positive for all cases, while velocities were particularly high for carnivore, large-sized, specialist, and threatened groups in both RCP 4.5 and RCP 8.5 scenarios (Supplementary Fig. S3; the results for RCP4.5 were similar and are not shown). Functional groups were more exposed to areas with higher velocities of change in temperature according to their range in comparison to precipitation velocity.
The results show that different assemblages of species will be affected differently by the magnitude and velocities of climate change. It also reveals that a large proportion of the current range of carnivore, desert, large-bodied functional groups and threatened species is currently jeopardized by the projected climate change, and it would be even more at risk in the next few decades. Contrary to previous studies, I found that lowland (rather than mountain areas; e.g., Maiorano et al. 2013) and large species (rather than small ones; e.g., Vale and Brito 2015) are more vulnerable to future climate change. Lastly, drought, rather than warming, is the factor that most likely will affect regional mammal biodiversity in the arid region of southwest Asia.

Conclusions
This study has advanced the understanding of mammalian biogeography in a complex region. The overall patterns found here underline the high vulnerability of large and threatened species to extinction and climate change while emphasizing the need for more research, in particular, to understand the various ways that direct killing and changing climate affect these declining species. There is therefore an urgent need for increasing conservation efforts, especially in the case of large and threatened species. Taken together, these findings provide a foundation for future biogeographic, systematics and ecological studies of Iranian mammals while simultaneously adding to the limited existing information on the bioregionalization process at regional scales, as well as highlighting the importance of incorporating threats and functional traits information in predicting extinction risk and in addressing the impacts of climate change, respectively.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb   Table S1. Description of variables used in the extinction risk analyses.  Submitted: 9 September 2020 First decision: 1 October 2020 Accepted: 11 February 2020 Edited by Janet Franklin