Genetic structure and ecological niche segregation of Indian gray mongoose (Urva edwardsii) in Iran

Abstract Combining genetic data with ecological niche models is an effective approach for exploring climatic and nonclimatic environmental variables affecting spatial patterns of intraspecific genetic variation. Here, we adopted this combined approach to evaluate genetic structure and ecological niche of the Indian gray mongoose (Urva edwardsii) in Iran, as the most western part of the species range. Using mtDNA, we confirmed the presence of two highly differentiated clades. Then, we incorporated ensemble of small models (ESMs) using climatic and nonclimatic variables with genetic data to assess whether genetic differentiation among clades was coupled with their ecological niche. Climate niche divergence was also examined based on a principal component analysis on climatic factors only. The relative habitat suitability values predicted by the ESMs for both clades revealed their niche separation. Between‐clade climate only niche comparison revealed that climate space occupied by clades is similar to some extent, but the niches that they utilize differ between the distribution ranges of clades. We found that in the absence of evidence for recent genetic exchanges, distribution models suggest the species occurs in different niches and that there are apparent areas of disconnection across the species range. The estimated divergence time between the two Iranian clades (4.9 Mya) coincides with the uplifting of the Zagros Mountains during the Early Pliocene. The Zagros mountain‐building event seems to have prevented the distribution of U. edwardsii populations between the western and eastern parts of the mountains as a result of vicariance events. Our findings indicated that the two U. edwardsii genetic clades in Iran can be considered as two conservation units and can be utilized to develop habitat‐specific and climate change‐integrated management strategies.


| INTRODUC TI ON
A number of biological and ecological processes are known to shape patterns of genetic structure in continuous or discrete populations of species such as geographic, ecological, or reproductive barriers (e.g., Barton, 2008;Bradburd et al., 2013). Landscape models have been developed to examine the roles of landscape, that is, isolation by resistance (McRae & Beier, 2007), and environmental niche, that is, isolation by environment (Wang & Summers, 2010) on gene flow.
Despite the importance of environmental niche dissimilarity as a motivator of genetic differentiations, its dynamism and strength are poorly known.
Spatially explicit environmental data in the framework of ecological niche models (ENMs) have been widely used to generate geographic predictions of a species' distribution (Guisan & Thuiller, 2005). However, in ENMs, a genetic uniformity throughout a species' range is assumed and the potential for local adaptation to specific environmental conditions is ignored (Gotelli & Stanton-Geddes, 2015;Grady et al., 2011Grady et al., , 2013. Research has been shown that incorporating information on genetic variability and local adaptation improves model performance (e.g., Ikeda et al., 2017;Marcer et al., 2016) and helps to identify processes that have shaped population structure of species in the past and present (e.g., Alvarado-Serrano & Knowles, 2014;Chang et al., 2012;Knowles et al., 2007;Marcer et al., 2016;May et al., 2011).
Ecological niche differences among species or populations can be analyzed to evaluate the possible ecological and evolutionary forces that shape geographic distributions, habitat preferences, and genetic structures (e.g., Raxworthy et al., 2007). A combination of ecological niche models and genetic analysis has potentially broad applications for defining conservation units, including evolutionary significant units (ESUs; e.g., Crandall et al., 2000), management units (MUs), and distinct population segments (DPSs; The Endangered Species Act, 1978) for a broad range of taxa.
Here, we employed an integrative approach to explore genetic and ecological differentiations within the Indian gray mongoose (Urva edwardsii, Geoffroy Saint-Hilaire 1818) distributed along the southern part of Iran. Southern Iran consists of the southern mountain ranges of Zagros, Central Iranian Range, Khuzestan Plain, and the northern coasts of Persian Gulf.
The southern Zagros mountainous ridge has been regarded as one of the most known natural barriers between central Iran and the Mesopotamian plain to the west, responsible for the vicariance event in the region for different taxa (e.g., Ghaedi et al., 2020;Nazarizadeh et al., 2016). The presence of unsuitable environmental conditions or geographic distance may limit gene flow among populations, resulting in different genetic clades; thus, the present distributional patterns of species may be due to vicariant events (Nazarizadeh et al., 2016;Nilson et al., 2003). Climate also act as the primary driver of adaptive divergence and shaping genetic evolution, by limiting individuals' establishment and dispersal (Pearson & Dawson, 2003). However, there is no information as to whether the Zagros Mountains have acted as a barrier for small carnivores such as mongoose species.
Smaller carnivores such as mongooses lack the attraction and attention that managers and conservationists seek in mega-carnivores as suitable flagship species and conservation tools. Thus scientific information on small carnivores such as the Indian gray mongoose is still scarce in many countries including Iran. Small carnivores, however, show stronger specialization and resource selectivity compared to large carnivores (Kalle et al., 2012), and then, they may serve as useful indicator species in the preservation of keystone habitats. Knowledge on spatial scale and landscape heterogeneity is integral to the understanding of species habitat associations. Several factors including trapping for meat or fur (used in shaving and paint brushes), and good luck charms threaten the viability of mongooses (Mallick, 2009). In addition, limited data are available on the taxonomic status, habitat requirements, ecological niche differentiations, and genetic distances among the populations of U. edwardsii across its distribution range, which hinders the effective management of the species. Therefore, the current study aimed to (i) provide insights into intraspecific genetic diversity and niche differentiations among U. edwardsii populations in Iran and (ii) assess whether the Zagros Mountain may have acted as a physical barrier between the Indian gray mongoose populations.

| Study area
The study area covers the distribution range of Indian gray mongoose in Southern Iran (Ziaie, 2008) (Figure 1). The species inhabits a variety of habitats including meadows, sparse woodlands, shrublands, croplands, habitats near to reedy wetland, and palm trees near villages (Karami et al., 2016). This region, which is mainly located in zoogeographic realm of Saharo-Arabian (Holt et al., 2013), is affected by subequatorial climate condition with cold and rainy winter and warm and dry summer (Heshmati, 2012). The annual precipitation ranges from 50 to 650 mm with an increasing gradient from east to west. Mean annual temperature ranges from 2 to 28°C, where January and February are the coldest, and July and August are the warmest months.
Urva edwardsii has recently been introduced into some parts of Central Iran such as Isfahan Province (Yusefi et al., 2019). While U. edwardsii is found from Central to Southern Iran, the small Indian mongoose (Urva auropunctata, Hodgson 1836) only reported from Sistan and Baluchestan Province (Sistan region), Hormozgan Province, and few localities in Khuzestan Province (Yusefi et al., 2019 and references therein). The species was also introduced to many islands and mainland to acts as a biological control of rats and snakes in plantations (Barun et al., 2013;Simberloff et al., 2000;Thulin et al., 2006), and it is considered as one of the worst invasive alien species (Lowe et al., 2000).
2.1.1 | Sampling, DNA extraction, and sequencing To assess historical relationships among U. edwardsii populations, 24 specimens were collected across the species range in Iran ( Figure 1).
In addition, we collected three samples of U. auropunctata from southwest of Iran. All samples were obtained from road kills or natural cause deaths between 2016 and 2018 and stored in 98% ethanol.
Also, 17 control region sequences of Urva javanica, U. edwardsii, and U. auropunctata were obtained from GenBank. The locality information, ID, and accession numbers are summarized in Table 1.

| Phylogenetic analysis
Iranian sequences of U. edwardsii (24 sequences) and U. auropunctata (3 sequences), and sequences from GenBank, including eight sequences of U. edwardsii, six sequences of U. auropunctata as well as three sequences of Urva javanica, were aligned using the Clustal W algorithm (Thompson et al., 1994) implemented in Mega v.7 (Kumar et al., 2016), and final adjustments were made by eye. In addition, the sequence of Herpestes brachyurus (KY117547) was used as the outgroup. The HKY+ G was selected as the best-fitting model of DNA substitution, using Bayesian information criterion (BIC) implemented in jModelTest (Posada, 2008). Alignments of noncoding regions of mtDNA frequently present gaps, and thus, gaps were treated as fifth base in the marker.
Bayesian phylogenetic analyses were carried out in MrBayes v.3.2 (Ronquist et al., 2012) with two independent runs of four Markov chains (one cold and three heated) over 10,000,000 generations and sampling every 1,000 generations. The first 25% of the sampled trees and estimated parameters were discarded as burn-in. Convergence of the model parameters was monitored using the program Tracer v.1.7.1 (Rambaut, Drummond, Xie, Baele, & Suchard, 2018 (Excoffier & Heckel, 2006), and a significance test of 10,000 permutations, of the obtained estimates of F ST , was performed. In addition, three methods for species delimitation were used to identify the specific boundaries in the species: (i) the Poisson tree process model (PTP; Zhang et al., 2013), The divergence time of the major clades was estimated using a Bayesian molecular clock approach implemented in the program BEAST v2.4.7 (Bouckaert et al., 2014). All molecular clocks need to be calibrated using independent information such as fossils. Fossils attributed to Urva are scarce (Peigné et al., 2005). The oldest specimen F I G U R E 1 Genetic samples and species presence localities of Indian gray mongoose in south-southeast (Clade 1, blue) and southwest of Iran (Clade 2, red) Markov chains were run for 3 x10 7 generations, sampling every 1,000 generations was employed to construct a maximum clade credibility tree. The log parameters and tree files from each run were combined, using the LOGCOMBINER v.1.8. (Rambaut & Drummond, 2015), and convergence diagnostics were assessed using TRACER v.1.6 (Rambaut et al., 2007).
Measures of DNA polymorphism were estimated using DnaSP v.5 (Librado & Rozas, 2009), including the number of haplotypes (H), haplotype diversity (h), and nucleotide diversity (π) for each clade separately and the whole samples in Iran. A median-joining (MJ) network was constructed for U. edwardsii using PopART v.1.7 (Leigh & Bryant, 2015) with the default settings.

| Ecological niche modeling
Species presence localities were compiled opportunistically from different sources including genetic specimens, camera-trap detections, trappings, and environmental guards' sightings across the species range (Figure 1). The presence localities of Indian gray mongoose were explored on biodiversity datasets (i.e., VertNet, iDigBio, and Arctos). However, the presence points on these datasets was not included in the ESM analyses for two reasons: (1) There was a significant overlap among presence points on these datasets and those that we compiled for the analysis, considering 1-km buffer around each point, and (2) the accuracy of the presence localities recorded in these databases was unknown, especially for countries with scarcity of data or studies that are conducted at a small scale.
We excluded the Indian gray mongoose occurrence localities in the contact zone of the two identified clades in the process of ecological niche modeling of clades as we did not have genetic samples from their regions, and therefore, their genetic clades were unknown. However, these localities were included in the ENM using all occupied localities clades together. Spatial autocorrelation in presence localities was assessed by Moran's I index for each variable separately in R package raster (Hijmans & van Etten, 2011 To create a reliable map of highly suitable habitats for each clade considering the low number of occurrence localities, we adopted ensembles of small models' (ESMs) (Breiner et al., 2015(Breiner et al., , 2018 in R packages ecospat (Broennimann et al., 2015) and biomod (Thuiller et al., 2021). ESMs represent a novel technique for predicting spatial niche of species with few observations (Breiner et al., 2018). By reducing the number of predictor variables and averaging simple small bivariate models to an ensemble, this approach avoids overfitting. For each clade, numerous bivariate maximum entropy models (MaxEnt) (Phillips et al., 2006) were calibrated and evaluated. Then, the final ESMs for each clade were calculated using the weighted average of the all resulting Somers' D (i.e., rescaled AUC) values of the bivariate models.
Models with a Somers' D lower than 0 were not included in the ESMs.
All bivariate models were calibrated with the threefold cross-validation and 10,000 background points and were calibrated using the 70% of occurrence localities as training and 30% as evaluation data.
We calculated AUC and continuous Boyce index (Hirzel et al., 2006) to evaluate the performance of each model separately and overall ESM. All distribution analyses were conducted at a spa-  The haplotype network provided further information regarding the haplotypes (Figure 4). Haplotype 1 was the most widespread haplotype of U. edwardsii and was found within 10 individuals in Clade 1 in four provinces of Fars, Kerman, Yazd, and Bushehr.

| Phylogenetic analysis and genetic structure
Haplotype network showed divergence between the Indian gray mongooses in Iran. In line with the results obtained from the phylogeny tree (Figure 2), the haplotypes derived from the U. edwardsii were divided into two separate groups (i.e., clades 1 and 2).
Only one unique haplotype was detected among three sequences of U. auropunctata, but nine unique haplotypes were detected among the sequences generated for CR of U. edwardsii with an overall haplotype and the nucleotide diversity h = 0.804 and π = 0.00627, respectively (  (Table 3).
The ESMs showed a notable consistency in predicting habitat suitability of both clades when compared with occurrence records.
The potential distribution of U. edwardsii for Clade 1 is almost two times larger than the genetic Clade 2 ( Figure 5). The relative habitat suitability values predicted by the ESM for the U. edwardsii in both clades revealed their niche separation ( Figure 5) of the variation in climatic variables, respectively. PC1 was mostly correlated with temperature-related factors (e.g., mean temperature of wettest quarter, temperature annual range, and mean temperature of warmest quarter), while PC2 was highly correlated with precipitation-related factors (precipitation of coldest quarter, and precipitation of wettest month; Table 4). comparing in the opposite direction, the observed overlap is higher (a more similar ecological niche), yet not significant. Randomization tests (i.e., niche identity and background similarity) showed that environmental space occupied by clades is similar to some extent, but the niches that they utilize (identity test) differ between the distribution ranges of clades.

| D ISCUSS I ON
The present study was the first to employ a combination of genetic data and ecological niche models to evaluate intraspecific genetic and ecological variation within a less studied carnivore, the Indian gray mongoose. Our findings provide evidence in supporting the idea that incorporating ecological niche models into phylogeographic analyses would improve the accuracy of intraspecific genetic information and predictive ability to population boundaries. The adopted methodology has great implications for conservation unit assessment for the species and can provide valuable information regarding the feasibility of this approach to many species around the world. We showed two main results, which are elaborated in the sections below. First, phylogenetic analysis suggested the existence of two genetic clades within the U. edwardsii populations in Iran, one in southwest and the other in south and southeast. Second, a genetic clades-specific combination of environmental variables among the clades was found, indicating the niches that each U. edwardsii clade utilizes differ between the distribution ranges of clades. Bangladesh (UeTC294-UeTC296), UAE (UeTC146), and Netherland (UeTC293) to another. Also, one specimen from Nepal (UeA163188) is placed as a sister group to the two Middle Eastern clades (see Veron and Jennings (2017) for accession numbers). Here, available CR sequences (Veron & Jennings, 2017), that is, UeTC294-UeTC296, UeTC144-UeTC146, UeTC148, and UeFj678485 were aligned with sequences of current study to reconstruct the phylogenetic tree.
The CR phylogenetic tree and species delimitation approaches have resulted in phenotypic divergences, genetic difference, and accelerating the evolution of reproductive isolation (Kozak et al., 2008;May et al., 2011;Ruiz-Sanchez & Sosa, 2010). Both biotic and abiotic factors are involved in determining spatial patterns of the species' genetic variation and promoting adaptive divergence and speciation (Peterson et al., 1999;Rissler & Apodaca, 2007 (Salah, 2019;Sborshchikov et al., 1981). The separation of the Arabian plate from Africa was accelerated in the beginning of the Pliocene (ca. 5 Mya) (Girdler, 1984) and led to mountain-building on the western margins of the Iranian Plateau (Macey et al., 1998). There are potentially genetic differences between the two regions based on mtDNA, and that the two regions have different environmental characteristics, which may lead to different adaptations.
Although ecological niche models are often used to map distribution range of threatened species and evaluate environmental factors affecting them (e.g., Ahmadi et al., 2017;Khosravi et al., 2019;Shahnaseri et al., 2019), hypothesis regarding ecological exchangeability (May et al., 2011), niche segregation between species, subspecies, or populations of same species in different regions or seasons, and intraspecific or interspecific taxonomic distinction (Nakazato et al., 2010) can be tested using these noninvasive and multifunctional tools. We performed niche differentiation tests to identify the environmental constraints for the distribution of U. edwardsii and explore the hypothesis that adaptation to different environmental condition are contributing to genetic structuring in Indian gray mongoose. In accordance to genetic findings, our ecological niche models showed little geographic overlap between two genetic clades of U. edwardsii, indicating that clades differ in their niches and occupy a part of the geographic range in south of Iran where environmental conditions are suitable ( Figure 5). Also, our environmental niche model showed that ensemble of small models of each genetic clade alone (Figure 5a and b) was better than models of two genetic clades together ( Figure 5c). However, the distribution models based on the all presence localities cover areas highly suitable for both genetic clades (Figure 5c), and the ESMs of the two clades separately showed that Clade 1 is commonly found at localities that are less suitable for Clade 2 and vice versa.
Although the distribution of both clade appears to be highly dependent on the herbaceous with sparse density of tree and shrub (Table 3), climatic and nonclimatic environmental variables contributed differently to each clade climate niche: "precipitation seasonality" and "human footprint" for the Indian gray mongooses in south-southeast and "mean temperature of warmest quarter" and "precipitation of coldest quarter" for the Indian gray mongooses in southwest of Iran. It is plausible that human footprint is import-

| CON CLUS ION
In conclusion, our results showed that two genetic clades of Indian gray mongoose occupy distinct environmental spaces that are similar, but not equivalent, relatively small niche overlap between these clades. These findings show that two clades of U. edwardsii underwent significant changes of their range in environmental niche space during the process of occupation of region within their present distribution. The estimated divergence time between the two Iranian clades (4.9 Mya) coincides with the uplifting of the Zagros Mountains during the Early Pliocene. The Zagros mountain-building event seems to have prevented the distribution of U. edwardsii populations between the western and eastern parts of the mountains as a result of vicariance events.

| Conservation implications
Our predicted ecological niche maps broadly estimate the most important ecological requirements that the species may need to establish a habitat. Also, the obtaining distribution data are informative to establish conservation priority and conservation efforts. As our genetic data confirmed the genetic structure among the populations of U. edwardsii, thus, a conservation status assessment should be conducted for both clades separately and may result in a subsequent definition of cryptic subspecies and conservation efforts. Also, it is important that conservation actions be taken in the immediate future to conserve the genetic variation of population in each region before it declines. Despite that the impact on native fauna of introduced Indian gray mongooses is not known, further introductions of any mongoose species to act as a biological control outside their native range are not recommended.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no competing interests.