Divergent selection along elevational gradients promotes genetic and phenotypic disparities among small mammal populations

Abstract Species distributed along mountain slopes, facing contrasting habitats in short geographic scale, are of particular interest to test how ecologically based divergent selection promotes phenotypic and genetic disparities as well as to assess isolation‐by‐environment mechanisms. Here, we conduct the first broad comparative study of phenotypic variation along elevational gradients, integrating a large array of ecological predictors and disentangling population genetic driver processes. The skull form of nine ecologically distinct species distributed over a large altitudinal range (100–4200 m) was compared to assess whether phenotypic divergence is a common phenomenon in small mammals and whether it shows parallel patterns. We also investigated the relative contribution of biotic (competition and predation) and abiotic parameters on phenotypic divergence via mixed models. Finally, we assessed the population genetic structure of a rodent species (Niviventer confucianus) via analysis of molecular variance and FST along three mountain slopes and tested the isolation‐by‐environment hypothesis using Mantel test and redundancy analysis. We found a consistent phenotypic divergence and marked genetic structure along elevational gradients; however, the species showed mixed patterns of size and skull shape trends across mountain zones. Individuals living at lower altitudes differed greatly in both phenotype and genotype from those living at high elevations, while middle‐elevation individuals showed more intermediate forms. The ecological parameters associated with phenotypic divergence along elevation gradients are partly related to species' ecological and evolutionary constraints. Fossorial and solitary animals are mainly affected by climatic factors, while terrestrial and more gregarious species are influenced by biotic and abiotic parameters. A novel finding of our study is that predator richness emerged as an important factor associated with the intraspecific diversification of the mammalian skull along elevational gradients, a previously overlooked parameter. Population genetic structure was mainly driven by environmental heterogeneity along mountain slopes, with no or a week spatial effect, fitting the isolation‐by‐environment scenario. Our study highlights the strong and multifaceted effects of heterogeneous steep habitats and ecologically based divergent selective forces in small mammal populations.

of phenotypic variation along elevational gradients, integrating a large array of ecological predictors and disentangling population genetic driver processes. The skull form of nine ecologically distinct species distributed over a large altitudinal range (100-4200 m) was compared to assess whether phenotypic divergence is a common phenomenon in small mammals and whether it shows parallel patterns. We also investigated the relative contribution of biotic (competition and predation) and abiotic parameters on phenotypic divergence via mixed models. Finally, we assessed the population genetic structure of a rodent species (Niviventer confucianus) via analysis of molecular variance and F ST along three mountain slopes and tested the isolationby-environment hypothesis using Mantel test and redundancy analysis. We found a consistent phenotypic divergence and marked genetic structure along elevational gradients; however, the species showed mixed patterns of size and skull shape trends across mountain zones. Individuals living at lower altitudes differed greatly in both phenotype and genotype from those living at high elevations, while middle-elevation individuals showed more intermediate forms. The ecological parameters associated with phenotypic divergence along elevation gradients are partly related to species' ecological and evolutionary constraints. Fossorial and solitary animals are mainly affected by climatic factors, while terrestrial and more gregarious species are influenced by biotic and abiotic parameters. A novel finding of our study is that predator richness emerged as an important factor associated with the intraspecific diversification of the mammalian skull along elevational gradients, a previously overlooked parameter. Population genetic structure was mainly driven by environmental heterogeneity along mountain slopes, with no or a week spatial effect, fitting the

| INTRODUC TI ON
In the era of rapid changes, predicting how species will respond to human-mediated habitat alterations has become a 21st-century scientific pursuit, uniting ecologists, evolutionary biologists, and conservationists (Fox, 2018). Anticipated climatic changes are expected to vary widely, leading to novel habitat-animal interactions and divergent selective forces (Bellard, Bertelsmeier, Leadley, Thuiller, & Courchamp, 2012;Keller, Alexander, Holderegger, & Edwards, 2013). Regarding this scenario, elevational gradients offer one of the best natural systems for understanding species performance under heterogeneous environmental constraints (Keller et al., 2013).
Species distributed along mountain slopes, which face conspicuous climatic and steep vegetation changes on a short geographic scale, are of particular interest for testing how ecologically based divergent selection promotes phenotypic divergence and shape genetic structure across populations (Rosenblum, 2006). Divergent natural selection, briefly defined as selection arising from uneven ecological conditions that promote phenotypic and genetic shifts among populations, is a key source of biological diversity (Nosil, Harmon, & Seehausen, 2008). As altitude increases, steep environmental changes promote spatially dynamic biotic and abiotic interactions. Such mixed effects of contrasting climatic variables on a short scale reflect marked heterogeneous vegetation zonation along mountain slopes, in some cases mirroring tropical-temperate forest transitions. For example, in southeastern Asian mountains, evergreen broadleaf forests occur at 1,000 m, and alpine landscapes occur at 3,500 m (Ohsawa, 1991(Ohsawa, , 1993. Individuals living in distinct elevation zones are thus expected to show phenotypic and genetic disparities owing to strong divergent ecological selection (Branch, Jahner, Kozlovsky, Parchman, & Pravosudov, 2017;Kawecki & Ebert, 2004;Keller et al., 2013;Ohsawa & Ide, 2008). In contrast, elevational gradients over a small geographic scale may facilitate migration, especially in species with a high dispersal capacity, thereby constraining morphological and genetic divergence (Branch et al., 2017;Kawecki & Ebert, 2004;Keller et al., 2013;Nosil et al., 2008;Ohsawa & Ide, 2008;Sexton, Hangartner, & Hoffmann, 2014;Waterhouse, Erb, Beever, & Russello, 2018). Thus, a tradeoff between the strength of divergent selective forces and dispersal patterns leads to distinct scenarios: a generalist phenotype and a panmictic population distributed across elevation zones or divergent phenotypes and structured genetic populations in each elevation zone (Allendorf & Luikart, 2007;Kawecki & Ebert, 2004;Nosil et al., 2008;Ohsawa & Ide, 2008).
Previous studies have pointed out that phenotypic and genetic divergence between animal populations along elevational gradients is not a rare phenomenon (reviewed by Keller et al., 2013).
Intraspecific morphological shifts have been detected in numerous species of arthropods, lizards, frogs, and even small birds and mammals (Branch et al., 2017;Caro, Caycedo-Rosales, Bowie, Slabbekoorn, & Cadena, 2013;Grieco & Rizk, 2010;Keller et al., 2013). However, whether this pattern is pervasive across groups and traits remains unclear because the majority of related studies have focused on low-vagility animals (Keller et al., 2013) and phenotypic features with clear adaptive value along altitudinal gradients (e.g., cold and desiccation tolerance).
While morphological shifts have been widely assessed, their associated factors remain mostly speculative. Studies on the effects of climatic variation have mainly focused on temperature-size correlations, testing the validity of Bergman's and Allen's rules along altitudinal gradients (Brehm & Fiedler, 2004;Du et al., 2017;Freeman, 2017;Gutiérrez-Pinto et al., 2014). How competition and predation shape species morphology across elevation zones remains an open question. Increasing altitude leads to a reduced area, enhanced spatial isolation, and relaxed predation risk, similar to the patterns observed in insular systems (Körner, 2007). Considering well-established island hypotheses, we would expect biotic interactions to play an important role in shaping morphological diversity in mountains as they do in islands (Heaney, 1978;Itescu et al., 2018;Michaux, Bellocq, Sarà, & Morand, 2002). To date, no study has investigated the combined effect of biotic and abiotic traits on phenotypic divergence on mammals along elevational gradients (Keller et al., 2013;Merilä & Hendry, 2014).
Further insights into divergent ecological selection can be derived from landscape genetic studies (Allendorf & Luikart, 2007;Kawecki & Ebert, 2004;Merilä & Hendry, 2014;Ohsawa & Ide, 2008). The presence of animals (particularly those with high vagility) in contrasting habitats on a short geographic scale makes such animals attractive candidates with which to test for ecologically induced genetic alterations. Isolation by environment (IBE) applies when the genetic variation among populations is mostly explained by climatic factors rather than spatial factors (Wang & Bradburd, 2014;Wang & Summers, 2010). Thus, the existence of limited gene flow between mammalian populations across mountain zones would indicate a strong role of heterogeneous landscapes in shaping the genetic structure of natural populations (Branch et al., 2017;Kierepka & Latch, 2015;Wang & Bradburd, 2014). isolation-by-environment scenario. Our study highlights the strong and multifaceted effects of heterogeneous steep habitats and ecologically based divergent selective forces in small mammal populations.

K E Y W O R D S
altitude-induced adaptation, isolation-by-environment, landscape ecology, mountain biota, phenotypic divergence, predation In this study, we aimed to obtain a comprehensive understanding of the effects of divergent selective forces along elevational gradients on genetic and phenotypic disparities in small mammals through three complementary steps. First, we analyzed nine species (shrews, moles, and rodents) distributed over a wide altitudinal range and with distinct ecological attributes to assess whether phenotypic divergence is a common phenomenon in small mammals and whether it shows consistent parallel patterns. We chose skull morphology (size and shape) as a surrogate of overall morphological diversity as it is associated with the main sensory systems related to environmental perception and reflects changes in the dietary and ecological interactions across populations (Klaczko, Sherratt, & Setz, 2016;Monteiro, Lessa, & Are, 1999;Nogueira, Peracchi, & Monteiro, 2009). Second, using a large array of biotic and abiotic parameters, we investigated the ecological factors associated with the intraspecific diversification of the mammalian skull along elevational gradients. Our goal here was to uncover the relative contribution of abiotic and biotic traits associated with phenotypic variation across ecologically distinct species. Finally, we assessed the population genetic structure of a widespread rodent species along three mountain slopes of southwest China and tested the IBE hypothesis, thereby disentangling the effects of climatic and spatially driven factors.
Given the dynamic habitat-animal interactions and strong divergent selection along elevational gradients, we hypothesized that species would show marked phenotypic disparity between populations inhabiting distinct altitude zones, regardless of their phylogenetic and ecological features. In contrast, considering species' inherent life-history attributes, we predicted inconsistent main drivers of morphologic divergence across taxa (e.g., that desert dwellers would be more affected by variation in precipitation than would forest dwellers). Because of the large home range and generalist habits of our focal species, we expected weak or absent genetic structure along mountain slopes.

| Focal species and phenotypic traits
Nine mammal species (six rodents, two shrews, and one mole) distributed over a large altitudinal range and representing distinct ecological niches were selected to assess phenotypic variation (size and shape of the skull) across elevational gradients (Table 1). All specimens are housed in the mammal collection of the Institute of Zoology of the Chinese Academy of Science, Beijing, China. Specimens were sampled from the same mountain (Soriculus) or along altitudinal gradients with similar climatic conditions. In the last case, to minimize putative effects of distinct evolutionary history, we selected individuals from closely genetic lineages Ge et al., 2018;He, Hu, Chen, Li, & Jiang, 2016;Liu et al., 2018;Lv et al., 2018).
The phenotype was abstracted from the shape and size of the skull of 494 adult individuals using geometric morphometric techniques. A list of the specimens studied is provided in Appendix S1.
The ventral view of the skull of adult specimens of a similar age TA B L E 1 Species of shrews, rodents, and moles used in this study, showing the sample size per elevation zone, species total elevation range, elevation range of the samples, and natural history attributes  (Mitteroecker, Gunz, Windhager, & Schaefer, 2013) and is used here as an indicator of individual size. No sexual dimorphism in skull shape or size was found in the nine species examined (Table S3); therefore, we pooled males and females in the following analyses.

| Ecological parameters
We quantified the effects of climate and biotic interactions on skull phenotypic changes of mammals using a large array of biotic and abiotic predictors. The topographical dataset was obtained from the WorlClim website and generated from NASA's Shuttle Radar Topography Mission (https ://www2.jpl.nasa.gov/srtm/). Based on the geographic coordinates of each individual, we extracted the altitude information via the raster R package (Hijmans, 2017). Five climate traits were used: annual mean temperature, temperature seasonality, annual precipitation, precipitation seasonality, and net primary productivity (NPP). Net primary productivity was obtained from the NASA Socioeconomic Data and Applications Center (Imhoff & Bounoua, 2006;Imhoff et al., 2004).
To assess the effects of competition (intra-and interspecific), we used the database of sampling collections of the mammal group of the Institute of Zoology of the Chinese Academy of Science. This unique database includes specimens of small mammals systematically collected from several mountains in China. Detailed descriptions of the sampling protocols can be found in Wu et al. (2013) and Wen et al. (2018). For each specimen's collection locality, we computed the number of individuals of a given selected species and the species richness recorded in a 1-km buffer. We used these two parameters as proxies of intraspecific and interspecific competition, respectively. Although this counting method (Dueser & Hallett, 1980;Krebs, 1966) likely underestimates the overall abundance and local richness of small mammals, the results from standardized sampling protocols applied by our group can still be used in a comparative way (Krebs et al., 2011;Pacheco et al., 2013).
To evaluate the influence of predation risk on phenotypic variation, we compiled a list of predators of small mammals (i.e., carnivorous mammals, snakes, and predatory birds) in China (Hoyo, Elliott, & Sargatal, 1999;Smith & Xie, 2013

| Molecular sampling
To explore the genetic structure of small mammal populations distributed along elevational gradients, we selected the Confucian white-bel-  . All sequences are available from GenBank (Table S5).
On Gongga Mountain, we obtained complete cytb sequences (1,132 bp) from 160 individuals comprising six populations covering an altitudinal range of 2,000 m ( Figure 2). The median Euclidian distance between sampling localities was 3.2 km (range: 2.69-6.9 km).
On Luoji Mountain, we obtained cytb sequences from 69 specimens representing five populations covering an altitudinal range of 1,600 m, and the median Euclidian distance between sampling localities was 1.6 km (range from 1.34-6.8 km). On Wolong Mountain, we obtained cytb sequences from 23 individuals comprising six populations distributed at an altitudinal range of 1,500 m with a median Euclidian distance between collecting localities of 11.2 km (range from 0.62-19 km; Table 2). For all mountains, we divided our sample into three elevation zones (low, middle, and high) as detailed in

| Phenotypic comparison across elevation zones
To evaluate the phenotypic divergence across elevational gradients, we divided our samples into three bioclimatic zones (low, middle, and high elevation). The delimitation of zones combined climatic and vegetation attributes and was based on Tang and Ohsawa (1997), Zhong, Zhang, and Luo (1999), and Li and Zhang (2010). The elevational range of each zone might vary among mountains as a reflection To test for phenotypic divergence among elevation zones, we conducted analyses of variance using Procrustes coordinates for shape and centroid size for size in the geomorph R package with 10,000 permutations. In addition, we carried out pairwise comparisons between mean shapes and performed Tukey's honestly significant difference (HSD) test to assess the significance of morphological differences between pairs of elevation zones. We also explored size and shape trends across species to determine whether there was a consistent overall pattern. The relationship between size (centroid size) and altitude was visualized with scatterplots highlighting trends in each elevation zone. Skull shape deformation between zones was assessed via between-group principal component analysis (BgPCA; Mitteroecker & Bookstein, 2011) using Morph R package (Schlager, 2017).

| Phenotypic-ecological interaction analyses
We assessed the combined effects of our ecological parameters (and quadratic terms when the relationship was nonlinear) on skull phenotypes by creating a set of competing models ranked by the secondorder corrected Akaike information criterion (AICc) via the MuMln R package (Barton, 2018) and tested them by using mixed hierarchical models in the lme4 R package (Bates, Maechler, Bolker, & Walker, 2015). In all models, site was included as a random effect to account for uneven sample sizes among collection localities (Harrison et al., 2018;Paterson, 2003;Schielzeth & Forstmeier, 2008). We prevented overparameterization by limiting the number of predictors in a single model based on the sample size (Burnham & Anderson, 2002;Grueber, Nakagawa, Laws, & Jamieson, 2011;Harrell, 2001;Harrison et al., 2018). The importance of each parameter in models with similar likelihoods was summarized using the full model-average approach (Burnham & Anderson, 2002;Grueber et al., 2011) considering only models with a ∆AICc ≤ 6 following Richards (2008).
Summarizing shape information in univariate terms might lead to a loss of information (Fontaneto et al., 2017). Therefore, we use two variables to provide a more complete description of shape in the mixed analyses. The first univariate term consisted of the score of regression of Procrustes coordinates on the log of centroid size, which has been shown to be an effective descriptor of intraspecific shape variation (Drake & Klingenberg, 2008;Maestri et al., 2016;Martinez et al., 2018). The second variable was the first principal component of the Procrustes coordinates, which summarized most of the shape variation. In this case, to account for the allometry effect, we added the scaled centroid size as an independent covariate in our models.

Spatial autocorrelation is an inherent issue in interpopulation
studies that can potentially bias statistical parameters and inflate the type I error (Chun & Griffith, 2011;Diniz-Filho, Bini, & Hawkins, 2003;Perez, Diniz-Filho, Bernal, & Gonzalez, 2010). To account for spatial dependence in our statistical models, we applied the eigenvector spatial filtering method (detailed in Griffith, 2003), which introduces a set of positive eigenvectors (spatial filters) extracted from the spatial weight matrix as fixed factors into regression models. The spatial filters were obtained via the spdep R package (Bivand & Piras, 2015) following the protocol described by Griffith and Chun (2014).
Finally, we ran Moran's I test to ensure that the final residuals of our mixed models lacked significant spatial autocorrelation.

| Population genetic analyses
DNA sequences were aligned using the ClustalW algorithm implemented in MEGA 7.0 (Kumar, Stecher, & Tamura, 2016). To assess the patterns of genetic diversity among elevation zones, we calculated the nucleotide diversity, the mean number of pairwise differences, the number of haplotypes, and the haplotype diversity for each population using DnaSP 6.11.01 (Rozas et al., 2017). We tested the genetic difference among populations based on haplotype frequencies via a chi-squared test with 1,000 permutations implemented in DnaSP. Population genetic structure along elevation zones was measured via analysis of molecular variance (AMOVA) based on haplotype diversity (Excoffier, Smouse, & Quattro, 1992) using Arlequin 3.5.2 software (Excoffier & Lischer, 2010), and the significance was assessed via 10,000 permutations. In a fine-scale approach, we calculated the F ST between pairs of populations ( Figure 2) to explore To investigate the correlation between genotype and spatial environmental components and assess whether the putative genetic structure was caused by random (isolation by distance) or nonrandom (IBE) mechanisms, we performed two complementary tests (Kierepka & Latch, 2015). First, we carried out a Mantel test with 1,000 permutations to explore the correlation between genetic and geographic distances for each mountain using the vegan R package (Oksanen et al., 2018). To disentangle the effects of spatial and climatic variables on the landscape genetic structure along mountain slopes, we carried out a redundancy analysis (

| Phenotypic comparison across elevation zones
Significant cranial differences among elevation zones were detected in all species examined (Table 3). Decomposition of skull form components revealed significant shape differences among the zones in eight species, while size diverged in five species. However, elevational zonation tended to explain a larger proportion of variance in size (9%-49%) than in shape (3%-19%; Table 3). Further inspection using pairwise and Tukey's HSD tests revealed that animals inhabiting low altitudes were differentiated from those living at middle and high altitudes to a greater extent than individuals from middle altitudes were differentiated from those at high elevations (Table S6).   (Figure 4).

| Ecological factors associated with phenotypic changes along elevational gradients
Hierarchical regressions of skull phenotype on ecological predictors uncovered an inconsistent combination of factors that best explained phenotypic variation across taxa (Table 4). The skull of shrews was mainly affected by environmental traits, whereas rodents and moles responded to both biotic and abiotic interactions. Size variation along elevational gradients was explained by equally important multiparameters in most cases, while skull shape variation was correlated with one variable or a few variables (Table 4). Altitude was among the top variables in the selected models in terms of explaining size variation in all species examined, and its effects varied from weak to strong across taxa (Table 4). Among the biotic variables, predator richness had the most widespread effect across species and, in most cases, had a positive relationship with skull size (Table 4).

| Population genetic structure along mountain slopes
On all three mountains, the frequency of haplotypes (Gongga Mt.:  Table S8). Based on comparisons of pairs of populations sampled along mountain slopes, middleelevation populations had the highest nucleotide diversity (Table S9).

F I G U R E 4
Between-group PCA of skull shape (Procrustes coordinates) per elevation zone with 95% confidence ellipses of shrews (a, b), moles (c), and rodents (d-i). Information on habitat, diet, and life mode of each species is shown in colored polygons. Deformation grids show skull shapes at the extremes of each axis TA B L E 4 Estimates (and percentage of relative importance) for each parameter from the full averaged models of biotic and climatic parameters on skull phenotype (size and shape) of small mammals along elevational gradients

*
Note: Bold face indicates parameters considered significant, when confidence intervals (95%) did not include zero. Full statistics are presented in Table S7. All traits were standardized prior to analyses. Spatial filters were included as fixed factors to account for spatial autocorrelation. Inter. Comp.: interspecific competition; Intra. Comp.: intraspecific competition; Main pred. richn.: main predator richness (species that have small terrestrial mammals as their primary diet); NPP: net primary productivity; Prec., precipitation; Pred. richn.: predator richness; Seas.: seasonality; Temp.: temperature. Asterisks indicate that the quadratic term of the parameter was used.  (Table S10). The genetic structure of N.
confucianus was mainly driven by climatic variables (Table 5). The amount of genetic variance explained by our set of environmental parameters (up to 37%) was significant across all mountains, even when the spatial component was controlled for. In contrast, spatial parameters had either no effect (in Luoji Mountain: Mantel r = 0.02, p = 0.37) or a weak effect on the landscape genetic structure of N.

| D ISCUSS I ON
By comparing diverse phylogenetic and ecological groups of small mammals, we found consistent phenotypic disparity in cranial traits and marked genetic structure along elevational gradients. This finding, combined with previous evidence of life-history and physiologically changes along mountain slopes (Keller et al., 2013), reinforces the strong and multifaceted effects of heterogeneous steep habitats on individuals' features.
The phenotypic pattern of individuals living at lower elevations differed greatly from that of individuals living at high elevations, which paralleled our genetic results. We can speculate that this morphogenetic divergence may be related to the pronounced ecological differences faced by populations at opposite ends of an elevational gradient (Caro et al., 2013;Kawecki & Ebert, 2004;Wang & Bradburd, 2014). This interpretation also explains the lower morphological and genetic differentiation of middle-elevation individuals as they inhabit less contrasting environmental conditions and thus are exposed to weaker divergent selective forces (Allendorf & Luikart, 2007;Kawecki & Ebert, 2004). In addition, animals at middle altitudes are in closer contact with both low-and high-altitude populations, which may favor intermediate forms able to exploit distinct elevation zones. In contrast, high-altitude animals with higher levels of differentiation are less likely to colonizing other zones (Kawecki & Ebert, 2004). Although further studies focusing on the genes linked to the observed phenotypic changes are still needed (Merilä & Hendry, 2014), the parallel results obtained from the genetic and morphological approaches suggest that evolutionary mechanisms rather than plasticity alone are likely to contribute to the intraspecific diversification of the mammalian skull along elevational gradients (Boutin & Lane, 2014;Merilä & Hendry, 2014).
The higher genetic diversity of middle-altitude populations of N. confucianus might reflect optimal ecological conditions faced by these individuals (Ohsawa & Ide, 2008). This hypothesis predicts that for species with elevational range centers at middle altitudes, as N. confucianus Wen et al., 2017), peripheral populations at low and high altitudes experience suboptimal conditions, which reduce their niche breadth, their population size, and ultimately their genetic variation (Heino & Gronroos, 2014;Ohsawa & Ide, 2008;Wen et al., 2018). Therefore, we expect that species with distinct elevational range centers will show different patterns of genetic diversity along mountain slopes. Nevertheless, the absence of an effect or a weak effect of spatial traits on the population genetic structure of N. confucianus along the three mountain slopes points to more general patterns for small mammals, where nonrandom factors are the main driver of genetic variation (Garant, Forde, & Hendry, 2007;Keller et al., 2013;Kierepka & Latch, 2015). Accordingly, RDA revealed that environmental variables explained up to 17% of the genetic variance when geography was controlled for. These results fit an IBE scenario, especially considering that the distance between the genetically sampled populations was shorter than the individual dispersal capacity (Kierepka & Latch, 2015;Sexton et al., 2014;Wang & Bradburd, 2014).
Skull shape divergence appears to be a more common phenomenon than size divergence in mammals distributed across a wide elevational range, although, when present, the difference in size tends to be more pronounced (Table 3). One possible explanation for this difference is that skull shape, as a proxy of dietary and life-mode adaptations (Klaczko et al., 2016;Monteiro et al., 1999;Nogueira et al., 2009), is under strong selection to enhance individual fitness in TA B L E 5 Results of redundancy analysis (RDA) models to assess the contribution of climatic and spatial variables on Niviventer confucianus genetic structure along the mountain slopes in China a given environment (Myers, Lundrigan, Gillespie, & Zelditch, 1996).
For example, intraspecific changes in cranial masticatory muscle attachments have been linked to varying levels of insectivory in armadillos as a morphogenetic response to the availability of food resources in heterogeneous habitats (Feijó, 2017). Myers et al. (1996) showed that changing only the consistency of the food (hard and soft diets) was enough to promote intraspecific cranial and dental divergence in rodents. Likewise, Fornel, Cordeiro-Estrela, and Freitas (2010) found that rodents inhabiting distinct habitats (sand dunes and fields) show significant skull shape differences. Alternatively, while skull shape is a three-dimensional structure and thus can display diverse phenotypes, size can change only by increasing or decreasing. The dynamic spatial interaction between abiotic and biotic forces along elevational gradients might generate competing pressures (see Table 4) on size changes and thereby minimize the differences between populations across elevation zones (Kawecki & Ebert, 2004;McCain & Grytnes, 2010).
The general change in size as a function of altitude exhibited a nonpositive (negative or neutral) pattern across the nine species. This finding is consistent with that of previous studies on animals along elevational gradients at intraspecific (Caro et al., 2013;Eastman, Morelli, Rowe, Conroy, & Moritz, 2012;Liao, Zhang, & Liu, 2006;Wasserman & Nash, 1979) and interspecific (Du et al., 2017;Geraghty, Dunn, & Sanders, 2007;Hu, Xie, Li, & Jiang, 2011) levels. Conversely, a positive relationship between size and altitude, as expected under Bergmann's rule, has also been reported in the literature (Keller et al., 2013). Such conflicting patterns might reflect idiosyncratic responses due to inherent ecological, physiological, and evolutionary constraints (Adams & Church, 2008;Itescu et al., 2018). Nonetheless, the nonpositive trend shared by moles, shrews, and rodents points to a more common scenario. Sampling part of a species' (altitudinal) range may fail to reveal the real trend, as shown by Adams and Church (2008)  At high elevations, the reduced species richness  and predation pressure (Fu et al., 2007;Kumar, Longino, Colwell, & O'Donnell, 2009; (Liao et al., 2006), while that of the high-altitude plateau pika (Ochotona curzoniae) presents a positive trend (Lin, Ci, Zhang, & Su, 2008). Regarding skull shape, we found a lack of consistency in changes across elevation zones, even between sympatric species, such as Allactaga sibirica and Dipus sagitta, revealing that different phenotypes may be equally effective in similar habitats. This manyto-one mapping scenario predicts that complex morphological structures, such as the skull, can perform similar functions with divergent forms (Maestri et al., 2017;Renaud et al., 2018;Wainwright, 2007;Wainwright, Alfaro, Bolnick, & Hulsey, 2005). Alternatively, selection for reduced niche overlap in areas with marked seasonality of resources, such as those at high altitudes, might favor the co-occurrence of distinct phenotypes.
The ecological parameters associated with phenotypic divergence along elevational gradients are expected to be tied to species ecological and evolutionary traits (Goodyear & Pianka, 2008;Itescu et al., 2018;Keller et al., 2013;Ohsawa & Ide, 2008). Consistent with this expectation, our models revealed somewhat inconsistent responses between shrews and rodents. The reduced exposure of shrews due to their fossorial habits (Healy et al., 2014) and asocial organization (Churchfield, 1990) might explain the reduced influence of biotic interactions on their phenotypic variability. Curiously, the size and shape of the mole Uropsilus soricipes, a close relative of shrews that lacks adaptation to fossorial life, were significantly associated with interspecific competition. In addition, phenotypic changes in the two semidesert dwellers (A. sibirica and D. sagitta) were mainly associated with seasonal temperature and precipitation, while forest dwellers tended to be affected by multiple parameters.
In mammals, predator richness is traditionally linked to body size variation in insular systems, where relaxed predation is correlated with an increase in size in small species (Heaney, 1978;Michaux et al., 2002). Surprisingly, we found a positive correlation between size and predation risk in most cases (Table 4). This relationship might be an indirect response to reduced prey species abundance at lower altitudes (Fu et al., 2007;Wen et al., 2018), which leads to relaxed competition for food (Heaney, 1978). Notably, our conclusions are derived from a coarse-grained approach using species' ranges to predict local predator richness. Although our two predator richness parameters generated similar results and agreed with results from previous studies showing higher predation risk at low elevations (Fu et al., 2007), future studies exploring predatorprey interactions along mountain gradients using locally collected data are necessary to test our hypothesis that predation risk is an important factor associated with the intraspecific diversification of the mammalian skull along elevational gradients, a previously neglected parameter.
An open question derived from our study is which main biological factor(s) drives phenotypic divergence at each elevation zone.
Answering this question could provide us with a better understanding of the dynamic ecological interactions along altitudinal gradients (Bolnick et al., 2011;Hart, Schreiber, & Levine, 2016;Violle et al., 2012). Because of our limited sample size in each zone, we could not test the effects of our parameters separately. Nevertheless, we anticipate that at lower elevations biotic interactions will play a stronger role in shaping phenotypic changes, whereas at high elevations due to a relaxed predation and lower competition, climatic factors will be the main drivers.

| CON CLUS ION
By comparing phylogenetically and ecologically diverse species of rodents, moles, and shrews, we found that phenotypic divergence in mammalian skulls along elevational gradients is a common phenomenon. The underlying causes of this phenomenon are partly related to species' life-history attributes. Fossorial and solitary animals are mainly affected by climatic factors, while terrestrial and more gregarious species are influenced by biotic and abiotic parameters. Animals living at lower elevations differ greatly in both phenotype and genotype from those living at high elevations, possibly due to the strong divergent selection acting upon populations at opposite ends of elevational gradients. Therefore, moderate selection faced by animals at middle altitudes seems to favor intermediate forms able to exploit distinct elevation zones. Skull size exhibits a nonpositive trend along elevational gradients, although mixed patterns were detected across elevation zones, possibly due to dynamic multi-interactions among selective drivers. Skull shape displayed contrasting changes between sympatric species across elevational gradients, suggesting equally effective phenotypes in similar habitats. Moreover, we detected consistent population genetic structure mirroring the phenotypic results, which was mainly driven by environmental heterogeneity along mountain slopes (with no or a weak spatial effect), fitting the IBE scenario.

ACK N OWLED G M ENTS
We are grateful to the people who helped us to collect the specimens used in this study. We wish to thank Peter Pearman and two anonymous reviewers for providing valuable comments. AF is sup-

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
AF conceived the idea; AF, ZW, and QY designed the research; ZW, JC, DG, LX, and QY conducted the field sampling; AF, ZW, JC, and DG curated the data; and AF analyzed the data and led the writing.

DATA ACCE SS I B I LIT Y
A list of all specimens examined including their museum collection numbers and localities is presented in Appendix S1.