Isolated alpine habitats reveal disparate ecological drivers of taxonomic and functional beta-diversity of small mammal assemblages.

The interpretation of patterns of biodiversity requires the disentanglement of geographical and environmental variables. Disjunct alpine communities are geographically isolated from one another but experience similar environmental impacts. Isolated homogenous habitats may promote speciation but constrain functional trait variation. In this study, we examined the hypothesis that dispersal limitation promotes taxonomic divergence, whereas habitat similarity in alpine mountains leads to functional convergence. We performed standardized field investigation to sample non-volant small mammals from 18 prominent alpine sites in the Three Parallel Rivers area. We estimated indices quantifying taxonomic and functional alpha- and beta-diversity, as well as beta-diversity components. We then assessed the respective importance of geographical and environmental predictors in explaining taxonomic and functional compositions. No evidence was found to show that species were more functionally similar than expected in local assemblages. However, the taxonomic turnover components were higher than functional ones (0.471±0.230 vs. 0.243±0.215), with nestedness components showing the opposite pattern (0.063±0.054 vs. 0.269±0.225). This indicated that differences in taxonomic compositions between sites occurred from replacement of functionally similar species. Geographical barriers were the key factor influencing both taxonomic total dissimilarity and turnover components, whereas functional beta-diversity was primarily explained by climatic factors such as minimum temperature of the coldest month. Our findings provide empirical evidence that taxonomic and functional diversity patterns can be independently driven by different ecological processes. Our results point to the importance of clarifying different components of beta-diversity to understand the underlying mechanisms of community assembly. These results also shed light on the assembly rules and ecological processes of terrestrial mammal communities in extreme environments.

INTRODUCTION

Biodiversity is commonly partitioned into alpha, beta, and gamma diversity levels (Whittaker, 1960, 1972), which represent local diversity, spatial variation, and regional diversity, respectively (Legendre & De Cáceres, 2013). These diversities can be measured from different dimensions, i.e., taxonomic, phylogenetic, and functional (Corbelli et al., 2015; Dreiss et al., 2015; Pavoine & Bonsall, 2011; Weinstein et al., 2014). Taxonomic diversity is either measured by traditional species concepts (Dreiss et al., 2015) or recently developed species delimitation approaches (Thormann et al., 2018). Phylogenetic diversity measures genetic relatedness, and thus represents evolutionary history (Cavender-Bares et al., 2009). Functional diversity explores how species adapt and respond to environments by assessing their trait metrics (Díaz et al., 2013).

A central concept in community ecology is the trait-environmental association (Keddy, 1992). Habitat filtering excludes maladapted species with inappropriate trait combinations for given abiotic conditions, which can lead to trait convergence, while competitive interactions prevent similar species from co-existing, i.e., “limiting similarity” (Fukami et al., 2005; Kraft et al., 2015; Mayfield et al., 2010). Such mechanisms have been examined via alpha-diversity measurements (Cornwell et al., 2006; De Arruda Almeida et al., 2019; Fichaux et al., 2019; Mouchet et al., 2010; Mouillot et al., 2007). However, habitat filtering may also work regionally as well as locally (Cadotte & Tucker, 2017; Kraft et al., 2015), and indeed may be more pronounced regionally as competition diminishes with coarsening scale (Mod et al., 2020). In isolated but homogenous habitats, species are expected to be taxonomically divergent but similar in their traits (Penone et al., 2016). Such mechanisms link to the explanation for beta-diversity. Nekola & White (1999) and Soininen et al. (2007) suggested that the distance decay of similarity, i.e., beta-diversity, could originate from three mechanisms, i.e., dispersal limitation, environmental variation, and neutral processes. However, they did not specify how these mechanisms function in different diversity dimensions.

Research in the last few decades has suggested that taxonomic beta-diversity is influenced by the joint effects of dispersal limitation and the environment (Gaston et al., 2007; Qian & Ricklefs, 2012; Wen et al., 2016a, 2016b), although they may have unequal influence. Melo et al. (2009) exclusively discussed the correlations between environments and taxonomic beta-diversity. Species are sorted along environmental gradients by their tolerance (Corbelli et al., 2015). However, a rarely mentioned but commonly adopted perspective is that environmental gradients do not sort species themselves but rather it is their trait combinations (Keddy, 1992; Kraft et al., 2015). Therefore, to understand ecological processes requires the separate examination of taxonomic and functional diversity (Carvalho et al., 2020; Villéger et al., 2012). Evaluating the mechanisms underlying different diversity dimensions is of increasing concern as these dimensions show distinctive patterns when they are differently determined by ecological drivers (e.g., Carvalho et al., 2020; Cilleros et al., 2016; Fluck et al., 2020; Peláez & Pavanelli, 2019). For example, strong habitat filtering can lead to high taxonomic but low functional beta-diversity (Pavoine & Bonsall, 2011; Weinstein et al., 2014). Furthermore, as beta-diversity can be decomposed into turnover and nestedness components (Baselga, 2010; Villéger et al., 2013), differences in specific components may help to better understand diversity patterns and the corresponding ecological processes. High taxonomic nestedness components indicate prevalent post-glacial recolonization (Dobrovolski et al., 2012; Wen et al., 2016b). High taxonomic turnover components imply strong dispersal limitation (Cilleros et al., 2016; Gianuca et al., 2017; Wen et al., 2016b). Low functional turnover components suggest that the regional species pool is functionally redundant (Villéger et al., 2013).

Alpine mountain areas are considered “sky islands” (Itescu, 2019; Spasojevic et al., 2014). Rates of allopatric speciation are especially high in these high elevation islands (He & Jiang, 2014; Knowles, 2001) as the probability of dispersing between mountain tops is generally low (Brown, 1971; Liu et al., 2012). Despite geographical isolation, however, these alpine habitats experience similar environmental conditions. Alpine ecosystems (i.e., high-elevation habitats above the climatic tree line) belong to a distinctive biogeographic unit represented across all continents and latitudes (Testolin et al., 2020). At the global scale, the mean temperature in the growing season at the tree line is ~6.7 °C (Körner & Paulsen, 2004) and vegetation shows a similar transition from dark coniferous forests to treeless scree (Körner, 2012; Testolin et al., 2020). Alpine habitats are generally considered to be harsh (Reymond et al., 2013; Spasojevic et al., 2014) due to their colder temperature, lower biomass, poorer oxygen, and higher UV levels with increasing elevation (Barry, 2008). These features raise the question of whether segregated alpine communities vary taxonomically but exhibit similar trait combinations when facing similar environmental pressures. However, our knowledge on the patterns and determinants of alpine assemblages remains limited and is generally restricted to plant species (e.g., Spasojevic et al., 2014; Spasojevic & Suding, 2012; Stark et al., 2017), with animal-based data, especially for mammals, which play critical roles in ecosystem functioning and biodiversity conservation, still scarce.

Non-volant small mammals (hereafter ‘small mammals’) are defined as small-sized mammals that are unable to fly (Barnett & Dutton, 1995). They are distributed widely across most terrestrial ecosystems and constitute more than 50% of species in the Class Mammalia (Burgin et al., 2018). Small mammals are naturally weak dispersers (Whitmee & Orme, 2013) as they have strong affinity to their habitats (Smith et al., 2008; Wells et al., 2004). More than 180 small mammal species have been recorded in the mountains of south-west China (including Order Scandentia, Carnivora, Eulipotyphla, Lagomorpha, and Rodentia), >10% of which are limited to areas higher than 3 000 m a.s.l. ( Smith et al., 2008). Thus, they are ideal models to study the effects of geographical isolation and habitat filtering on community assemblages. Although diversity patterns of small mammals and the underlying ecological mechanisms have been widely studied along elevational gradients (e.g., Chen et al., 2017; Dreiss et al., 2015; Du et al., 2017; Heaney, 2001; McCain, 2005; Musila et al., 2019; Ramírez-Bautista & Williams, 2019), few studies have exclusively focused on assemblages in isolated alpine habitats (but see Brown, 1971; Frey et al., 2007), which may be related to a distinct ecological process rather than normal elevational distribution patterns.

The present study aimed to evaluate the relative importance of geographical and environmental factors in structuring taxonomic and functional compositions in sky island topography, i.e., Three Parallel Rivers (TPR) area (He & Jiang, 2014). We tested two main hypotheses: 1) dispersal limitation would facilitate species divergence, while habitat filtering would constrain functional dissimilarity (Fukami et al., 2005; Kraft et al., 2015). Thus, isolated assemblages would vary in taxonomic compositions but exhibit similar functional trait combinations, indicated by a slow increasing rate of functional richness compared to species richness (Mouillot et al., 2007), low functional beta-diversity (Pavoine & Bonsall, 2011; Weinstein et al., 2014), and low functional turnover components (Villéger et al., 2013). 2) Spatial and environmental factors would play different roles in shaping site-to-site taxonomic and functional variations. Taxonomic diversity would be correlated with spatial configuration (Gómez-Rodríguez & Baselga, 2018; Wen et al., 2016b), while functional variation would not be explained by geographical predictors (Spasojevic et al., 2014).

MATERIALS AND METHODS

Ethics approval

The field surveys were permitted by the Forestry Departments of Yunnan and Sichuan Provinces. All methods and procedures were performed in accordance with the guidelines and regulations approved by the internal review board of the Kunming Institute of Zoology, Chinese Academy of Sciences (approval ID: SMKX-2012023 & SMKX-20161115-01).

Study area

The TPR is located in north-west Yunnan, bordering Tibet and Sichuan (Sherman et al., 2008), with representative topography of the Hengduan Mountains (Wen et al., 2016b). Three south-north flowing rivers, i.e., the Salween (Nujiang), Mekong (Lancangjiang), and Yangtze (Jinshajiang), have carved deep parallel canyons and divided the area into four major mountains, i.e., the Gaoligong, Nushan, Yunling, and Shaluli (Sherman et al., 2008; Wang, 2014; Zhang, 1997). It has a minimum east-west distance of 65 km at N27°30′. Elevations of the mountain summits slowly rise from 4 100 m a.s.l. in the south to 6 700 m a.s.l. in the north (Figure 1A). The annual tree line temperature is ~2.5 °C (Wang, 2014). Attributed to the longitudinal range-gorge barrier of the south-west monsoon and accompanying moisture, there is a moderate west-east decline in annual precipitation (Wang et al., 2013; Zhang, 1997). In the high elevation zones, vegetation shows similar patterns, with a rapid transition from dark coniferous forests and dwarf shrubs to meadows and scree (Figures 1B, 2; Sherman et al., 2008).

Figure 1

Study area and field work information

Figure 2

Typical alpine habitats in TPR area

A: Topography and major river systems in Three Parallel Rivers (TPR) area. Triangles mark localities of 18 alpine sites. Respective sites of each major mountain are represented in different colors. Black horizontal lines refer to division of latitudinal groups. B: Typical alpine vegetation transition in our study area (photo from May 2017 at Cangshan Mountain, 3 800 m a.s.l.). Photo by Wen-Yu Song. C: Field sampling protocols of alpine elevation gradients, red bands indicate transects.

A: Langdu Mountain at 4 400 m a.s.l. (photo from October 2017). B: Baima Snow Mountain at 4 400 m a.s.l. (photo from June 2017). C: Shika Snow Mountain at 4 200 m a.s.l. (photo from October 2017). D: Meili Snow Mountain at 4 200 m a.s.l. (photo from June 2017). E: Biluo Snow Mountain at 4 000 m a.s.l. (photo from May 2017). F: Gawagapu at 3 800 m a.s.l. (photo from September 2018). Photos by Wen-Yu Song.

Small mammal sampling

The tree line is considered the best-defined biome boundary in alpine studies (Körner, 2012; Körner & Paulsen, 2004; Testolin et al., 2020). In the TPR, the tree line ranges from 3 800 to 4 300 m a.s.l. in different locations (Sherman et al., 2008; Wang, 2014). Therefore, we restrained the targeted alpine zone to elevations ranging from 200 m above and below the tree line. Specifically, three transects were established, i.e., 200 m below, along, and 200 m above the tree line. At six sites, the summits were at the tree line level or under 3 800 m a.s.l., therefore only two transects were established, i.e., summit and 200 m below (Figure 1C). The transects at different sites were grouped into lower, middle, and higher based on their relative elevational positions. For sites with only two transects, the higher elevational transect was considered to be missing. Thus, we sampled small mammals along 48 transects in 18 prominent alpine sites (Table S1), representing the upper distribution limits of small mammals in the study area. Fieldwork was conducted in September 2013 in Yulong and from December 2016 to November 2018 at the other 17 sites.

We used ~1 000 trap nights per transect to obtain approximate sampling effort. To reach this number, sub-transects were set along each main transect. Three types of sampling tools were used to trap small mammals: i.e., 1) 7.62 cm×8.89 cm×22.86 cm Sherman traps (HB Sherman Traps Inc., Tallahassee, FL, USA); 2) 10 cm×6 cm museum snaps; and, 3) Φ20 cm×20 cm buckets. We defined one Sherman trap and one snap, set 1–2 m apart, as a “trap station” (Chen et al., 2017). Thirty trap stations were placed every 10–15 m, and 10 buckets were buried on potential small mammal trails along each sub-transect. Sherman traps were baited with sugar-free oatmeal, snaps were baited with fresh peanuts, and buckets were unbaited. The capture equipment was first positioned, then adjusted and re-baited the next morning, and moved to another sub-transect after two nights. After deducting missing nights, a total of 48 751 trapping nights, averaging 1 015.6 (min: 940, max: 1 152) for each transect, were set during the field surveys. For every small mammal captured, body weight (BW), head-body length (HB), tail length (TL), hindfoot length (HF), and ear length (EL) measurements were recorded in the field following Pan et al. (2007) for taxonomic identification and functional traits. At least one individual of each morphologically discriminant taxon from each site was prepared as a stuffed voucher specimen, with skulls cleaned for taxonomic identification. All collections were deposited in the Kunming Institute of Zoology, Chinese Academy of Sciences.

Organization of community and trait data

The community data were initially built at the transect scale. As the between-transect distance within a site was neglectable compared with the between-site distance, we considered each site as a sky island. Therefore, the transect-scale community data were combined to the site scale. As such, species in each site were considered as one community. To ensure adequate sampling effort in each site, we first employed coverage-based rarefaction (Chao et al., 2014). Sampling completeness (or “coverage”) is defined as the relative abundance of the observed species represented in the sample (Chao et al., 2014). Here, we considered a site to be sufficiently sampled if sampling completeness was ≥90% (Moreno & Halffter, 2000; Ramírez-Bautista & Williams, 2019). Rarefaction curves and species coverage were calculated using the iNext function in the R package ‘iNext’ (Hsieh et al., 2019). To measure alpha- and beta-diversities based on incidences of species and traits, we transformed the abundance-based community data into presence-absence data.

We used both morphological and behavioral features as functional traits. Five body measurements (i.e., BW, HB, TL, HF, and EL) and the tail-body ratio (TB: proportion of TL to HB) were obtained in the field. Diet preference, activity cycles, and life mode were obtained from the combined MammalDIET and MammalDIET 2 datasets (Gainsbury et al., 2018; Kissling et al., 2014) and other published literature (Pan et al., 2007; Smith et al., 2008). These morphological and behavioral features can comprehensively address the adaptive strategies of such animals to alpine habitats. From a heat retention perspective, Bergmann’s rule states that body size (represented by BW and HB) tends to be larger in colder environments (Bergmann, 1847), whereas Allen’s rule argues that extremities (represented by TL, HF, and EL) tend to be smaller (Allen, 1877). In addition, HB and HF are considered relevant to mobility (Forsman et al., 2011; Whitmee & Orme, 2013), whereas TL and TB are indicative of arboreality (Alroy, 2019; Du et al., 2017). Diet preference, activity cycle, and life mode essentially impact resource acquisition and usage – i.e., what, when, and where (Dreiss et al., 2015). We listed continuous variables by the mean value of individuals of each species, with behavioral characteristics then categorized (Table 1, Supplementary Table S3).

Table 1

List and respective formats of morphological and behavioral traits used for functional alpha- and beta-diversity estimates

Type	Functional component	Attribute	Value
Mensural	Morphological	Body weight	Mean (g)
		Head-body length	Mean (mm)
		Tail length	Mean (mm)
		Hindfoot length	Mean (mm)
		Ear length	Mean (mm)
		Tail-body ratio	Proportion (%)
Categorical	Diet	Herbivore
		Omnivore
		Carnivore
	Activity cycle	Diurnal
		Nocturnal
		Both
	Life mode	Fossorial
		Terrestrial
		Semi-fossorial
		Arboreal
		Semi-aquatic

Quantification of different levels of taxonomic and functional diversity

Alpha-diversity

Taxonomic alpha-diversity was measured as observed species richness and functional alpha-diversity was measured as functional richness (FRic; Villéger et al., 2008). FRic represents the amount of functional space occupied by an assemblage (Mouchet et al., 2010) and is calculated as the minimum convex hull volume of all species in a site (Cornwell et al., 2006). It can either remain unchanged or increase as species richness increases (Mason et al., 2005) and reflects assembly rules by different accumulation rates compared to species richness (Mouchet et al., 2010). Therefore, it is frequently used to distinguish different ecological mechanisms (e.g., De Arruda Almeida et al., 2019; Fichaux et al., 2019). FRic was calculated with the dbFD function in the R package ‘FD’ (Laliberté et al., 2014). Specifically, a Gower’s pairwise distance matrix was computed, after which principal coordinate analysis (PCoA) was applied to the matrix. The first four PCoA axes, which explained 59.5% of trait variations, were then retained to calculate the FRic (Villéger et al., 2014). Species richness and FRic were calculated for both transect-scale (n=48) and site-scale (n=18) community data. FRic results were standardized between 0 and 1 (Villéger et al., 2008). We also calculated functional dispersion (FDis; Laliberté & Legendre 2010), measured as the mean distance of each species to the centroid of all species in the community. However, the distribution of FDis was random and uncorrelated with other metrics (Supplementary Text S1, Figures S2, S3, Table S6), therefore we only reported the FRic results here. In addition, as the alpha-diversity metrics were not normally distributed (Shapiro-Wilk normality test), we used Spearman’s coefficients to evaluate the correlation between site-scale taxonomic and functional alpha-diversity (log-transformed).

Beta-diversities and their components

The pairwise indices of taxonomic and functional beta-diversity decomposition were calculated based on Sørensen dissimilarity coefficients (Baselga, 2010; Villéger et al., 2013). For taxonomic beta-diversity, the first dissimilarity index refers the total dissimilarity among assemblages (TDsor), the turnover component (TDsim) refers to compositional changes due to species replacement, and the nestedness component (TDsne) refers for the dissimilarity caused by sites with lower species richness being subsets of richer sites (Baselga, 2010). Functional beta-diversity decomposition is based on the volume of convex hull intersections in a multi-dimensional functional space (Villéger et al., 2013). Following FRic, the first four PCoA axes (representing 59.5% of trait variations) were used to compute the three indices representing functional total dissimilarity (FDsor), functional turnover component (FDsim), and functional nestedness component (FDsne). We used the means of pairwise indices to evaluate the total dissimilarities and respective components (Fichaux et al., 2019; Si, 2016). Beta-diversity decomposition was carried out with the beta.pair and functional.beta.pair functions in the R package ‘betapart’ (Baselga et al., 2018). Mantel tests (999 permutations) were then applied to evaluate the correlations among taxonomic and functional total dissimilarity, turnover, and nestedness components. Mantel tests were performed with the mantel function in the R package ‘vegan’ (Oksanen et al., 2019) based on Spearman’s coefficients.

Null models

Our null hypothesis was that species are able to disperse freely among the sky islands (no dispersal limitation). We randomized the site-scale community matrix with the “independent swap” algorithm (Gotelli, 2000). This constrained randomization maintains species richness and occurrence frequency, but not spatial contagion or dispersal limitation (Swenson, 2014). For local assembly (alpha-diversity), FRic was calculated with 1 000 randomizations. For regional assembly (beta-diversity), we calculated functional total dissimilarity, turnover, and nestedness components with 1 000 randomizations. All null models were established using site-scale community data and the first four PCoA axes of Gower’s distance of species traits (same as alpha- and beta-diversity quantifications). The null values were then compared with observed values by the standardized effect size (SES), calculated as: SES=(observed-mean(null))/SD(null). For FRic, communities with negative SES are considered functionally convergent, which points to the effects of habitat filtering; whereas those with positive SES are considered functionally complementary, which points to the effects of competition (Fichaux et al., 2019; Mouchet et al., 2010). For functional beta-diversity and its components, positive/negative SES suggests that the observed pairwise site differences are higher/lower than those species distributions not deterministically structured (Villéger et al., 2013). The significance (P) of SES values was calculated following Swenson (2014) and two-tailed tests were used to determine whether the observed values rejected the null hypothesis. We considered the observed values to be significantly lower or higher than null expectations when P<0.025 orP>0.975.

Contributions of spatial and environmental factors in structuring beta-diversity

Spatial and environmental variables

To test the effects of geographical barriers, i.e., the main rivers, we grouped the studied sites into the four major mountains (Zhang, 1997) and coded them west-east directionally (i.e., Gaoligong=1, Nushan=2, Yunling=3, Shaluli=4). As biodiversity patterns are widely consistent with latitudinal gradients (Lamanna et al., 2014; Willig & Presley, 2013), we grouped the study sites into 0.5° north-south bands. We then coded the latitudinal groups from 1 (southernmost) to 5 (northernmost). We did not use spatial distance because geographical groups better describe the unique topography in the study area (Supplementary Table S7). We used temperature, precipitation, potential evapotranspiration (PET), and the normalized difference vegetation index (NDVI) as environmental variables, which are commonly applied to explain diversity patterns in vertebrates (Melo et al., 2009; Qian & Ricklefs, 2012; Wen et al., 2016a). Elevation was excluded as an explanatory variable as it has a limited influence on habitat variations and diversity patterns in our case (Supplementary Figures S4–S6). We downloaded the 19 bioclimatic variables containing temperature and precipitation information from the Worldclim dataset 2.0 (http://worldclim.org/version2, 30 s spatial resolution). Annual PET data were downloaded from https://cgiarcsi.community. The NDVI data were downloaded from https://earthdata.nasa.gov (Product: MOD13A1) at 500 m spatial resolution and averaged for 16 day layers for 2018. We extracted environmental values from the coordinates of start points for each sub-transect and acquired the mean values of each site. We standardized environmental variables to a zero mean and unit variance, then applied principle component analysis (PCA) to reduce dimensionality (Carvalho et al., 2020; Qian & Ricklefs, 2012). The first four principal components, which represented 94.6% of environmental variation, were retained for further analyses.

Explanatory power of different factors for beta-diversities

We first conducted PCoA for the six dissimilarity matrices (i.e., two dimensions: taxonomic and functional, each including three indices: total dissimilarity, turnover, and nestedness). We then adopted distance-based redundancy analysis (db-RDA) with site scores of the PCoA results and obtained adjusted R2 values for all candidate variables. The adjusted R2 values were used as the upper scope in the following model selections. Stepwise model selection was performed using the “forward” method with 999 permutations and maximum adjusted R2 at every step; the initial explanatory variables were constant (i.e., 1). Selection stopped when the (1) adjusted R2 started to decrease, (2) adjusted R2 of the scope was exceeded, or (3) selected permutation P-value exceeded the alpha significance of 0.05 (Blanchet et al., 2008). The variables selected by the best models were used as explanatory variables and the pairwise matrices for corresponding beta-diversity components were used as response variables in PERMANOVA (omitted if no explanatory variables were selected). Significance tests for PERMANOVA were based on 999 Monte Carlo randomizations. The PERMANOVA R2 values were used as estimates of the magnitude of the effect of the explanatory factors on the response variables. The ordiR2step and adonis functions in the R package ‘vegan’ (Oksanen et al., 2019) were used to conduct model selection and PERMANOVA, respectively.

RESULTS

A total of 5 686 small mammals, comprised of Rodentia (70.75%), Eulipotyphla (26.17%), and Lagomorpha (3.08%), were captured, averaging a capture success rate of 11.74%. Eothenomys custos (19.66% of total collection) was the most frequently captured species, followed by Apodemus ilex (12.59%), Sorex bedfordiae (11.68%), and Eothenomys cachinus (7.95%). Ochotona thibetana (2.41%) was most abundant in Lagomorpha. Total observed species richness was 44 (Supplementary Table S2), ranging from 8 to 19 per site. Coverage-based rarefactions showed that species coverage was higher than 90% in all studied sites, indicating that the communities were adequately sampled (Supplementary Table S1 and Figure S1).

Correlations and comparisons of different levels and dimensions of diversity

When alpha-diversity was evaluated at the transect scale, species richness and FRic tended to decrease towards higher, treeless habitats (Figure 3). Average species richness was 11.06±2.92 at lower, 9.12±1.95 at middle, and 7.58±1.56 at higher elevational positions. Average FRic values were 0.23±0.15, 0.16±0.12, and 0.08±0.07 from lower to higher elevations, respectively. The differences were significant (for species richness, ANOVA P=0.0007, df=2; for FRic, ANOVA P=0.0103, df=2).

Figure 3

Vertical distributions of species richness (A) and functional richness (FRic) (B)

Transects were grouped depending on their relative elevational positions regarding the tree line. Horizontal solid lines inside boxes are medians. Notches are 95% confidence intervals of medians. Box heights represent value ranges from 25% to 75% (middle 50% of scores) for groups. Upper and lower whiskers represent scores outside middle 50%. Black dots mark outliers.

When alpha-diversity was evaluated at the site scale, average species richness was 13.22±2.82 and FRic was 0.33±0.17. These two indices were positively correlated (r=0.856, P<0.001;Figure 4). Mantel tests showed significant positive correlations for the total dissimilarity, turnover, and nestedness components of taxonomic and functional beta-diversities. However, the coefficients were lower than that for alpha-diversity (total dissimilarity: r=0.631, P=0.001; turnover: r=0.367, P=0.001; nestedness: r=0.586, P=0.001; Figure 4). Taxonomic and functional total dissimilarity had similar average values (TDsor=0.535±0.208 and FDsor=0.512±0.230). However, the components in taxonomic beta-diversity showed clear dominance of turnover components (TDsim=0.471±0.230 and TDsne=0.063±0.054), whereas nestedness components were higher than turnover components in functional beta-diversity (FDsim=0.243±0.215 and FDsne=0.269±0.225; Figure 5).

Figure 4

Relationship between taxonomic and functional alpha-diversity (A), beta-diversity (B), turnover(C), and nestedness (D) components of small mammal assemblages from 18 alpine sites in TPR

Figure 5

Comparison of taxonomic (A) and functional (B) pairwise (n=153) Sorensen dissimilarity, turnover, and nestedness components of small mammal assemblages from 18 alpine sites in TPR

For alpha-diversity, functional richness (FRic) and species richness were log-transferred. Gray areas indicate 95% confidence intervals.

Horizontal solid lines inside boxes are medians. Notches are 95% confidence intervals of medians. Box heights represent value ranges from 25% to 75% (middle 50% of scores) for groups. Upper and lower whiskers represent scores outside middle 50%. Black dots mark outliers.

For alpha-diversity, the observed FRic values in seven alpine communities were lower (SES<0) than the null expectations, though none were significant; in contrast, values in 11 sites were higher (SES>0) than the null expectations, one of which was significant (Supplementary Table S4). For beta-diversity, among the 153 pairs of site-to-site indices, nine, ten, and two observed values deviated significantly from the null distributions for functional total dissimilarity, turnover components, and nestedness components, respectively (Supplementary Table S5).

Explanatory power of spatial and environmental factors on beta-diversity

Based on the PCA results of environmental variables, the first (env1) and second (env2) axes explained 69.76% and 13.24% of site-to-site environmental variations, respectively. The 10 most important variables in env1 and env2 are shown in Figure 6. The five most important variables in env1 were Bio06 (min temperature of coldest month), Bio09 (mean temperature of driest quarter), Bio01 (annual mean temperature), Bio12 (annual precipitation), and Bio11 (mean temperature of coldest quarter).

Figure 6

Contributions of environmental variables to first and second axes of PCA results

First 10 variables are shown for each axis. Variables with a contribution above the red dashed line were considered important contributors to the axis.

For taxonomic total dissimilarity and turnover components, only “mountain range” was chosen by the forward model selections (adjusted R2=0.250 and 0.211). Based on PERMANOVA, the contributions of mountain range to total taxonomic dissimilarity and turnover components were significant (Table 2). Model selection showed that no variable had a significant effect on taxonomic nestedness components (adjusted R2=0). Latitudinal gradients and all environmental factors were excluded by model selection as explanatory variables for taxonomic beta-diversity and its components.

Table 2

Explanatory factors selected by best models and PERMANOVA results

	Explanatory factor	Pseudo-F	R2
Factors were added sequentially following “forward” stepwise model selection results using permutation tests (999). R2 values based on Monte Carlo randomizations (999) for each selected factor are presented (*: P<0.05;**: P<0.01;***: P < 0.001). No candidate variable was selected to explain taxonomic and functional nestedness components.
Taxonomic beta-diversity
TDsor	Mountain range	11.434	0.417***
	Residual		0.583
TDsim	Mountain range	16.779	0.512***
	Residual		0.488
TDsne	None
Functional beta-diversity
FDsor	env1	4.092	0.182**
	env2	3.347	0.149**
	Residual		0.668
FDsim	env1	26.616	0.550**
	Mountain range	6.756	0.140*
	Residual		0.310
FDsne	None

For functional total dissimilarity, the first and second PCA axes of environmental variables were selected when spatial variables were excluded by forward model selection (adjusted R2=0.297). PERMANOVA was performed to show significant contributions (Table 2). For functional turnover components, model selection retained env1 and mountain range (adjusted R2=0.05), with their contributions found to be significant (Table 2). No variables were selected to explain functional nestedness components (adjusted R2=0). Latitudinal gradient was excluded by model selection as an explanatory variable for functional beta-diversity and its components.

DISCUSSION

Taxonomic and functional diversity patterns: correlations and distinctions

At the transect scale, results showed that both taxonomic and functional richness tended to decrease with higher elevation. Diversity patterns of small mammals usually monotonically decrease at high elevations, despite various patterns occurring at low-mid elevations (McCain, 2005; see also Chen et al., 2017; Heaney, 2001; Musila et al., 2019). Although some studies have not found this decreasing pattern (Dreiss et al., 2015; Ramírez-Bautista & Williams 2019), their “high elevation” levels are not “true” alpine as their summits did not reach the forest edge. Near the tree line (Körner, 2012; Körner & Paulsen, 2004), environmental conditions become harsher with increasing elevation due to the decline in temperature, oxygen, and productivity, and increase in sun exposure (Barry, 2008). This means that available niche spaces for small mammals in alpine habitats are considerably narrow. Food resources for small mammals also become scarcer when productivity is low (Heaney, 2001; Ramírez-Bautista & Williams, 2019). Moreover, meadow and scree habitats may be uninhabitable for species with certain strategies and requirements. For example, arboreal species (e.g., squirrels) are unlikely to survive in treeless habitats (Wells et al., 2004). Analogously, increased rock cover is not suitable for fossorial species (e.g., moles) that rely on humus and worms (Gainsbury et al., 2018; Pan et al., 2007; Smith et al., 2008). Furthermore, typical or unique species traits, e.g., arboreal species with long tails (Alroy, 2019) and high tail-body ratio (Du et al., 2017) and fossorial species without ears (Pan et al., 2007), are less likely to appear.

Against our hypothesis, the alpha-diversity results suggested that taxonomic and functional diversity had similar patterns, i.e., sites with higher species richness also exhibited higher functional richness. These results are in contrast with previous studies, which reported a plateau in functional richness with increased species richness when species are functionally redundant (Oliveira et al., 2016; Villéger et al., 2013). Our null models provided further insight into this pattern. Specifically, the observed FRic values were not significantly lower than the null expectations. Therefore, no signs of habitat filtering were detected, at least within each site. These results do not conform to the concept that environmental stress, e.g., cold temperatures (Körner & Paulsen, 2004; Reymond et al., 2013), will filter species with inappropriate traits and thereby increase functional similarity (Cornwell et al., 2006; Denelle et al., 2019; Fukami et al., 2005; Reymond et al., 2013). On the contrary, many study sites showed positive SES values for FRic, one of which was significant. These results may be explained by harshness, i.e., stress, driving competition and limiting functional similarity in co-existing species (Cavender-Bares et al., 2009; Mayfield et al., 2010; Spasojevic & Suding, 2012; Villalobos et al., 2017). These findings indicate that local assembly mechanisms in alpine habitats are sophisticated, with habitat filtering and competition occurring simultaneously under environmental pressures (Helmus et al., 2007). For example, recent study found that high-elevation pikas exhibit convergence in resource-use traits but striking divergence in craniums, indicating multi-level adaptations to extreme environments (Feijó et al., 2020). As a result, assembly mechanisms may be blurred (Helmus et al., 2007).

Unlike our expectation, taxonomic and functional total dissimilarity demonstrated similar levels. However, the correlations of beta-diversity components, especially turnover components, between the two dimensions were lower than that for alpha-diversity. The taxonomic turnover components were approximately two times that of the functional components, whereas the taxonomic nestedness components were approximately a quarter of that of the functional components. These results indicate that many sites acted as functional subsets to the others, despite different species compositions. These findings also provide new implications for regional functional convergence (Carvalho et al., 2020; de Bello et al., 2009; Rodríguez & Ojeda, 2014). Traditionally, species in isolated but homogenous habitats are expected to have high taxonomic but low functional beta-diversity (Pavoine & Bonsall, 2011; Penone et al., 2016; Weinstein et al., 2014). However, our results indicated that regional functional convergence can also be revealed by high taxonomic but low functional turnover components. In other words, differences in taxonomic compositions between sites are increased by replacement with functionally similar species (Si, 2016; Villéger et al., 2013). The observed increase in levels of functional nestedness components may be due to different niche filtering intensity among communities (Villéger et al., 2013). This regional functional convergence was supported by the null models for functional beta-diversity and its components, with only a few pairs of sites departing from the null distributions. These findings highlight the need to partition different beta-diversity components. This is because, when different dimensions of beta-diversity have similar levels, traditional measurements may overlook the differences that arise in their respective components.

Underlying drivers for taxonomic and functional beta-diversity patterns: relative roles of geographical and environmental factors

Results showed contrasting drivers of taxonomic and functional beta-diversity, thus supporting our hypothesis about the roles of spatial and environmental factors in determining diversity dimensions. Considering functional traits represent adaptive strategies of species to survive in given habitat conditions (Díaz et al., 2013; Kraft et al., 2015), theoretically, functional beta-diversity can be explained by environmental factors independently from dispersal (Cadotte & Tucker, 2017; Díaz et al., 2007). However, it is difficult to examine this idea directly. For example, previous studies have reported on the joint effects of dispersal and niche processes in functional beta-diversity (Fluck et al., 2020; Peláez & Pavanelli, 2019; Siefert et al., 2013). This is due to the simple fact that, in spatial terms, environmental variations cannot exist independently. However, if the spatial effects are maximized (i.e., geographical segregation) and environmental influences are minimized (i.e., homogenous habitats), taxonomic beta-diversity should be explained by geographical factors only, rather than environmental factors. In contrast, functional beta-diversity should not be explained by geographical factors. Spasojevic et al. (2014) provided relevant evidence after studying tundra communities in multiple sky islands. They found that tundra communities in colder high elevation sites were functionally similar to each other, and thus geographic isolation did not influence functional beta-diversity. In the present study, we extended these patterns to animals and found that taxonomic beta-diversity was solely explained by geographical barriers. Thus, these results suggest that taxonomic and functional compositions can be independently driven by different ecological processes.

When looking at the specific patterns, there is a rationale behind the relationship between taxonomic beta-diversity and geographical barriers. Across global terrestrial topography, species turnover is especially high in mountain areas, which is usually attributed to strong dispersal limitations (Melo et al., 2009; Qian & Ricklefs, 2012; Wen et al., 2016b). Unsuitable habitats, such as lowlands, canyons, or rivers, are well-known dispersal barriers for alpine species (Brown, 1971; Itescu, 2019). In the TPR, the main rivers running longitudinally are physical barriers for small alpine mammals (He et al., 2019; He & Jiang, 2014; Liu et al., 2012). Due to limited dispersal corridors, allopatric speciation is promoted in sky islands (He & Jiang, 2014; Knowles, 2001). This process can yield a high proportion of turnover components in taxonomic beta-diversity (Cilleros et al., 2016; Gianuca et al., 2017; Wen et al., 2016b), as observed in the present study. Therefore, our results indicate that, when environmental variations are minimized, taxonomic beta-diversity in the TPR sky islands mainly originates from allopatric speciation. We note that the taxonomic nestedness component was not explained. The glacial extinction and post-glacial recolonization likely caused prominent nested patterns in taxonomic compositions (Baselga, 2010; Dobrovolski et al., 2012; Wen et al., 2016b). The absence of latitudinal patterns in the taxonomic nestedness component may be due to the topography of the TPR providing a refugium for cold-adapted species (He et al., 2016; Qu et al., 2011). As a result, the post-glacial recolonization in the alpine area mainly occurred altitudinally (He & Jiang, 2014; Wen et al., 2016b).

In the present study, we made no assumptions about which ecological factors may explain functional beta-diversity and its components. This is because the small mammals were sampled along a unitary biome boundary, i.e., the tree line (Körner 2012; Testolin et al., 2020; Wang, 2014). However, we still found that the functional total dissimilarity and turnover components were primarily explained by the first environmental axis. This axis generally explains joint variation of temperature and precipitation. Notably, the coldest scenario (represented by Bio06, Bio09, and Bio11) accounted for three of the five most important variables within it. This result supports that high-latitude/altitude range limits are formed by stressful abiotic environments (e.g., cold climates; Darwin, 1859; Normand et al., 2009) from a functional perspective. The functional turnover component was primarily explained by the first environmental axis and marginally by geographical barriers. This may be because environments still slightly co-vary with spatial structures in the study area (Sherman et al., 2008; Wang et al., 2013; Zhang, 1997). Finally, functional nestedness was not explained by geographical or environmental factors. We assume that the gain/loss of functional space may be attributed to undocumented factors, such as microtopography and microclimate (Lampei et al., 2019; Opedal et al., 2015; Stark et al., 2017). For example, some sites contained more cliffs or rocks, and were thus less favored by species that live in tunnels. Such information is not represented in the available data.

CONCLUSIONS

Due to their wide distribution but limited dispersal ability, diversity patterns of small mammals have been studied extensively; however, communities in alpine habitats have received less attention. We used field data from broad taxa to study taxonomic and functional diversity at site and regional scales. Results showed that habitat filtering was not prominent in local assembly. Thus, the viewpoint that species evolve similar traits to cope with environmental stress may overly simplify local assembly processes in extreme habitats. However, habitat filtering was detected when communities were examined between distant sky islands. We propose that partitioning beta-diversity into different components could provide further insight into vague patterns. We found contrasting ecological drivers for taxonomic and functional compositions. This finding highlights the importance of studying different diversity dimensions separately. We propose that further study should incorporate phylogenetic data as such information can bridge species and trait concepts (Cavender-Bares et al., 2009; Corbelli et al., 2015; Kraft et al., 2015).

SCIENTIFIC FIELD SURVEY PERMISSION INFORMATION

Permission for field surveys in Yunnan Province was granted by the Yunnan Provincial Forestry and Grassland Bureau, No. [2017]084 & No. [2018]40.

Supplementary data to this article can be found online.

Click here for additional data file.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORS’ CONTRIBUTIONS

W.Y.S., X.Y.L., and X.L.J. conceived and designed the study. W.Y.S., Z.Z.C., and Q.L. performed the fieldwork. W.Y.S., K.O.O., S.W.H., and X.L.J. carried out species identification. W.Y.S. and X.Y.L. conducted analyses and wrote the manuscript. All authors read and approved the final version of the manuscript.