Morphological diversity and ecological niche divergence in goitered and sand gazelles

Abstract The phylogeny and species boundaries of Gazella subgutturosa and G. marica have been long debated. The achievements of past conservation efforts have been compromised by a lack of knowledge about the phylogeny and taxonomic status of different populations. We integrated the recent genetic findings by previous studies with morphometric analyses and ecological niche modeling (ENM) to assess discreteness among populations of these gazelle species in Asia. Taxonomic diversity of gazelles was investigated by using principal components analysis (PCA) based on 14 cranial measures of male skulls. Ecological niche divergence was examined based on a PCA on climatic factors and a species distribution modeling (SDM) with environmental variables. Morphometric results indicated substantial differentiation in size between skulls of the western Zagros Mountains including west and south‐western Iran and Arabian Peninsula from all other samples east of the Zagros Mountains from Iran to China. ENM also revealed that gazelles in the east and west of Zagros Mountains occupy distinct niches and that there are apparent areas of disconnection across the goitered gazelle suitable range. A complete divergent niche occupation was also observed between goitered gazelles of northern Mongolia and other populations of the species, except those in China. Taking the inferences from ENM and morphology together with previous genetics results, we conclude that gazelles in the west and south‐west of Iran may represent G. marica. Also, our combined analyses revealed divergence among gazelles of Iran, Central Asia, and Mongolia/China. These results may pave the way for future studies and have conservation implications particularly for reintroduction/supplementation programs.


| INTRODUC TI ON
Gazelles (genus Gazella) are widely distributed through mainly desert country from North Africa via the Middle East to Mongolia and India, although in most places, they have been depleted by overhunting. The G. subgutturosa is the most widespread species group in the genus Gazella (Groves & Grubb, 2011), ranging from Oman across the Arabian Peninsula to southern Turkey, and following the steppes of Central Asia eastwards into Central Mongolia, China, and Pakistan (Groves & Grubb, 2011). This wide geographic distribution encompassing diverse biogeographic zones has formed divergent populations in the species.
The taxonomy of sand gazelle (G. marica) and goitered gazelle (G. Subgutturosa) has been the subject of much debate. Sand gazelle was long considered to be a subspecies of G. subgutturosa on the basis of morphology and karyology (Groves & Harrison, 1967;Kingswood, Rebholz, Vassart, & Kumamoto, 1997;Kingswood, Vassart, & Williamson, 1996). Thomas (1897) described the Arabian sand gazelle as a distinct species, G. marica (Thomas, 1897); Ellerman and Morrison-Scott (1951) considered it as a subspecies of the Saharan G. leptoceros (Cuvier, 1842); and it was later suggested by Groves and Harrison (1967) that it is associated more closely with G. subgutturosa, probably as a subspecies of it. Wacher et al. (2011) provided a phylogenetic framework based on the analysis of mtDNA sequences of a number of wild and captive individuals throughout the species' natural range and restored G. marica to full species which is distributed in Arabian Peninsula, Iraq, Jordan, Syria, and Turkey. Other authors (Durmuş, 2012;Groves & Harrison, 1967;Kingswood & Kumamoto, 1988; indicated that gazelles in a broad geographic area between Euphrates-Tigris basin and the Zagros Mountains of Iran have morphologically mix characters. More information on the taxonomy of goitered gazelle is provided by Lerp et al. (2016), Abduriyim, Zibibulla, Eli, Ismayil, and Halik (2018), Silva (2019), and Fadakar et al. (2020).
There is a paucity of data on the phylogenetic status of gazelle populations in the northern geographic range of G. marica . Fadakar et al. (2019) found mitochondrial haplotypes of G. marica in south-western Iran. They attributed their finding to either hybridization between G. subgutturosa and G. marica in their contact zone, or the existence of G. marica populations in Iran. Accordingly, the dunes and desert areas of south-west and west of Iran at the border of Iran and Iraq are probably the historical habitat of sand gazelles. The resemblance of gazelles in south-western Iran to G. marica has hitherto been noted by Groves and Harrison (1967). They described the size and morphology of gazelles from south-western Iran as an intermediation between G. subgutturosa and G. marica. Groves (1997) suggested that gazelles in eastern Iraq and south-western Iran have probably been arisen through secondary intergradation between G. s. marica and G. s. subgutturosa. Groves and Grubb (2011) recognized four phylogenetic species for G. subgutturosa group including G. subgutturosa, G. gracilicornis, G. yarkandensis, and G. marica. Subsequent studies suggested the elevation of the latter to full species (Wacher et al., 2011), but the three others are generally considered as subspecies of G. subgutturosa (Abduriyim et al., 2018;Sorokin, Soldatova, Lukarevskiy, & Kholodova, 2011). However, based on mitochondrial genes, Fadakar et al. (2020) suggested that G. s. gracilicornis is synonymous with G. s. yarkandensis. Fadakar et al. (2020) recognized Lut and Kavir deserts in Central Iran as the main geographic barriers between western (Central Iranian; G. s. subgutturosa) and eastern (Asiatic; G. s. yarkandensis) populations of goitered gazelle and suggested that the distribution of G. s. subgutturosa extends to the west of the Zagros Mountains.
Another case of taxonomic debate about gazelles in south-western Iran is the status of G. karamii, a taxon initially described by Groves (1969) as a subspecies of G. bennettii (Sykes, 1831) from Borazjan, south-western Iran near the Persian Gulf. The specimen was later transferred to the G. gazella group, mainly due to its dark pelage (Groves & Grubb, 2011). However, based on morphometric analysis, Bärmann et al. (2013) found that the type of G. karamii was close to G. marica, and they suggested that it was possibly a synonym of it. Until now, no molecular data exist for the only specimen of G. karamii preserved in the MfN Berlin museum (ZMB MAM 41400).
There has also been taxonomic uncertainty in gazelles of Khark (or Kharg) Island in the south of Iran. Based on cranial and morphometric analyses, Karami, Hemami, and Groves (2002) suggested that Khark Island population along with those in the west of the Zagros Mountains in Bushehr and Khuzistan provinces should be considered as a distinct subspecies of G. subgutturosa. Later, Mirzakhah, Naderi, Rezaei, Fadakar, and Naseri (2015) indicated that there is no evidence of G. marica haplotypes in this population.
Effective biodiversity conservation requires knowledge of ecological and phenotypic variation among taxa and their evolutionary relationships (Noguerales, Cordero, & Ortego, 2018). For many gazelle taxa, the achievements of past conservation efforts have been compromised by a lack of knowledge and confusion about the phylogeny and taxonomic status of different populations, and even species boundaries have not been certain (Groves & Grubb, 2011).
Combining multiple sources of data, including genetic, ecological, and phenotypic data, can provide greater insight into the species boundaries (Ahmadi et al., 2018;Leaché et al., 2009). In this regard Lerp et al. (2016), Silva et al. (2017), andFadakar et al. (2020) suggested the use of morphological and nuclear genetic studies to investigate further the relatedness of gazelle populations on the two sides of the Zagros Mountains.
We used phenotypic (cranial) and ecological (niche) data to examine whether gazelle populations on the two sides of the Zagros Mountains exhibit distinct phenotypic and ecological traits. Our study also includes skull measurements and occurrence locations collected across the entire range of G. subgutturosa, and accordingly, can assist in better understanding of the continent-wide intraspecific variations among the species populations.

| Species concept
Since the 1990s, there has been much discussion about "what is a species"; de Queiroz (2007) cogently argued that a species is an evolutionary lineage, and that traits such as reproductive isolation, specific mate recognition systems, an ecological niche, increased genetic differentiation, and subjective "enough difference," which have often in the past been used actually to define a species, develop along the lineage over time; the initial marker of the individuated lineage that is a species is likely to be diagnosability. In this paper, we follow de Queiroz's unified concept of species defined as "separately evolving metapopulation-level lineages" (de Queiroz, 2007).
Based on morphometric analyses and as suggested by Wacher et al. (2011), we considered G. marica a distinct species from G. subgutturosa. We followed Fadakar et al. (2020) for the intraspecific classification of G. subgutturosa.  Figure 1). Only adult male skulls were included in the study as not enough female skulls were available.

| Data collection
Adult specimens were recognized by the complete eruption and full occlusion of premolars and molars (Angelici & Luiselli, 1999).
Specimens from across the range of the group were studied between 1980 and 2010 in several museums and institutions. The measurements were taken by C. P. Groves and M. R. Hemami. Specimens kept in the Iranian National Museum of Natural History (MMTT) were measured by both authors; the difference in measurements was <0.2 mm. The list of skull specimens and their numbers in PCA plots is presented in Table S1 and Figure S1 in Online Resource 1.
For cranial evaluation, fourteen linear skull parameters were chosen (Table S2; Figure S2). Mandibular measurements were excluded because the mandible was missing in many skulls. All measurements were taken with vernier callipers to the nearest 0.1 mm.

| Data transformation
Because multivariate analyses are sensitive to missing data (Daneri, Esponda, De Santis, & Pla, 2005), specimens were excluded if data were missing for more than one measurement. To estimate the missing value for each specimen with only one missing data point, we used stepwise regression from the remaining subset of cranial variables available for that specimen (Reig, 1992;Reig & Ruprecht, 1989).
In order to separate the effects of size and shape on morphology, we followed the approach of Jungers, Falsetti, and Wall (1995) and Klein, Franciscus, and Steele (2010). We used three derived variables from the original raw measurements to investigate the craniometric differences among specimens: (a) log 10 of each measurement on each specimen was taken to allow for normal distribution and homogeneity of variances (log size and shape variables; Lewontin, 1966).

F I G U R E 1
Locations where gazelle skulls (large colored circles) and presence points of Gazella subgutturosa and G. marica (small colored circles) were collected. See Table S1 and Figure S1 for geographic names and locality abbreviations, and gazelle skulls number respectively To obtain an index for overall size (Darroch & Mosimann, 1985;Mosimann & James, 1979), the geometric mean of all cranial dimensions for each individual was calculated (Klein et al., 2010). (c) The geometric mean of the specimen was subtracted from logged measurements. This difference is a measure of log-shape and ratios between each raw measurement and overall size or geometric mean.
In log-shape variables, the variance is partitioned due to the specimen's size and shape; therefore, they are more readily justified than arbitrarily selected ratios between untransformed measurements for multivariate analysis (Klein et al., 2010). Correlation coefficients between transformed variables and greatest skull length were calculated to check if the data transformation was effective in removing the effect of size (Khosravi, Kaboli, Imani, & Nourani, 2012).

| Data analysis
To identify the combination of variables that best separate skull samples and to define morphologically similar population groups, a principal component analysis (PCA) was performed applying both log size and shape and log-shape variables. We calculated the correlation coefficient between cranial scores and the geometric means for log size and shape variables and log-shape variables in the first four components. The relationship between size and shape was inspected by calculating the correlation between each of the log-shape measurements and the geometric mean for all samples. Finally given the geographic distribution of the specimens, we assigned them into nine groups and projected observations position on two first PCs using the geographic prior. These geographic groups include Arabian P, WZIran, EZIran, Khark, Borazjan, Caucasus, CAsia, China, and Mongolia. Assigning specimens to the nine groups was done regarding the biogeographic structuring, geographic morphology of the region, and recent findings of Fadakar et al. (2020). We performed PCA and the associated graphs in R using "ade4" and "ggplot2" packages, respectively.

| Ecological niche modeling procedure
To examine ecological divergence between gazelle groups, we adopted a between-taxa species distribution modeling (SDM) based on environmental variables and a between-group niche comparison based on a PCA analysis on climatic variables. For the first procedure, we conducted a maximum entropy (MaxEnt) model to predict highly suitable habitats for the three taxa of G. marica (Wacher et al., 2011), G. s. subgutturosa, and G. s. yarkandensis (Fadakar et al., 2020). To do so, four categories of environmental variables including land cover, anthropogenic, topographic, and climatic were compiled with occurrence points to conduct gazelles habitat suitability model. First, a total of 280 presence points were obtained from a variety of sources including direct field surveys from 2014 to 2018 by the authors, museum-based collections, extraction from Global Biodiversity Information Facility (GBIF), and scientific literature review (Table S3). For those with no coordinates but exact locality names, records were georeferenced using the global gazetteer version 2.3 (http://www.falli ngrain.com/world). The reliability of all records was assessed by mapping them in Google Earth version 7.1. We did not include gazelle populations in the contact zone of G. subgutturosa and G. marica in western Iran and Iraq as either we did not have skull specimen from them, or their genetic status was unknown. Moreover, we excluded Khark population and Borazjan gazelle from the ENM analysis, as the origin (native or introduced) of gazelles in Khark, and the locality of Borazjan gazelle was unknown.
To prepare environmental variables for MaxEnt analysis, we extracted four land cover types of the mosaic of herbaceous with sparse tree and shrub, sparse herbaceous, consolidated land (e.g., hardpans, gravels, and rocks), and unconsolidated land (e.g., bare soils and shifting sands) from the Global Land Cover by National Mapping Organization (GLCNMO) version 3 (Kobayashi et al., 2017).
We then calculated the proportion of each cover type in a 7 × 7 grid cell moving window neighborhood analysis. To incorporate anthropogenic impact, we used human footprint as a measure of human influence on the land surface based on the data derived from the 2009 Human Footprint (Venter et al., 2018). Using the Shuttle Radar Topography Mission (SRTM) elevation model (http://srtm.csi.cgiar. org), the two most important topographic variables were compiled: elevation and topographic roughness, that is, standard deviation (SD) of the elevation of all raster cells included in the 7 × 7 grid cells.
Modeling was repeated 10 times based on a cross-validation method, and the predictive performance of models was evaluated based on the area under the receiver operating characteristic curve (AUC). The pairwise spatial overlap between derived habitat suitability models was then calculated based on Schoener's D niche overlap using ENMTools (Warren, Glor, & Turelli, 2010). Schoener's D ranges from 0 to 1 and represents a gradient of complete dissimilarity to fully overlapping niches. were used to test the hypotheses of niche divergence. The available niche space for the gazelle groups was defined as all pixels of the 19 climatic variables within a buffer of 50-km enclosing the species presence points.

| Morphometric analysis
The number of specimens for each location and some descriptive statistics of 14 skull measurements are presented in Table S4.

| Principal components analysis based on log size and shape variables
In of braincase were most highly correlated to PC1, and greatest width across horns, horn length, and distance between tips of horns were most highly correlated to PC2. We address only the first two components because these components contain most of the information on how well log size and shape variables separate skulls from different regions. Figure 2 plots the first and second principal component scores against each other; those groups that do not overlap on PC1 differ significantly in size (mainly of the skull), while those that do not overlap on PC2 differ significantly in shape (especially of the horns). Figure 2 shows that size and shape, singly or together,

| Principal components analysis based on logshape
In

| The interaction of size and shape
The results of the correlation between the principal component scores and the geometric means for log size and shape variables in the four first components showed that there is a significant correlation between the first component and the geometric mean (Table 2).
Since in PCA, the first component reflects the size, it was expected that the correlation between component 1 and the geometric mean would be tight ( Figure S3). The correlation between the scores of the samples based on cranial variables and the second component was low and insignificant ( Figure S3; Table 2). Correlation between scores of individuals' log-shape variables and the geometric mean showed that no components were correlated significantly with size (the geometric mean), implying that skull size has no impact on skull shape (Jungers et al., 1995;Klein et al., 2010). Table 3 presents the correlation coefficients between the individual log-shape variables and corresponding geometric means. A coefficient near 0 indicates isometry, a positive coefficient indicates positive allometry, and a negative coefficient indicates negative allometry (Klein et al., 2010;Mosimann & James, 1979). The measured variables tend to increase more slowly than the overall size (Table 3).
Horn variables (DBTH and GWH) showed a significant positive correlation, but HL also showed a non-significant negative correlation, meaning that horn variables tend to increase more rapidly than the overall size. Skull length showed a significant negative correlation and appears to increase particularly slowly.

| Ecological niche divergence
Notable consistency was detected in the predicted habitat suitability of gazelles when compared with occurrence records ( Figure 3)  Considering response curves of the environmental variables, we found that G. marica prefers lower-elevated habitats mostly covered by unconsolidated lands with higher temperature and lower temperature seasonality compared to G. s. subgutturosa and G. s. yarkandensis ( Figure S4).
For between-group niche comparisons, PC1 and PC2 of the climatic variables explained 49.6% and 20.8% of the variation in climatic variables, respectively. PC1 was mostly correlated with mean variables of temperature and precipitation while PC2 was highly correlated with standard deviation of temperature (temperature seasonality;  Table 4).
Background similarity test also showed a significant niche difference between WZIran and EZIran. Arabian P-WZIran niche comparison indicated that there is a significant difference (p < .05) in niche overlap between WZIran observed niche and niche selected randomly from the Arabian P. For pairwise comparisons of EZIran with other eastern gazelles (Caucasus, CAsia, China), we found no significant niche differences. While niche comparison of the Caucasus-China was significantly different, we found no significant background difference in the Caucasus-CAsia and CAsia-China niche comparisons (Table 4). Results also showed a complete divergent niche occupation (D = 0) between gazelles of Mongolia and other Asiatic gazelles, except for the China-Mongolia niche comparison for which niche overlap was very low (D = 0.08) and their background niche was significantly different (p < .05; Figure 5 and Table 4).
Some populations occurring in EZIran and CAsia are somehow separated either in size or shape from each other. According to Groves and Grubb (2011), gazelles in western Turkmenistan confirming the difference between gazelles of these two regions.
Our morphometric analyses separated Mongolian, CAsia, and EZIran populations from each other based on shape, but the latter two showed considerable overlap in size. The implication is that  (Table S1) Our analyses also separate the Khark Island population from other populations in size and somehow shape. A recent genetic study has suggested that the Khark Island gazelle is closer to G. subgutturosa than to G. marica (Mirzakhah et al., 2015). Nevertheless, the Khark Island population showed a strong connection with the closest population on the mainland both in geographic distance and size, that is, the presumably extinct Borazjan gazelle population named G. karamii, a taxon that its affiliation to G. marica has previously been verified (Bärmann et al., 2013). Male-biased gene introgression between these two species (Murtskhvaladze, Gurielidze, Kopaliani, & Tarkhnishvili, 2012) (1912), and so brighter than in the gazelle in south-western Iran.

| Ecological niche divergence
Our results indicated a significant niche divergence between gazelles of EZ and WZ. This high and long mountainous belt has been acting as a serious barrier toward dispersal of species with Sahara-Arabian origin into the central parts of Asia (Lerp, Wronski, Butynski, & Plath, 2013;Mouthereau, 2011). However, the penetration of sand gazelle's genes into EZIran has been revealed by detecting two haplotypes of G. marica in central and north of Iran (Fadakar et al., 2020).
Our result showed that the path may have been from south-west to south and south-east Iran into the central and northern parts of the country. Nevertheless, the rate of gene flow from WZ to central Iran has been very low, as despite the existence of suitable habitat in central Iran for G. marica, no population of this species has been established. In addition, western Pakistan and south-west of Afghnistan appeared suitable for G. marica; hence, the haplotypes of this species may also be detectable in these two countries as there is considerable habitat connectivity for the species in southern Iran.
The same paths may have allowed the sand gazelle's ancestors to bypass the Zagros Mountains into Arabia and North Africa, and subsequently diverge into new species . In EZ, a strong disconnectivity in suitable ranges of Central Asia and China was detected along the Pamir and Tian Shan Mountains. This finding implies that there might be cryptic diversification patterns in the region, which is also supported by our morphometric results.
Our habitat suitability models showed that all the studied gazelle groups are associated with arid steppes covered by consolidated and or unconsolidated lands with sparse vegetation, as also reported by Lerp et al. (2016). However, we found that they show an apparent disparity in their climatic niche space. Annual, monthly, and quarterly temperature and precipitation showed strong negative relationships with PC1. We found that along with this PC, the ecological niche of WZ gazelles (i.e., G. marica) is limited to warmer and drier condi-  (Lerp et al., 2016). Generally, species with greater dispersal abilities occupy larger ranges which in turn allow them better expansion into available space (Pavlicev & Mayer, 2009). Due to the high dispersal ability of large to medium-sized species, in this case gazelles, gene flow might yet be maintaining, hampering the formation of distinct genetic structures (Ashrafzadeh, Khosravi, Ahmadi, & Kaboli, 2018), and facilitating hybridizations in contact zones (Fadakar et al., 2019). This pattern has been suggested for gazelles of Euphrates-Tigris basin to western borders of Iran, where populations of G. subgutturosa and G. marica represent a mix of morphometric characteristics (Wacher et al., 2011). Moreover, the lack of significant topographic heterogeneity in the landscape of Central Asia might have led to a non-significant divergence in the ecological niche.  (Fadakar et al., 2020).

| CON CLUS ION
Intraspecific classification of goitered gazelle has been prone to change as enough genetic, morphometric, and ecological data have not yet been available for a final conclusion. For most of the studied populations, female specimens were not available; therefore, we only used male skulls to investigate taxon differentiation. In contrast, the supportive genetic data available were primarily based on mitochondrial DNA, which is inherited maternally. We therefore incorporated additional ecological analyses (climatic niche separation and habitat suitability) to mitigate this deficiency and strengthen the morphometric results.
Nevertheless, the integration of morphometric and ecological analyses of this study with those obtained from previous genetic studies is still unconvincing and does not allow for definitive classification. The whole picture may be compromised by incomplete and opportunistic sampling, small samples sizes (particularly for our morphometric data), and the unreliability of mitochondrial markers compared to DNA fingerprinting for phylogenetic studies at the intraspecific level (Ingman, Kaessmann, Pääbo, & Gyllensten, 2000;Wan, Wu, Fujihara, & Fang, 2004). However, the notable barriers restricting gene flow among geographic population groups of goi-

ACK N OWLED G M ENTS
This paper is dedicated to the memory of our great co-author, Professor Colin Groves. We are grateful to all institutions that provided us with skull samples, to the National Iranian Oil Company for permission for MRH to enter Khark Island, and to Davoud Fadakar for providing helpful information. We are also grateful to Hannes Lerp and an anonymous reviewer for their helpful comments on an earlier version of the manuscript.

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