Altitudinal patterns of plant diversity on the Jade Dragon Snow Mountain, southwestern China

Background Understanding altitudinal patterns of biological diversity and their underlying mechanisms is critically important for biodiversity conservation in mountainous regions. The contribution of area to plant diversity patterns is widely acknowledged and may mask the effects of other determinant factors. In this context, it is important to examine altitudinal patterns of corrected taxon richness by eliminating the area effect. Here we adopt two methods to correct observed taxon richness: a power-law relationship between richness and area, hereafter “method 1”; and richness counted in equal-area altitudinal bands, hereafter “method 2”. We compare these two methods on the Jade Dragon Snow Mountain, which is the nearest large-scale altitudinal gradient to the Equator in the Northern Hemisphere. Results We find that seed plant species richness, genus richness, family richness, and species richness of trees, shrubs, herbs and Groups I–III (species with elevational range size <150, between 150 and 500, and >500 m, respectively) display distinct hump-shaped patterns along the equal-elevation altitudinal gradient. The corrected taxon richness based on method 2 (TRcor2) also shows hump-shaped patterns for all plant groups, while the one based on method 1 (TRcor1) does not. As for the abiotic factors influencing the patterns, mean annual temperature, mean annual precipitation, and mid-domain effect explain a larger part of the variation in TRcor2 than in TRcor1. Conclusions In conclusion, for biodiversity patterns on the Jade Dragon Snow Mountain, method 2 preserves the significant influences of abiotic factors to the greatest degree while eliminating the area effect. Our results thus reveal that although the classical method 1 has earned more attention and approval in previous research, method 2 can perform better under certain circumstances. We not only confirm the essential contribution of method 1 in community ecology, but also highlight the significant role of method 2 in eliminating the area effect, and call for more application of method 2 in further macroecological studies. Electronic supplementary material The online version of this article (doi:10.1186/s40064-016-3052-1) contains supplementary material, which is available to authorized users.

patterns along altitudinal gradients, i.e., taxon richness peaks at mid-elevations, for a broad spectrum of taxa (Rahbek 1995(Rahbek , 2005. Roughly 25 % of all published investigations have shown a monotonic decrease in taxon richness with elevation (Rahbek 1995;Kessler 2002). A few studies have observed approximately constant values from lowlands to mid-elevations, followed by a pronounced fall (Rahbek 2005;Kessler 2001). The mechanisms determining altitudinal patterns of taxon richness are still under debate (Currie 1991;Jetz and Rahbek 2002). It has been broadly recognized that area is a decisive factor shaping altitudinal taxon richness patterns (Rahbek 1997) and that the area of altitudinal bands alone could account for a large percentage of the variation in taxon richness (Bachman et al. 2004;Fu et al. 2004;Kattan and Franco 2004;McCain 2007).
Although area is admittedly important in determining altitudinal patterns of plant diversity (Rosenzweig 1995), it is also widely accepted that area is not the sole factor influencing the variation of taxon richness along the altitudinal gradient (Bachman et al. 2004;Wang et al. 2007a, b). The crucial mechanisms behind altitudinal patterns of taxon richness include the combined effects of area, climate, and the mid-domain effect (MDE;McCain 2009). Mean annual temperature (MAT) and mean annual precipitation (MAP) could also be responsible for a large part of the variance in plant taxon richness (O'Brien et al. 1998;Baudena et al. 2015). The MDE has been proposed as one important determinant for hump-shaped patterns of taxon richness along the altitudinal gradient (Colwell et al. 2004;Kluge et al. 2006). According to the MDE, species' ranges are randomly placed within a geographical domain under the constraint that no species can exceed the hard boundaries of this domain (Colwell and Lees 2000;Colwell et al. 2004). We can simulate this stochastic procedure multiple times, and the mean number of species generated by these multiple simulations is considered to be a prediction of MDE. Some findings have shown that the combined influences of area and MDE account for virtually all variance in taxon richness along altitudinal gradients (Sanders 2002;Bachman et al. 2004). Thus, the key to understanding the underlying mechanism is disentangling the relative importance of area and MAT, MAP and MDE (Rahbek 2005;Sanders et al. 2007). However, the practical difficulty is that area, MAT, MAP, and the prediction from MDE usually covary along the altitudinal gradient (Wang et al. 2007a, b). Such collinearity might mask the possible importance of MAT, MAP, and MDE. It is therefore necessary to accurately quantify the effect of area on altitudinal patterns of plant diversity and to evaluate the effects of MAT, MAP, and MDE on the patterns after eliminating the area effect.
Species-area relationships are often applied, and are generally considered to account for area effects on altitudinal patterns of taxon richness (Rosenzweig 1995). Meza-Joya and Torres (2016) found that a power-law curve was the best-fit species-area model on the Tropical Andes and its domains. A large-scale study of vascular plant richness in North America (Qian 1998) used an exponential relationship between richness and area to eliminate the area effect on taxon richness. The species-area relationship is clearly a powerful tool for adjusting taxon richness to account for the area effect, and we refer to it as method 1 in our analysis.
Another novel yet concise method was first proposed by Bachman et al. (2004), accounting for the area effect through the use of equal-area altitudinal bands. We refer to this as method 2. Bachman et al. (2004) clarified that when assessed in equal-elevation bands, palm species richness appeared to drop monotonically with elevation, while if evaluated in equal-area bands, species richness showed a noticeable hump-shaped pattern. The reason for this was the huge percentage of lowlands in New Guinea (Bachman et al. 2004). Although this equal-area band methodology is effective in directly eliminating the area effect, it has seldom been applied to account for the area effect so far (Zhu et al. 2007;VanDerWal et al. 2008;Xing et al. 2011). Zhu et al. (2007 demonstrated that area along the equal-elevation gradient of Helan Mountain declined monotonically with elevation, and vascular plant species richness displayed hump-shaped patterns along both equal-elevation and equal-area gradients. VanDerWal et al. (2008) evaluated the effect of MDE on amphibian, bird, mammal and tree species richness along equal-area altitudinal bands in North America, and found that MDE could explain the observed pattern well for any taxa considered.
To our knowledge, the above two methods have not been implemented simultaneously for the same altitudinal gradient to date. In this study, we examine altitudinal patterns of plant diversity on the Jade Dragon Snow Mountain in southwestern China and use both methods to eliminate the area effect. The Jade Dragon Snow Mountain encompasses a broad altitudinal gradient with rich flora, diverse climate, and extensive field surveys of plant distribution over recent years (Wang et al. 2007a, b). It thus provides a valuable opportunity to explore seed plant diversity patterns along an altitudinal gradient, and the impacts of MAT, MAP, and MDE on the patterns with and without the area effect.

Study area and data collection
The Jade Dragon Snow Mountain and its adjacent region (26°35′N-27°45′N, 99°22′E-100°32′E) are located in the northwest of Yunnan Province, China, bordering the Tibetan Plateau (Fig. 1). This region covers a total area of 6127 km 2 and has a large altitudinal gradient from 1350 to 5050 m. The climate in this region is under the control of the southwest monsoon from the Indian Ocean, with richer precipitation in the eastern part than in the western part (Wang et al. 2007a, b). Seventy percent of the MAP occurs between May and September, with winter precipitation only contributing 30 % of the annual total. The Jade Dragon Snow Mountain, with its adjacent region, is a global biodiversity hotspot (Myers et al. 2000), where the flora is rich, including a total of 2028 native seed plant species from 625 genera and 136 families. More than half of the species are herbaceous (66.4 %; 82 families; 426 genera; and 1346 species), and 33.6 % are woody species (290 tree species from 120 genera and 56 families, and 392 shrub species from 133 genera and 62 families; Wang et al. 2007a, b). The vegetation types on the Jade Dragon Snow Mountain from low to high elevations are tropical forest, subtropical forest, alpine meadow and alpine tundra, which corresponds to the general altitudinal pattern of plant diversity in northwestern Yunnan Province.
We generated a species database from the published book Checklist of Seed Plants of Lijiang Alpine Botanic Garden (Wang et al. 2007a, b), which is based on substantial surveys on the Jade Dragon Snow Mountain and in its adjacent region. This species database includes minimum and maximum elevations of occurrence for each species, and lifeform (tree, shrub, and herb) information on each species. We extracted topographical data from ASTER GDEM V2 (Global Digital Elevation Model Version 2, DOI:10.5067/ ASTER/ASTGTM.002) with 30 m × 30 m resolution. We obtained MAT and MAP data from the WorldClim database (Hijmans et al. 2005), which is frequently used in ecological studies (e.g., Sommer et al. 2010).

Equal-elevation altitudinal gradient
We divided the Jade Dragon Snow Mountain into 37 100-m altitudinal bands from 1350 to 5050 m. Along this equal-elevation gradient, we calculated the area of each altitudinal band by multiplying 900 m 2 by the number of digital elevation model (DEM) grids in each band. The area increases steeply with increasing elevation, and then decreases above the 2650-2750 m altitudinal band, showing a hump-shaped pattern (Fig. 2a). MAT declined monotonically with elevation ( Fig. 2b). MAP illustrated a much more complex pattern with elevation, increasing below 1900 m, then declining steeply to 912.0 mm at 4200 m and afterward increasing to the top of the Jade Dragon Snow Mountain (Fig. 2c).
We assumed that each species had a continuous distribution range between its recorded minimum and maximum elevations, as widely used in previous studies (e.g., Rahbek 1997;Vetaas and Grytnes 2002;Sanders 2002). However, among the 2028 seed plant species, there are 717 species that have been recorded only in a single elevational band. We counted the number of seed plant species, genera, and families present at each altitudinal band as observed species, genus and family richness (Sobs, Gobs, and Fobs). To explore altitudinal patterns of plant diversity for different life-forms, we counted the number of tree, shrub, and herb species occurring at each altitudinal band as the species richness for trees, shrubs, and herbs (TSobs, SSobs, and HSobs). Additionally, in order to evaluate whether the altitudinal biodiversity patterns and their determinant predictors were dependent on species' range size, we divided all species into three groups and 500 m (Group II, 554 species), and >500 m (Group III, 644 species). These range size limits were chosen to make the number of species in each group comparable. In the same way, we obtained the species richness for Groups I-III (ISobs, IISobs, IIISobs) by counting the number of species in Groups I, II, and III present in each altitudinal band. We thus generated nine observed taxon richness variables (hereafter TRobs): Sobs, Gobs, Fobs, TSobs, SSobs, HSobs, ISobs, IISobs, and IIISobs.

Effect of interpolation
Like many previous studies, we considered a species to be present at all elevations within its recorded elevation limits (Kluge et al. 2006). However, this approach may produce artificially elevated species richness at mid-elevations, since such interpolated data are disproportionately added to mid-elevations as opposed to edges of the gradient (Karger et al. 2011). In order to find out whether using interpolated species richness masks its real pattern along the altitudinal gradient, we evaluated the relationship between elevation and species richness of a particular plant species set (717 species that have been recorded only in a single elevational band) without interpolation (Vetaas and Grytnes 2002). We checked whether this species richness pattern generated without interpolation manifested a hump-shaped curve. If it indeed shows a hump-shaped pattern, we can conclude that the effect of interpolation may not be essential for the hump-shaped pattern of plant species on the Jade Dragon Snow Mountain.

Species-area relationship
The species-area relationship has been universally acknowledged, but the exact structure of the relationship is still under discussion (Connor and McCoy 1979;Crawley and Harral 2001). Three common versions of the species-area relationships are: untransformed (species richness versus area), semi-log (species richness versus log area), and log-log (log species richness versus log area) transformed (Matthews et al. 2014). We conducted linear fittings to all three versions using ordinary least square (OLS) regression models, which was commonly used for studies on the species-area relationships (Matthews et al. 2014). To select the version (i.e., untransformed, semi-log or log-log) with the best performance, we calculated the modified Akaike information criterion (AIC C ) corrected for small samples as follows: where n is the number of samples and K the number of parameters in the model. We considered the model with the lowest AIC C score to be the best model, which is consistent with previous studies (Matthews et al. 2014). Then we compared the difference between the AIC C of each model and the minimum AIC C found, and we refer to this difference as ∆(AIC C ). Any model with ∆(AIC C ) < 2 is reported to be as good as the best model (Burnham and Anderson 2002). This analysis was performed using the function ' AICc' within the ' AICcmodavg' package (Mazerolle 2011) in R 2.14.2 (R Core Team 2014). To further assess the performances of linear fits, we also measured the adjusted coefficients of determination and conducted the significant F-test (Crawley 2002).
The log-log transformed species-area relationship was found to be the best-fit model (Δ(AIC C ) > 2) for all plant groups (Additional file 1: Table S1; Figures S1-S3). We thus chose the log-log transformed version (i.e., power-law) to correct TRobs. This powerlaw model, S = cA z , where S is TRobs, z is a constant describing the slope of the speciesarea relationship in the log-log transformed space (Additional file 1: Fig. S3), A is the area of elevational bands along the equal-elevation gradient and c is the area-corrected taxon richness (hereafter TRcor 1 ). In order to rescale TRcor 1 to similar values as TRobs, TRcor 1 were calculated as 100(TRobs/A z ).

Equal-area altitudinal gradient
Another approach to account for area effect was introduced by Bachman et al. (2004), and we refer to this as method 2 in our analysis. The original DEM cell values are integers and were added a random number between −0.5 and +0.5, as done in previous studies (Bachman et al. 2004). This produced DEM cells are easily classified into 37 equal-area altitudinal bands from 1350 to 5050 m using software ArcGIS version 9 (Environmental Systems Research Institute, Redlands, CA, USA). The middle elevations of the equalarea bands were not uniformly distributed along the altitudinal gradient (Fig. 3). The altitudinal band width first decreased then increased as elevation rose (Fig. 4a). The band with the smallest width (35.4 m) was the 15th, while the largest width was 1014.1 m for the 37th band. MAT still declined monotonically, and MAP showed a similar but simpler pattern along the equal-area gradient compared with the equal-elevation gradient (Figs. 2b,c,4b,c). Along the equal-area altitudinal gradient, we counted the number of (1) seed plant species at each band as Scor 2 . In the same way, we got Gcor 2 , Fcor 2 , TScor 2 , SScor 2 , HScor 2 , IScor 2 , IIScor 2 , and IIIScor 2 for other plant groups. All these nine areacorrected taxon richness variables were collectively called TRcor 2 .

Mid-domain effect
We used RangeModel (Colwell 2008) to generate the simulated taxon richness, which are the mid-domain null model predictions. This software placed empirical species ranges within the domain (i.e. the mountain range from the lowest elevation to the peak)  randomly and under the constraint that no species extended beyond domain boundaries (Colwell and Lees 2000;Colwell et al. 2004). Then the number of species was counted within each elevational band. We conducted 1000 simulations, and the mean of these simulations was called a prediction from MDE. This procedure was carried out along the equal-elevation and equal-area gradient.

Statistical analysis
We first conducted OLS regressions between taxon richness (TRobs, TRcor 1 , and TRcor 2 ) and each explanatory variable (elevation, MAT, MAP, and prediction from MDE), respectively (Additional file 1: Figures S4-S6). In order to compare the performances of first-and second-order polynomial regressions, we calculated ∆(AIC C ) according to Eq. (1). Model with the minimum AIC C and ∆(AIC C ) > 2, was selected as the best one. Two models with ∆(AIC C ) < 2 were considered to have the same good performance, and in this case, the first-order polynomial regression was selected as the best model for simplicity (Additional file 1: Tables S2-S4). Secondly, we examined the spatial autocorrelation in the residuals of the best OLS model using Moran's I coefficient. We regard OLS model residuals in the 37 sequential elevational bands (from 1350 to 5050 m) as 37 observations along a single geographic axis. Moran's I coefficients are computed from pairs of observations found at preselected distances: distance = 100 m, 200 m, 300 m, etc. (Legendre and Legendre 1998). Moran's I coefficient is defined as follows: where n is the number of elevational bands, x i and x j represent observations in elevational bands i and j, x is the mean of all x, and w ij is an element in the (n × n) weighting matrix W. It can be given as follows: S represents the sum of the weights w ij (i.e., the number of connections in the matrix W) as follows: Moran's I coefficient varies between −1.0 and 1.0 for maximum negative and positive spatial autocorrelation, respectively. Non-zero values of Moran's I coefficient indicate that observations in elevational bands connected at a given distance are more similar (positive autocorrelation) or less similar (negative autocorrelation) than expected by chance. Through plotting Moran's I coefficients against the preselected distances, we constructed spatial correlograms and detected considerable spatial autocorrelation in OLS model residuals (Additional file 1: Figures S7-S18). This analysis was performed using function 'correlog' within the package 'ncf ' (Bjørnstad 2006).
However, spatial autocorrelation has been reported to inflate Type I errors and thus lead to the biased model comparison and poor parameter estimates, through violating assumptions of independence and identical distribution of model residuals (Dormann et al. 2007). Therefore, thirdly we ran simultaneous autoregressive (SAR) models of the error type  to correlate taxon richness (TRobs, TRcor 1 , and TRcor 2 ) and each explanatory variable (elevation, MAT, MAP, and prediction from MDE), respectively. Spatial weights matrices in SAR were based on row standardization and neighborhood distance of 100 m. Pseudo-R 2 values of SAR were calculated as the squared Pearson product-moment correlation coefficient between observed and predicted values . We implemented SAR models with package 'spdep' (Bivand 2014). We also performed first-and second-order polynomial SAR models and compared their performances using ∆(AIC C ). The selection criterion of the best model was exactly the same as OLS analysis. The best SAR models reduced spatial autocorrelation in the model residuals to a lower level, especially at the distance of 100 m (Additional file 1: Figures S7-S18).
There is no consensus on which area-correction method is better. Each method has its advantages and drawbacks. In this study, we try to find which area-correction method preserves the influences of other factors to a larger degree after eliminating the area effect. Therefore, fourthly we compared SAR pseudo-R 2 values of method 1 and 2, and consider the method with the larger pseudo-R 2 values to have the better performance.

Altitudinal patterns of observed taxon richness
Without accounting for area, Sobs, Gobs, and Hobs showed hump-shaped patterns along the equal-elevation altitudinal gradient: richness increased steeply with elevation at low elevations, and then decreased at high elevations after peaking at intermediate elevations (between 2750 and 2850 m for both species richness and genus richness, and between 2550 and 2650 m for family richness) ( Fig. 5a; Table 1). TSobs, SSobs, HSobs also showed hump-shaped curves along the equal-elevation gradient, with maxima at approximately the same altitude (2550-2850 m; Fig. 5b). Moreover, ISobs, IISobs, and IIISobs also illustrated hump-shaped patterns along the equal-elevation gradient (Fig. 5c). Maximum species richness for Group I and Group II appeared at the 2750-2850 m elevation interval, whereas Group III maximum species richness occurred at 2950-3050 m.

Effect of interpolation
According to Vetaas and Grytnes (2002), we should rule out the possibility that the interpolation approach generates an artificial hump-shaped pattern in species richness. In this study, a particular species set (717 species that have been recorded only in a single elevational band) provides an appropriate opportunity to check whether the humpshaped species richness pattern found in Jade Dragon Snow Mountain is real. Without interpolation, this particular species set indeed showed a hump-shaped pattern in species richness with a peak at 2800 m and changes in elevation could explain 40.2 % of the variation in species richness (Fig. 6). This implies that the hump-shaped pattern for plant species in the Jade Dragon Snow Mountain may not due solely to interpolation.

Table 2 Simultaneous autoregressive models fit area-corrected taxon richness achieved by method 1 against the first-and second-order polynomials of four variables (elevation, mean annual temperature, mean annual precipitation, and prediction from mid-domain effect)
Area-corrected taxon richness 1 n First-order Second-order
Once we excluded the effect of the interpolation method in producing the humpshaped patterns in taxon richness as suggested by Vetaas and Grytnes (2002), we could explicitly assess the effect of abiotic factors on the patterns. Variations in taxon richness are rarely wholly due to a single factor (Oommen and Shanker 2005). Area has a strong influence on the pattern of species richness at all scales (Whittaker et al. 2001). Larger areas contain more individuals and thus more species ('passive-sampling hypothesis' , Connor and McCoy 1979). The analyses on the Jade Dragon Snow Mountain also yielded further evidence for the significant influence of area on these patterns. A power-law species-area relationship explained a good deal (≥53.1 %) of the variance in TRobs (Additional file 1: Table S1; Figures S1-S3). While accounting for area, TRcor 2 still displayed Fig. 8 Altitudinal patterns of area-corrected a seed plant species richness (Scor 2 ), genus richness (Gcor 2 ), family richness (Fcor 2 ), b tree species richness (TScor 2 ), shrub species richness (SScor 2 ), herb species richness (HScor 2 ), c Group I species richness (IScor 2 ), Group II species richness (IIScor 2 ), and Group III species richness (IIIScor 2 ) along the equal-area gradient. The dashed, dotted, and dash-dotted lines represent the elevations at maximum richness. For abbreviations, see Fig. 5

Table 3 Simultaneous autoregressive models fit area-corrected taxon richness achieved by method 2 against the first-and second-order polynomials of four variables (elevation, mean annual temperature, mean annual precipitation and prediction from mid-domain effect)
Area-corrected taxon richness 2 n First-order Second-order
Altitudinal gradients in taxon richness have usually been attributed to a linear relationship with MDE predictions (Colwell and Lees 2000;Colwell et al. 2004), which is confirmed by our findings (Table 1; Additional file 1: Fig. S4). While accounting for the area effect, MDE accounts for less of the variation in area-corrected taxon richness compared with observed taxon richness, revealing collinearity between area and MDE. This is inconsistent with the case of palms in New Guinea (Bachman et al. 2004). Bachman et al. clarified that MDE predictions conditionally explain up to 98 % of the variance in palm taxon richness after removing the effect of area, since in New Guinea both area and observed taxon richness of palms along the equal-elevation gradient decrease monotonically with elevation (Bachman et al. 2004). However, taxon richness along the equal-area gradient showed a hump-shaped pattern, which was more consistent with the changes in MDE with elevation. Thus an area correction approach increased the explanatory power of MDE for palms in New Guinea. However, on the Jade Dragon Snow Mountain, observed taxon richness and simulated taxon richness predicted by MDE all show hump-shaped patterns along the equal-elevation altitudinal gradient. After removing the area effect by method 2, i.e., along the equal-area altitudinal gradient, both areacorrected taxon richness and simulated taxon richness evaluated by MDE still exhibit hump-shaped patterns. The strength of the relationship between TRcor 2 and MDE is weakened when the area effect is eliminated. Except for trees, the prediction from MDE accounted for a larger part of the variation in TRcor 2 than in TRcor 1 . This indicates that method 2 retains the influence of MDE to the maximum large extent as well as eliminating the area effect.
For patterns of plant diversity along the altitudinal gradient on the Jade Dragon Snow Mountain, our analysis implies that method 2 performs better than method 1 in terms of preserving the significant effects of MAT, MAP, and MDE. There is no doubt that the species-area relationship is one of ecology's few laws (Rosenzweig 1995;Rahbek 1997;McCain 2007). Thus method 1, based on this law, should perform well in accounting for the area effect, and this is indeed confirmed by our results ( Fig. 9; Tables 1, 2, 3). In this paper, we not only appreciate the essential contribution of method 1 to community ecology, but also highlight the significant role of method 2 in macroecological studies conducted in mountainous regions.

Conclusions
In this study, we investigated the altitudinal patterns of taxon richness for seed plants on the Jade Dragon Snow Mountain, southwestern China. We evaluated the effect of area on taxon richness patterns (TRobs) and invoked two methods to calculate areacorrected taxon richness (TRcor 1 and TRcor 2 ). We calculated the number of species in the smallest sampling area based on a power-law species-area relationship implemented in method 1. According to method 2, we counted the number of taxa along the equalarea altitudinal gradient. We assessed the influences of MAT, MAP, and MDE on taxon richness before and after eliminating area effect. Our results reveal that both TRobs and TRcor 2 show hump-shaped patterns along the altitudinal gradient, while TRcor 1 dose not. Elevation, MAT, MAP, and MDE explain a smaller proportion of the variance in area-corrected taxon richness than in observed taxon richness. Moreover, they were all responsible for a larger percentage of the variance in TRcor 2 than TRcor 1 . These findings indicate that area effects should be taken into consideration when assessing the influences of other abiotic factors on taxon richness patterns along altitudinal gradients. Method 2 performs better in controlling area effect than method 1 for seed plants on the Jade Dragon Snow Mountain, in terms of preserving the significant effects of abiotic factors. We not only strongly acknowledge the significance of method 1 in community ecology but also highlight the remarkable contribution of method 2 in eliminating area effects and delineating the explanatory power of abiotic factors. Therefore, we appeal for more application of method 2 in macroecological studies in mountainous regions.