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).
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).
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.
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).
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
Table 1. List and respective formats of morphological and behavioral traits used for functional alpha- and beta-diversity estimates
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, Figure 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).
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.
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 or P>0.975.
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.
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.
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, Figure S1).
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).
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, 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).
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).
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.
Explanatory factor Pseudo-F R2 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 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.
Table 2. Explanatory factors selected by best models and PERMANOVA results
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.
Isolated alpine habitats reveal disparate ecological drivers of taxonomic and functional beta-diversity of small mammal assemblages
- Received Date: 2020-04-16
- Accepted Date: 2020-08-18
- Available Online: 2020-08-30
- Beta-diversity partitioning /
- Community assembly /
- Environmental stress /
- Habitat homogeneity /
- Hengduan Mountains /
- River barriers /
- Sky islands /
- Tree line
Abstract: 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.
|Citation:||Wen-Yu Song, Xue-You Li, Zhong-Zheng Chen, Quan Li, Kenneth Otieno Onditi, Shui-Wang He, Xue-Long Jiang. Isolated alpine habitats reveal disparate ecological drivers of taxonomic and functional beta-diversity of small mammal assemblages. Zoological Research. doi: 10.24272/j.issn.2095-8137.2020.085|