Energy determines broad pattern of plant distribution in Western Himalaya

Abstract Several factors describe the broad pattern of diversity in plant species distribution. We explore these determinants of species richness in Western Himalayas using high‐resolution species data available for the area to energy, water, physiography and anthropogenic disturbance. The floral data involves 1279 species from 1178 spatial locations and 738 sample plots of a national database. We evaluated their correlation with 8‐environmental variables, selected on the basis of correlation coefficients and principal component loadings, using both linear (structural equation model) and nonlinear (generalised additive model) techniques. There were 645 genera and 176 families including 815 herbs, 213 shrubs, 190 trees, and 61 lianas. The nonlinear model explained the maximum deviance of 67.4% and showed the dominant contribution of climate on species richness with a 59% share. Energy variables (potential evapotranspiration and temperature seasonality) explained the deviance better than did water variables (aridity index and precipitation of the driest quarter). Temperature seasonality had the maximum impact on the species richness. The structural equation model confirmed the results of the nonlinear model but less efficiently. The mutual influences of the climatic variables were found to affect the predictions of the model significantly. To our knowledge, the 67.4% deviance found in the species richness pattern is one of the highest values reported in mountain studies. Broadly, climate described by water–energy dynamics provides the best explanation for the species richness pattern. Both modeling approaches supported the same conclusion that energy is the best predictor of species richness. The dry and cold conditions of the region account for the dominant contribution of energy on species richness.


| 10851
PANDA et Al. (Turner, Gatehouse, & Corey, 1987). In general, there is a proven high positive correlation between energy and species richness in cold climates. However, there are variations in these relationships in warm conditions (Francis & Currie, 2003;Kreft & Jetz, 2007). Water is an essential solvent for all physiological activities in plants and determines species patterns in the tropics, subtropics, and warm temperate zones (Hawkins et al., 2003). The scarcity of water adversely affects plant species richness in arid climates (Li et al., 2013). Water, energy, and their dynamics determine plant richness at high latitudes and explain more than 60% of the variations found in plant/animal species richness (Hawkins et al., 2003). However, the individual contributions of water and energy vary within regions. Vetaas and FerrerCastán (2008) observed that energy exerts greater control than does precipitation over the woody plants of the Iberian Peninsula. In contrast, Hawkins, Diniz-Filho, Jaramillo, and Soeller (2006) demonstrated a stronger influence of precipitation relative to energy in high to mid-latitudes.
Physiography describes the geographic complexity and regulates species diversity at both local and regional scales (Moeslund, Arge, Bøcher, Dalgaard, & Svenning, 2013). Several physiographic variables are used to quantify its impacts on species richness. Elevation regulates species richness by controlling the effects of climate and soil (Day & Monk, 1974), and it strongly influences the vegetation of most mountain ecosystems (Zhang et al., 2009). The terrain ruggedness index obtained by the differences between the elevation values of adjacent cells relative to a central cell (Riley, Hoppa, Greenberg, Tufts, & Geissler, 2000) describes the impacts of heterogeneity in species' niche differentiation (Whittaker, Levin, & Root, 1973). Aspect distributes the solar radiation affecting the microclimate and vegetation at a local scale (Kirkby et al., 1990). Slope correlates with the spatial pattern of tree species by controlling solar radiation (Bianchini, Garcia, Pimenta, & Torezan, 2010).
Disturbance is a critical factor for community composition and diversity in the species-rich landscapes of Western Himalaya (Kharkwal, 2009). Roy et al. (2002) used patch density, porosity, fragmentation, and juxtaposition to assess the disturbance regime of the region. Sharma, Gairola, Ghildiyal, and Suyal (2009) found a positive correlation between poor socio-ecological status of villagers and fuelwood collection in the temperate forests of the Garhwal Himalaya.
Overgrazing is reported to be another factor in the alpine grasslands of the Tibetan plateau (Yu et al., 2012). Both the human appropriations of net primary productivity and the global human footprint index, which assess the intensity of human intervention in ecosystems and its sustainability, have been used as socio-ecological indicators. Crowther et al. (2015) studied the effects of these variables to quantify human interference when mapping tree density at a global scale. Here, we emphasize their influence to monitor impacts of human disturbance on species richness. This has never been attempted previously in Western Himalayan studies.
Western Himalaya is geomorphologically complex with an altitudinal extent of >8,000 m. Its rich species diversity is broadly characterized by physiography, climate, soil, and anthropogenic disturbance (Shah et al., 2011). Its biogeography is climatologically distinct. It lies in a rain shadow region and receives little rainfall (Singh & Singh, 1987).
Western Himalaya gets rainfall during winter due to nonmonsoonal precipitation of westerlies. Overall, a temperate climate prevails.
However, the region experiences a broad range of temperature and precipitation anomalies with a mean annual temperature of ca. 5°C and annual precipitation of 2,500 mm (Chitale, 2014). The southern parts of the region are species-rich. For example, the angiosperm diversity of Himachal Pradesh comprises about 19,395 species or 7% of the world total (Karthikeyan, 2000). In contrast, the northern part is species poor due to the very low rainfall it receives (10-70 mm) and extremely cold temperature (≤45°C). The published literature focuses on the species diversity, community structure, and distribution pattern of its different forest types (Khan, Page, Hahmad, & Harper, 2013;Shaheen, Khan, Harper, Ullah, & Allem Qureshi, 2012;Sharma, Rana, Devi, Randhawa, & Kumar, 2014;Singh, 2008). Some studies have been carried out on species-environment relationships in adjacent areas (Oommen & Shanker, 2005;Wang, Tang, & Fang, 2007;Yan, Yang, & Tang, 2013). However, no specific study has been performed in Western Himalaya due to nonavailability of floral data. In this study, we took advantage of a scientifically designed national database to investigate the influence of the abiotic environment on species richness.
Understanding high species richness in a dry and cold climate of the understudied Himalaya remains interesting, and little research has focused on disentangling the effects of abiotic and human impacts. This study may provide better insights into ecologists and planners dealing with plant distribution patterns vis-à-vis climate change projections.
In this study, we explored the determinants of species richness in Western Himalayas. We used high-resolution species data available for the area to energy, water, physiography, and anthropogenic disturbance to explain distribution of species richness. Because of geomorphological complexity, we expected a greater control of physiography on the species richness. The general dry and cold conditions of the area may posit significant impacts of water and energy on species distribution. We speculate anthropogenic disturbance has some control on the plant richness because of the proximity of the region to densely populated areas. We applied predictive models to assess the influences of environmental variables and human impact on the species richness.

| METHODS
We selected two Indian Himalayan states, Himachal Pradesh and Uttarakhand (28.5-33°N and 75.5-82°E), for the study (Figure 1). We obtained plant species richness data from a national database developed during the project "Biodiversity Characterization at Landscape Level" (bis.iirs.gov.in). The field sampling involving all the forest classes was carried out using a stratified random sampling approach for the vascular plant species richness (Roy et al., 2012). A nested quadrat of 20 × 20 m was laid for trees or lianas, which accommodated two 5 × 5 m quadrats for shrubs/saplings and five 1 × 1 m quadrats for herbs/seedlings ( Figure 1). We gathered data pertaining to 27 variables under three categories, viz. climate (21 variables), physiography (four variables), and human disturbance (two variables) from different public domain sources (Appendix S1). We obtained 19 climate variables from the Worldclim site, along with two variables from CGIAR_ CSI (http://www.cgiar-csi.org/), all computed for the period from 1950 to 2000 (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005;Zomer, Antonio, Deborah, & Louis, 2008). We obtained physiographic variables (elevation, slope, and aspect) from the GMTED2010 database (https://lta.cr.usgs.gov). We derived the terrain ruggedness index from the differences of the elevation values of adjacent cells relative to a central cell. We procured the human appropriation of net primary productivity (Imhoff et al., 2004) and global human footprint (WCS-CIESIN, 2005) values from the SEDAC database (http://sedac.ciesin. columbia.edu/). The spatial resolution of each of the 27 variables was approximately 1 km 2 . The data conformed to WGS'84 projection. We extracted the data corresponding to 1,178 species locations using ArcGIS 10 for regression analysis (see Appendix S1).
Spatial autocorrelation describes the dependencies between observed samples and bias due to clustering, which can be quantified using Moran's I index (Moran, 1950). Principal component analysis (PCA) was used to select the best among correlated predictors (Xu, Wang, Rahbek, Sanders, & Fang, 2016). We computed Moran's I index of the floral data using the package "ape" in RStudio (Paradis, Claude, & Strimmer, 2004). We performed a multicollinearity test through hierarchal clustering analysis using the package "corrplot" in RStudio (Wei & Simko, 2016). We used PCA to examine the loadings of each collinear variable of the different categories separately and to pick a few least-correlated variables from each category for further analysis. We performed simultaneous centering and scaling to nullify the skewness effect. We performed PCA using the package "caret" in RStudio (Abdi & Williams, 2010). We calculated the percentage absolute weight of each variable corresponding to each principal component (PC) axis by multiplying the percentage variance explained by each axis with the absolute weight of each variable corresponding to that PC axis. We added the values of the first three PC axes, picked up the variable with the maximum weight, and removed its collinear partners (R > 0.8).
We continued this process and finally selected the top eight variables.
These include four climate (two energy and two water variables), two physiography and two disturbance variables. As two disturbance variables were noncollinear, we selected them directly without PCA. We examined the multicollinearity between the selected variables further to cross-check whether they satisfied correlation criteria (R < 0.8).
Several species distribution models (SDMs) deal with linear or nonlinear species-environment relationships. We selected the generalized additive model (GAM) for its robustness and proficiency in predicting species richness pattern (Hastie & Tibshirani, 1990). GAM fits smooth functions to establish nonlinear relations between species and the environment (Vetaas & FerrerCastán, 2008). The cubic spline criteria apply appropriate smoothness to model terms to minimize the GCV or GACV score, and, thereby, maximize the performance of the models. GAM quantifies complex species-environment relationships and reliable to fit local species occurrences (Aguirre- Gutiérrez et al., 2013). On the other hand, the structural equation model (SEM), a linear regression technique is often used to validate the output of another model such as GAM (Wu, Wurst, & Zhang, 2016). It establishes causal relationships between multiple predictors and response variables within a single framework (Grace, 2006).
We used GAM to find relationships between species richness and the eight selected variables both at independent and at cumulative levels with/without interaction(s). We fitted GAM with a Poisson error distribution with a "log" link function using the package "mgcv" in RStudio (Wood, 2016). We fitted the cubic spline (s) smoother to variables without interactive terms and the tensor product (te) smoother to interactive terms. We selected GCV (GCV.Cp) as prediction error criterion, which automatically selects models with the least error and improves the explanatory ability of models. The significance of each predictor was quantified using the deviance (%) [calculated as (null deviance-residual deviance/null deviance) × 100]. The data were split into four groups, and then, 1 group was held out and the model fit to the other three groups combined, that model then test on the held-out group to get the performance. This process was repeated, with each of the four groups acting as the held-out group, and overall R 2 taken as an average across all held-out groups. We used the package "lavaan" in RStudio for SEM to examine the performance of the GAM predictions (Rosseel, 2012). We derived structural equations for a different set of variables and categories and plotted them in figures wherein different colors were assigned to represent positive/negative correlation, and the thickness of arrows explained the impact of independent variables on the dependent variable ( Figure   3). The physiography, climate, and disturbance variables were grouped, and a separate SEM was fitted for each group. Further, we plotted the response curves of four climate variables to analyze the multidimensional regression through a two-dimensional representation ( Figure 4).

| RESULTS
The The first three PC axes of the climate variables explained 61.77%, 17.95%, and 9.08% variance, respectively. They cumulatively explained >91.8% variance, and the percentage absolute weight of each variable varies between 11.39% and 20.49%. We found at least three sets of collinear variables (Appendix S3). One set of variables was led by potential evapotranspiration (PET) with the maximum percentage absolute weight (20.48%), which showed a very strong correlation (R > 0.9) with the mean annual temperature, temperature of the coldest month, temperature of the coldest quarter, temperature of the driest quarter, temperature of the wettest month, temperature of the wettest quarter, and temperature of the warmest quarter, and a strong correlation (R > 0.8) with precipitation seasonality and temperature of the wettest quarter (Appendix S4).
The second set of variables was led by the aridity index, which showed a strong correlation (R > 0.8) with the mean annual precipitation, precipitation of the wettest month, precipitation of the wettest quarter, precipitation of the warmest quarter. Temperature seasonality led the next set of variables, which exhibited a very strong correlation with the temperature annual range and a strong correlation with the annual precipitation, precipitation of the warmest month, precipitation of the warmest quarter, precipitation of the wettest quarter, temperature of the coldest month, and temperature of the coldest quarter. were strongly correlated with each other (Appendix S4). However, the mean elevation was not considered for its strong correlation with selected climatic variables ( Figure S1; Appendix S4). The percentage absolute weights of slope and terrain ruggedness index were next to those of the mean elevation. They showed no strong correlation with the selected variables of other categories of variables and, therefore, selected for modeling. The eight least-correlated variables two each from energy, water, physiography, and disturbance: aridity index, human appropriation of net primary productivity, global human footprint, precipitation of the driest quarter, PET, slope, terrain ruggedness index, and temperature seasonality were finally selected for modeling ( Figure 2).
The study area was found to have a mean aridity of 0.13. The dryness of the northwest is greater than that of the southeast of the region. Structural equation model predicted negative relationships between the aridity index and temperature seasonality with species richness.
The precipitation of the driest quarter and PET were found to have positive correlations with species richness. The global human footprint and human appropriation of net primary productivity showed weak positive and weak negative correlations with species richness. Slope and terrain ruggedness index were found to have weak negative and weak positive correlations with species richness, respectively ( Figure 3). The response curves of GAM models showed trends similar to those of the SEM predictions. The piecewise polynomial curves were closely fitted with low degrees of freedom (7-9), that is, they neither straight nor wiggled. The aridity index had a linear and upward trend, but it became more or less straight after an index value of 0.10. Its variation from the mean zero level was found to be insignificant (Figure 4a). In contrast, the precipitation of the driest quarter exhibited an irregular curve pattern and differed greatly from the mean zero level. It showed a general positive trend with species richness (Figure 4b). PET showed a greater variation from the mean zero level and a positive relationship with species richness (Figure 4c). Temperature seasonality showed a complex curve pattern and was negatively correlated with species richness after approximately 625 CofV from the mean zero level (Figure 4d).

Structural equation model showed better goodness of fit with
interactions. The latent variables (groups of common predictors) were found to influence the species richness better compared with the combined effect of independent variables of the same groups (Table 1). The best model involving latent variables described 41% variance, the maximum predicted by any SEM model. The influences of physiography and disturbance with groups of common climatic predictors fitted better than their influences with separate groups of water and energy predictors ( Figure 3). GAM exhibited similar relationships, that is, enhanced performance with inclusion of interactive terms in the models. The performance of models improved with complexities.
The full model with a combination of all selected variables described the maximum deviance of 65.7% (Table 2). With climate, physiography showed greater significance than did disturbance. The deviance explained by common predictors of any two sets of variables with inclusion of interactive terms was between 52.8% and 62.1%, and the base climate model explained the deviance to the extent of 62.1%. In general, the inclusion of interactive terms of common predictors with noninteractive terms enhanced the percentage of the explained deviance.
With inclusion of interactive terms, the deviance explained by energy was 54.7%, a 5% greater value than the deviance explained by water.
Comparatively, the interactions between the disturbance variables explained more deviance than did the physiographic variables (Table 2).

| DISCUSSION
The dominant presence of herbaceous plant families, specifically Poaceae, Asteraceae, and Rosaceae, indicates the prevalence of grassland vegetation in the area. Similar familial relationships have also been reported by Rashid et al. (2012) in biodiversity-rich areas of the northern Himalaya with predominantly grassland vegetation.
(1) Water + Disturbance reported that Himalayan ecology is broadly defined by climate. Several studies describe strong relationships between climate and species richness (Clark et al., 1999;Curie et al., 2004;Hawkins et al., 2003;Ricklefs, 1987;Wang et al., 2007). Our results show the key roles of the water-energy dynamics in climatic control in explaining the capacity of species richness (O'Brien, 1993;O'Brien, Whittaker, & Field, 1998). The energy variables, that is, temperature seasonality and PET, are the most significant determinants of species richness patterns. The prevailing cold and dry conditions of the area might have led to the strong biotic dependency on energy. Many ecologists agree that energy has greater control in cold climates (Francis & Currie, 2003;Kreft & Jetz, 2007). The energy hypothesis proposes that greater availability of energy in a usable form enhances the species richness (Turner et al., 1987). The water variables show a strong association with species richness, but with a relatively lower magnitude compared with the energy variables. Water is an essential solvent for richness patterns in the tropics, subtropics, and warm temperate zones (Hawkins et al., 2003).
However, our results contradict with the findings of Hawkins et al. (2006), who demonstrated a stronger influence of precipitation compared with energy in high to mid-latitudes. The geographic position and evolutionary history of the area might have shaped this water and energy transitions with reference to species richness (Xu et al., 2016).
The species-energy affinities may be associated with Himalaya's youthful physiography and unstable geology (Mani, 1974). Although water and energy contribute differently to the species richness pattern, their synergy defines the mechanism by which the contemporary climate affects the species richness pattern. Water stress negatively affects species richness in spite of availability of sufficient energy, and with optimum water availability, plants would exploit the maximum photon flux essential for their physiological activities.
Water-energy variables and their interactions have primacy over anthropogenic disturbance and physiography. Independently, energy and water variables have primacy over nonclimatic variables.
The temperature seasonality is the greatest contributor, which indicates a significant influence of seasonal variation in temperature on the species richness. This variable exhibits complex, but negative correlations with species richness. This explains why the ecological stability of the region is balanced by fluctuations in temperature, and a global rise in temperature would be significant in controlling plant richness in future climates. PET is the second-best predictor of plant richness. It is a crucial factor in the Himalayan region (Chitale, 2014). PET has been reported to be the best predictor of the species richness pattern (O'Brien, 1993). A low level of PET shows positive relationships with the species richness. However, PET saturates at above 900 mm with no significant variations indicating the indirect influence of elevation on the species richness pattern. The SR, Species richness; PET, potential evapotranspiration; AI, aridity index; PDR, precipitation of the driest quarter; TS, temperature seasonality; SLP, slope; TRI, terrain ruggedness index; HNP, human appropriation of net primary productivity; HF, global human footprint. In SEM, some variables predicted insignificant (p > .05) are highlighted. HF and PET of M21 were significant at p < .05 and p < .01, respectively, for the same model; in GAM, HF of M20, M21, and M24 models was significant at p < .05, p < .05, and p < .01, respectively. "nr" indicates "no results"; The R 2 values the proportion of variance explained in the held-out group in the cross-validation.
T A B L E 1 Regression statistics for variable combinations without interactive terms using both structural equation model (SEM) and generalized additive model (GAM); results of GAM were computed using repeated cross-validation; except in few cases, the p value of each variable within the models was significant at p < .001. Both approaches supported the same conclusion that energy (temperature seasonality and potential evapotranspiration) was the best predictor of species richness, and the mutual influence of common predictors is more significant than their cumulative effects monotonous decrease in PET along elevation gradient probably explains this conjecture. The aridity index is the best-contributing water variable, which explains the impact of water stress on the species richness. This index was a crucial factor in an arid climate (Li et al., 2013). To a certain extent, dryness may facilitate species richness, but an index value >0.10 is probably counterproductive for plant richness. In general, a positive relation between the precipitation of the driest quarter and species richness is predicted. It is likely that species in places of high mean dryness face greater stress due to water deficiency during the driest quarter. In contrast, dry northwestern parts of the study area get the maximum precipitation during the driest quarter, which might have reduced the adverse effect of dryness on species richness. Anthropogenic disturbance has a key role in determining the species richness pattern of Western Himalaya (Gupta, 1978;Kharkwal, 2009;Negi et al., 2012;Shrestha et al., 2012). The two disturbance variables could explain the species richness almost equivalently. However, human appropriation of net primary productivity has a positive correlation with species richness, whereas global human footprint has a negative correlation.
Slope is reported to be a key factor in determining the vegetation at a local scale (Kirkby et al., 1990), but it may not be a significant factor at the landscape or regional scale.
The explanatory power of the models improved significantly with inclusion of interactive terms (

| CONCLUDING REMARKS
This study used a newly available national biodiversity database for India to explore species-environment relationships in two states in Western Himalaya on the basis of environmental variables and human impact. With GAM and SEM models including climatic, physiographic, and anthropogenic variables, energy variables, that is, temperature seasonality and PET, were most significant in explaining the species richness pattern. Water-energy variables and their interactions had primacy over anthropogenic disturbance and physiography. Disturbance is a critical factor in the southwestern region and places of low elevations with large populations. Physiography is less significant compared with disturbance. The exclusion of elevation from the list of selected variables may be significant in this low performance of the physiographic variables with respect to species richness. Further investigations may improve the findings if elevation is included in the list of predictors. Additionally, the better correlation of species richness with energy than with water needs to be verified in different climatic regions. Although this study focuses on disentangling the effects of abiotic predictors on species richness, the relationships between the biotic factors, that is, dispersion, competition, and isolation, are likely to improve the reliability and performance of the models. Nevertheless, the present study provides better insights into ecologists and planners dealing with plant distribution patterns vis-à-vis climate change projections.

ACKNOWLEDGEMENTS
The floral data utilized from Indian National Project on "Biodiversity Characterisation at Landscape Level" are thankfully acknowledged.
We owe our thanks to all four reviewers for their constructive suggestions to improve the quality of this manuscript.