Birds in the Himalayas: What drives beta diversity patterns along an elevational gradient?

Abstract Beta diversity patterns along elevational gradients have become a hot topic in the study of biogeography and can help illuminate the processes structuring mountain ecosystems. Although elevational species richness patterns have been well documented, there remains much uncertainty over the causes of beta diversity patterns across elevational gradients. We conducted bird surveys and obtained high‐resolution climatic data along an elevational gradient in Gyirong Valley in the central Himalayas, China, between 1,800 and 5,400 m elevation. In total, we recorded 182 bird species (including 169 breeding birds). We simulated beta diversity patterns with the mid‐domain effect (MDE) null model and conducted distance‐based redundancy analyses (db‐RDA) to relate beta diversity to dispersal limitations, spatial constraints, habitat complexity, contemporary climate, and historical climate. Mantel tests and variation partitioning were employed to identify the magnitude of independent statistical associations of environmental factors with beta diversity. Patterns of empirical and simulated beta diversity were both hump‐shaped, peaking at intermediate elevations. The db‐RDA indicated that beta diversity was correlated with changes in spatially structured environmental factors, especially with contemporary climate and habitat complexity. Mantel tests and variation partitioning also suggested that climate dissimilarity was the major independent correlate of beta diversity. The random community structure and spatial constraints may also contribute to the overall hump‐shaped pattern. Beta diversity of bird communities in Gyirong Valley could be explained by the combination of different factors but is mainly shaped by the spatially structured environmental factors, especially contemporary climate.


| INTRODUC TI ON
Uncovering the mechanisms responsible for variation in biodiversity across space and time remains one of the central aims of ecology.
Although elevational species richness patterns have been well documented (Rahbek, 1995), there remains much uncertainty over the patterns of beta diversity across elevational gradients (McCain & Beck, 2016). McCain and Beck (2016) argued that elevational beta diversity patterns were changeable and difficult to use as a general explanation of elevational richness patterns. However, the highest beta diversity frequently appears at intermediate elevations, commonly accompanied by high biodiversity or associated with ecotones (Clements, 1916;Levanoni, Levin, Pe'er, Turbe, & Kark, 2011;Mena & Vazquez-Dominguez, 2005;Whittaker, 1975).
Stochastic processes (such as random community structure) are also considered to influence patterns of beta diversity (Chase, 2010;Chisholm & Pacala, 2011). Inspired by studies of elevational richness patterns, mid-domain effect null models (MDE; Colwell & Hurtt, 1994) have been used to test whether species ranges vary individualistically or whether species cluster into structured zonal communities (Herzog, Kessler, & Bach, 2005;McCain & Beck, 2016;Mena & Vazquez-Dominguez, 2005). Moreover, it has been argued that an interacting framework of different drivers rather than any single one should explain the beta diversity patterns (Baselga & Jimenez-Valverde, 2007;Qian & Ricklefs, 2007).
The Himalaya Mountains, the greatest mountain range in the world, offer a unique environment to conduct studies on the mechanisms of beta diversity patterns and the ecological theories of species distribution. However, numerous studies in this area mainly focused on the elevational patterns of species richness (Acharya, Vetaas, & Birks, 2011;Grytnes & Vetaas, 2002;Hu et al., 2017Hu et al., , 2018Hu et al., , 2014Joshi & Bhatt, 2015;Pan et al., 2016). Studies on beta diversity in the Himalayas have focused on plants in northwest Himalayan region (de Bello, Dolezal, Ricotta, & Klimesova, 2011;Saha, Rajwar, & Kumar, 2016) and trans-Himalayan region (Paudel & Vetaas, 2014).
To identify comprehensive beta diversity patterns, more rigorous studies in different regions for different taxa are needed.
Identifying underlying drivers of beta diversity can help to predict how species ranges or community composition shift with abiotic factors and biotic interactions (Jankowski et al., 2013;Kissling, Field, Korntheuer, Heyde, & Bohning-Gaese, 2010). To better understand the mechanisms of beta diversity and help identify the necessary management actions toward preservation of biodiversity in the Himalayas, we examine beta diversity patterns of birds in the Himalayas by sampling continuously across an elevational gradient and using high-resolution environmental or spatial variables to analyze diversity patterns.
We test the hypotheses that (a) beta diversity is higher at intermediate elevations and related to species richness; (b) random community structure contributes to the beta diversity patterns; and (c) beta diversity is determined mainly by environmental conditions.

| Study area
Gyirong Valley (28°15′-29°0′N, 85°6′-85°41′E) is located in the central Himalayas, China, and spans more than 4,000 m in elevation, from Resuo Village (1,700 m above sea level, asl) to the summit of Mt. Mala (5,770 m asl). The distance from the bottom of the valley to the summit of Mt. Mala is 79 km. Gyirong Valley is in a transition zone between the Oriental and Palearctic realms (Hu et al., 2014;Zhang, 2011 (Hu et al., 2017). The wet season approximately lasts from May to October, whereas the dry season is from November to April. The vegetation outside Gyirong Valley (alpine tundra with sparse grasses) is different from the vegetation within the valley (evergreen broadleaf forest, broadleaf mixed forest, coniferous forest, shrub and grasses). We used topography to define the study region, with the highest ridge lines and the major watercourses surrounding Gyirong Valley as the boundaries (Figure 1).

| Sampling and sampling effort
Bird censuses were conducted from 1,800 to 5,400 m asl along the elevational gradient in Gyirong Valley. Field sampling could not be done at lower or higher elevations because the lowest elevation in this area is 1,800 m asl and the habitats above 5,400 m asl are inaccessible cliffs and glaciers. Four of the five climate zones were sampled, excluding the Alpine Frigid Zone. We divided the elevational gradient into twelve consecutive 300-m-wide elevational bands. There were three 300-m elevational bands in each climate zone, and we placed three transect lines in each 300-m elevational band to cover most of the available habitat types (for a total of 36 transect lines, see Figure 1). We standardized sampling effort across the gradient by limiting the total length of the three transect lines in each 300-m elevational band to 7.5 km, to reduce sampling bias (Rahbek, 2005); the length of each transect was between 2 and 3 km. Field surveys were carried out between 20 min after dawn and 10:00 hours and between 16:00 hours and 20 min before sunset (local time) by the same proficient observers (J. J. Li and H. F. Cao). Bird species within 50 m of the observers were recorded during each census. The movement rate of the observers depended on bird activity. To reduce the temporal autocorrelation among transects, the transects in the different elevational bands were sampled in a random order. Replicated bird censuses were made four times during the wet season, from May to June 2012, August 2012, from September to October 2012, and from July to August 2013. We used Zheng (2011) as the taxonomic system of birds in this study. Only breeding birds (including residents and summer migrants) were used for statistical analyses, as migratory birds could cause potential bias (Pan et al., 2016).
Because of the strong movement ability of birds and the relatively small study area, we interpolated the presence of species to all elevational bands between the lowest and highest observed presences. In this study, the elevational ranges of species were then transformed by "n × 300" m (where "n" is the number of elevational bands that a species occurs). This interpolation method can reduce bias of the underestimation of species richness caused by insufficient sampling, especially in areas with high biodiversity (Wu et al., 2013). This interpolation method is widely used in elevational research and allows for comparisons with other relevant studies (Brehm, Colwell, & Kluge, 2007;Rowe, 2009;Wu et al., 2013).
We used three methods to assess how well the species diversity was sampled. Species accumulation curves have often been used to evaluate whether species diversity was sampled adequately; an adequate survey of species was assumed if the species accumulation curve reach an asymptote (Magurran & McGill, 2011). In the first method (method 1), we randomized the accumulation order of individuals 50 times by EstimateS 9.10 (Colwell, 2013; https://purl. oclc.org/estimates/) and obtained the individual-based rarefaction curves (cumulative species number as a function of individual number). Individual-based rarefaction is sensitive to biases in the quantification of the number of individuals per species (Gotelli & Colwell, 2001), while species richness estimation based on samples is less F I G U R E 1 Location of the study area. The study area encompasses 12 elevational sampling bands. Sample sites display the midpoints of each transect sensitive to this problem (Herzog et al., 2005). Thus, the second method (method 2) is sample-based rarefaction. We used a modified version of the "m-species-list method" (Poulsen, Krabbe, Frølander, Hinojosa, & Quiroga, 1997) to generate small samples in each 300-m elevational band (each m-species list was treated as a separate sample). We combined all field surveys of each 300-m elevational band into a master list of bird observations (each transect was surveyed separately, and records were joined according to the sequence of the survey). We then divided the master list of each 300-m elevational band into lists of 10 species: The first 10-species list consists of the first 10 species observed, and the second 10-species list may include repetitions or new species compared with the first 10-species list. We randomized the sample accumulation order 50 times using EstimateS and obtained the sample-based rarefaction curves (cumulative species number as a function of list number).
However, it is unlikely to detect all species in natural communities even after thorough and extensive sampling. Therefore, in the third method (method 3), estimators (MMMean and Chao2) were used to compute estimated species richness. For species-rich bird data sets, MMMean was the least biased, but MMMean could not reflect the statistically variance (Herzog, Kessler, & Cahill, 2002).
Chao2 could be used as supplementary, as it is the least biased estimates of species richness for small numbers of samples (Colwell & Hurtt, 1994). The MMMean and Chao2 statistic of each 300-m elevational band were obtained from the sample-based rarefaction in EstimateS. To compare survey effort in different 300-m elevational bands, observed richness was expressed as the proportion of the respective MMMean statistic (Herzog et al., 2005).

| Measuring beta diversity
The Simpson dissimilarity index (Simpson, 1943), which reflects species compositional differences (or describes spatial turnover) without the influence of richness gradients (Leprieur et al., 2011), has been used as a measurement of beta diversity (or turnover) in elevational studies (Lennon, Koleff, Greenwood, & Gaston, 2001;McCain & Beck, 2016;Mena & Vazquez-Dominguez, 2005). The index is defined as: where a is the number of species shared by the two elevational bands, and b and c are the number of species that only appear in each of the different elevational bands. In this study, we use the β sim between all pairs of 300-m elevational bands to reflect the change in species composition along the elevational gradient in Gyirong Valley.

| Environmental factors
We collected data related to dispersal limitations, geometric constraints, habitat complexity, contemporary climate, and historical climate for Gyirong Valley. Each of the environmental factors is described in detail below:

| Area and MDE
The planimetric area of each 300-m elevational band was calculated in ArcGIS 10.4 (ESRI, Redlands, CA, USA) using 30-m digital elevation model (DEM). The DEM data were acquired from the Geospatial Data Cloud website (GDC; https://www.gscloud.cn). The values of the area were much higher than the values of other variables, and then we log(x)-transformed the values of area to alleviate heteroscedasticity.
We randomized (without replacement) the empirical species ranges within the bounded domains (from 1,800 to 5,400 m asl) to generate a predicted species richness (R MDE ) under geometric constraints using RANGEMODEL 5 (https://purl.oclc.org/rangemodel; Colwell, 2008). R MDE and their 95% confidence intervals were computed for each 300-m band based on 10,000 simulations of the MDE rang-model. We also used the Beta Simulation program (https://spot.colorado.edu/∼mccainc/simulation_programs.

| Contemporary climate
To catch the real characteristics of climate in such fine-scale mountainous region, we established six mini weather stations in Gyirong Valley (2,457,2,792,3,368,3,740,4,140, and 5,230 m asl; Figure 1).
Each mini weather station consists of three data loggers (HoBo Pro-RH/Temp, HoBo Pro-Precipitation/Temp, and HoBo Pro-PAR) and was surrounded by fences to prevent interference from wild animals.
Mean daily temperatures (MDT), precipitation (P), relative humidity, and photo-synthetically active radiation were measured and recorded from September 2015 to July 2016. Only temperature and precipitation data which show more consistent spatial structured signals were used as the climatic factors in this study. Temperature data were recorded every 10 min and were averaged afterward to 1 hr to minimize the impact of possible outliers from September 2015 to July 2016 (Brehm et al., 2007). Precipitation data were accumulated also from September 2015 to July 2016. Then, the MDT and P recorded by the mini weather stations were extrapolated to all elevational bands in the study area using ArcGIS 10.4.

| Habitat complexity
We combined a 30-m DEM and the 300-m GlobCover landcover data of the study to obtain the landcover type in each 300-m elevational band. We extracted the cells of the landcover raster data that correspond to the areas defined by the 30-m DEM for each 300-m elevational band by using the extract by mask tool in ArcGIS 10.4.
We set the cell size of the result as the inputted DEM data (30-m).
Excluding the artificial areas, 21 landcover types are defined following the United Nations Land Cover Classification System (LCCS). The , landcover data were from Globcover2009 (https://www.gscloud. cn/, data were accessed in 25th October 2015). We calculated the habitat heterogeneity (HH) using the Shannon diversity index.
Plant species richness represented structural complexity of the habitat (Hu, et al., 2017

| Historical climate
The Quaternary Period covers the last 1.8 million years, characterized with frequent glacial advances and retreats (Leprieur et al., 2011

| Statistical analysis
Polynomial regression analyses were performed to assess the form of elevational patterns of β sim . We compared β sim within the two adjacent 300-m elevational bands with the combined species richness of those two elevational bands by using linear regressions. We also compared β sim with β MDE to test whether random species community structure contribute to the beta diversity patterns by using linear regressions and assessing whether the elevational gradient had any β sim outside the confidence intervals of the MDE model.
To be consistent, we preserved the same dissimilarity measure (β sim ) used in the remaining statistical analyses. We conducted distance-based redundancy analysis (db-RDA) to explain the beta diversity by relating the bird species composition to the eight environmental factors. A site × species matrix and a site × variables (environmental factors) matrix were used in db-RDA. Two subsequent approaches were conducted to search for parsimony and reduce the possible strong linear dependencies (collinearity) among the predictor variables. In the first approach, we used Pearson correlation to examine the relationships among the eight environmental factors and excluded the highly correlated factors based on the correlation coefficients. MDT were highly correlated with P (r = 0.916, p < 0.001) and Area (r = −0.944, p < 0.001); HH were highly correlated with MDE (r = 0.957, p < 0.001); and PC were highly correlated with Area (r = −0.871, p < 0.001) and TC (r = 0.851, p < 0.001; Supporting information Appendix S1: Table S1). After re- To analyze variation in beta diversity, Mantel tests were conducted to explain the correlation between species dissimilarity and dissimilarity of the variable subsets. We quantified the dissimilarity in species (β sim ) and dissimilarity in four variable subsets (variable-subset.diff) between all pairs of 300-m elevational bands and arranged the dissimilarities into sites × sites triangular "dissimilarity" matrixes (Keil et al., 2012;Winter et al., 2010

| Fauna and sampling effort
We recorded a total of 182 bird species (belonging to 12 orders, 43 families, and 105 genera), including 169 breeding bird species (belonging to 11 orders, 41 families, and 100 genera; Supporting information Appendix S1: Table S2). All individual-based rarefaction curves for each 300-m elevational band reached a plateau or an asymptote (Figure 2a). The sample-based rarefaction curves were steeper compared to the individual-based curves, but most still reached an asymptote (except the 12th elevational band ranging from 5,100 to 5,400 m asl, presumably caused by the limited numbers of species or 10-species lists in this band; Figure 2b; Supporting information Appendix S1: Fig. S1). Although empirical species richness values were lower than the estimators (MMMean and Chao2), the survey effort (observed richness/ MMMean) of each 300-m elevational band was still higher than 70% ( and as such, we conclude that our data processing approaches were appropriate in this study and the richness data can reflect F I G U R E 2 Individual-based (a) and sample-based (b) rarefaction curves for bird data sets for each elevational band in Gyirong Valley. The numbers from 1 to 12 represent the twelve 300-m elevational bands from 1,800 to 5,400 m asl; for example, "1" is the 300-m elevational band ranging from 1,800 m to 2,100 m asl. Accumulation order of all curves was randomized 50 times

| Diversity patterns of birds
Species richness showed a hump-shaped pattern along the elevational gradient in Gyirong Valley, with a peak at 2,400-3,000 m asl (76 species, Figure 4). Qualitatively, the elevational patterns of β sim between adjacent 300-m elevational bands and combined species richness of those pairs of adjacent bands both peaked at intermediate elevations (with a second peak; Table 2). Quantitatively, the elevational pattern of β sim was hump-shaped (r 2 first-order = 0.0334, second-order = 0.536, p < 0.05; r 2 third-order = 0.606, p > 0.05) and the correlation between β sim and species richness was significant (r 2 = 0.513, p < 0.01).

F I G U R E 3 Environmental patterns across elevation in
The elevational pattern of β MDE predicted by the MDE model was also hump-shaped, peaking at 3,300-3,900 m asl (Figure 4 and Table 2; r 2 first-order < 0.01, p > 0.05; r 2 second-order = 0.837, p < 0.01; r 2 third-order = 0.837, p < 0.01). The correlation between the simulated β MDE pattern and empirical β sim pattern was moderate but significant (r 2 = 0.446, p < 0.05). Most of the empirical β sim values were higher than the simulated values with five empirical β sim points (two peaks) falling outside the 95% confidence intervals of the MDE model ( Figure 4).

| Factors explaining beta diversity patterns
The results of the db-RDA with all eight environmental variables r 2 HH = 0.561, p < 0.05; r 2 Area = 0.913, p < 0.01; r 2 MDE = 0.593, p < 0.05; r 2 TC = 0.500, p < 0.05; r 2 PC = 0.808, p < 0.01). After removing those highly correlated factors, db-RDA with five variables (P, PR, Area, MDE, and TC; db-RDA.select) explained 96% of the variation in species composition (λ Total = 1.977, λ constrained = 1.904). Axis No. 1 (λ 1 = 1.703) explained 86% of the total variation in the data set and 89% of the variation explained by the five constrained axes; axis No. 2 (λ 2 = 0.185) explained 9% of the total variation and 10% of the variation explained by the constrained axes (Figure 5b) We performed variance partitioning based on two db-RDA models (db-RDA.select and db-RDA.step) and partitioned the overall explained variance into components of independent (I) and joint spatial constraints (SFs). For db-RDA.select, the overall explained variance was partitioned into components of independent and joint effects of EF 1 (P, PR, and TC) and SF 1 (Area and MDE). The results of variation partitioning based on db-RDA.select showed that EF 1 explained more of the overall variance of species dissimilarity than SF 1 (I EF1 = 0.0932, I SF1 = 0.0839, and J = 0.756; Figure 7a). For db-RDA.
step, variance was partitioned into components of independent and joint effects of P (EF 2 ) and Area (SF 2 ). The results of variation partitioning based on db-RDA.step showed that P explained more of the overall variance of species dissimilarity than Area (I EF2 = 0.243, I SF1 = 0.104, and J = 0.581; Figure 7b).

| Elevational patterns of beta diversity
In this study, we found that beta diversity along the elevational gradient showed a hump-shaped pattern, peaking at intermediate ele- vations (and with a second peak at higher elevations, Figure 4). There is a long-standing hypothesis that ecosystems with high beta diversity (or turnover) are accompanied by high biodiversity (Clements, 1916;Stevens, 1992;Whittaker, 1975;Wilson & Shmida, 1984), and beta diversity has been confirmed to contribute to overall species richness patterns (Buckley & Jetz, 2008;Davies, Scholtz, & Chown, 1999;Fattorini, 2014;Herzog et al., 2005;McCain & Beck, 2016;Naniwadekar & Vasudevan, 2007). Our study provides evidence to support the hypothesis that beta diversity is a driver of species richness along elevational gradients: beta diversity was consistent with species richness (r 2 = 0.626, p < 0.01).
Beta diversity patterns predicted by the MDE null model (β MDE ) were also hump-shaped, consistent with the empirical β sim patterns ( Figure 4). However, the predictions of mid-domain null models were also not unified, with the predicted beta diversity (or turnover) patterns showing shallow-unimodal patterns, bimodal patterns, and U-shaped patterns due to the difference of method, scale, and study area (McCain & Beck, 2016). Although the predictions of the null models were inconsistent, most of the studies could not reject the null prediction (few empirical beta diversity values fell outside the 95% confidence intervals of the MDE null model predictions; McCain & Beck, 2016;Mena & Vazquez-Dominguez, 2005). Only one exception rejected the random structuring assumption (Herzog et al., 2005). In our study, five (of 12) empirical β sim points fell outside the 95% confidence intervals of the MDE null model, but with a moderate (significant) correlation between the empirical and MDE-modeled patterns. These indicated that overall shape of the patterns may be correlated with the MDE, but that the peaks of beta diversity were not explained by the MDE model. We should note that spatial constraints and stochastic processes may play some role in shaping the beta diversity pattern for birds along the elevational gradient in Gyirong Valley.
Some studies have found that beta diversity is related to ecotones (Jankowski et al., 2009;Vazquez & Givnish, 1998). In Gyirong Valley, beta diversity peaked at 2,400-3,000 m asl (the transition zone between evergreen broadleaf forest and broadleaf mixed forest), and 3,600-4,200 m asl (the transition zone between dark coniferous forest and shrub and grass habitat). Moreover, the region from

| Effects of environmental factors on beta diversity
In this study, results of the db-RDA models suggested that species were distributed along a climate-habitat gradient. Additionally, our study found high rates of beta diversity in the ecotones with more precipitation and higher temperature (at 2,400-2,700 m:  (Keil et al., 2012). This should call for caution that we should use higher resolution data to cope with ecological studies in fine-scales, especially in mountainous areas with high biodiversity, complex landscapes, and variable climate.
According to previous studies, effects of environmental filtering are stronger at the local scale, while dispersal limitations play an important role at the larger regional scale (Rejmanek, 2000). In Gyirong Valley, the reason why the main driver for beta diversity patterns was environmental filtering (rather than dispersal limitations or other) may simply be that elevational bands along the elevational gradient were not sufficiently isolated, and the mobility of avian species is quite strong. The distances along the elevational gradient in Gyirong Valley (the planimetric distance from the bottom of the valley to the summit of Mt. Mala is only 79 km) might not be great enough to prevent the birds from dispersing. Still, another reason might be that the environment in Gyirong Valley is heterogeneous and variable. Habitats with different environments should harbor different sets of species, and the more different the environment, the greater the beta diversity (turnover) would be.

| CON CLUS ION
In general, environmental filters, habitat heterogeneity, spatial constraints, and stochastic processes influence beta diversity

DATA ACCE SS I B I LIT Y
The raw data have been supplied as Supplementary Files.