The response of treeline-forming species to global climate change is uncertain. While numerous treeline species have recently experienced range advance along their upper elevational boundary, this has been species- and region-dependent. Making an accurate prediction of how taxa will respond is essential for conservation and land management, as treeline advance is likely to result in a loss of alpine biodiversity through habitat change and fragmentation. Predicting any species response requires an understanding of the current physical and climatic determinants of its distribution. We used the Maxent species distribution modeling software to predict the likelihood of treeline advance in the Nepalese Himalayas by modeling the extent of suitable habitats for 3 dominant treeline species—Abies spectabilis, Betula utilis, and Pinus wallichiana—under present and projected climate conditions. Temperature-related climatic variables and elevation explained the greatest amount of variance in the distribution of the study species. Under projected climate conditions, we found a regional increase in suitable habitat for all 3 treeline species, predicting a potential for northward and upslope advance.
The ranges of mountain plant communities are sensitive to climate change (Telwala et al 2013). As treeline ecotones invade high elevations, alpine tundra biomes shrink or shift upslope. Previous work has shown treelines are distributionally linked to temperature (Gehrig-Fasel et al 2007; Körner 2012), with new trees prevented from establishment above the upper treeline limit by unfavorable conditions. As climate has warmed, many treelines around the planet have shifted upward (Harsch et al 2009; Gaire et al 2014; Greenwood and Jump 2014). Nevertheless, not all sites have responded in the same way, with some species and regions remaining stable or actually retreating (Harsch et al 2009). Within regions already threatened by habitat change, estimating the degree to which future climatic shifts may influence habitat suitability for treeline-forming species is essential for both scientific understanding and conservation planning. Loss of alpine tundra to treeline advance (Moen et al 2004) reduces the habitat for endangered tundra species. For instance, the endangered snow leopard in the Himalayas could lose as much as 30% of its habitat under future warming scenarios and experience an increase in competition from other cat species (Forrest et al 2012; Lovari et al 2013). This combination of factors (loss of habitat and increased competition) has the potential to severely reduce Himalayan biodiversity.
The Himalayas contain an expansive area of understudied treeline habitat that is likely to respond to ongoing climate change. They also contain a steep elevational gradient, complex topography, high warming rate, biodiversity hotspots, and proximity to indigenous communities (Xu and Grumbine 2014). These features make the Himalayas a complex environment where treeline-related studies carried out elsewhere cannot be generalized. Therefore, studies focusing on the distribution of treeline species in the Himalayas are needed to understand and conserve Himalayan ecosystems (Singh et al 2013). Previous studies of treelines in this area have indicated mixed treeline response to climate change over a limited area (Gaire et al 2014; Shrestha et al 2014; Chhetri and Cairns 2015; Suwal et al 2016). The degree to which these trends scale up to broader patterns remains unresolved. Few modeling efforts have been conducted to identify contemporary climatic limits of Himalayan treeline-species ranges and predict distributions under potential future climate regimes (Schickhoff et al 2015). Recently, modeling studies of the environmental niches of Betula utilis in Uttarakhand, India (Singh et al 2013) and the Himalayan range (Schickhoff et al 2015; Bobrowski et al 2017) have been carried out, both indicating potential habitat shift. We wish to add to this knowledge by examining 3 representative Himalayan treeline-forming species from Nepal to gauge the future of this region.
In this study, we used species distribution modeling to determine if the distributions of these 3 species are defined by climate and predict whether their ranges are likely to expand or contract under projected climate conditions through an increase or decrease in suitable habitat. The study area is the country of Nepal, which covers the majority of the Himalayan range and has undergone significant shifts in climate within the last 50 years (Shrestha and Devakota 2010). The 3 study species, Abies spectabilis (D. Don), Betula utilis (D. Don), and Pinus wallichiana (A. B. Jacks), are important components of the Nepalese subalpine ecosystem and are dominant within the Nepalese Himalayan treeline (Chhetri et al 2017). We posed 2 questions:
Which topographic-climatic variables best explain the distributions of Himalayan treeline-forming species in Nepal? To answer this question, we examined the relationship between the contemporary distributions of our study species and the variables, both climatic (19) and topographic (5), that define each location.
Will the distribution of Nepalese treeline-forming species likely expand or contract under future climatic change? To explore this question, we projected how suitable habitats for the 3 study species would shift under 3 future climate scenarios and compared the elevational range of those habitats to their current range.
Material and methods
Nepal is a mountainous country and occupies the central part of the Himalayas (Figure 1). Analysis of climate records has indicated that temperature has increased at a higher rate in the Nepalese Himalayas than in other mountain areas, especially since the 1950s (Shrestha and Devkota 2010). This warming has progressed at a steady rate since the mid-1970s and is more pronounced at higher elevations. In addition, precipitation (rain, snow, sleet, or hail) has increased by 13 mm per year on average, while the number of rainy days has decreased by 0.8 days per year, suggesting more intense rainfall (Shrestha et al 2000). Several climate models have predicted a steady increase in temperature throughout Nepal accompanied by a decrease in monsoon rainfall in the north and an increase in the south (NCVST 2009; Shrestha and Aryal 2011).
Species selected for modeling
A. spectabilis, B. utilis, and P. wallichiana are dominant treeline species in the Nepalese Himalayas (Chhetri et al 2017). A. spectabilis (Himalayan silver fir) is a tall pyramidal evergreen tree growing in the subalpine forests of the Himalayas. It is found at 2800–4000 m above sea level (masl) in Nepal (Stainton 1972; Ghimire and Lekhak 2007). B. utilis (Himalayan birch) is native to the Himalayan region and is found at 2700–4500 masl. Blue pine (P. wallichiana) is an evergreen conifer found in the Himalayan region at 1800–4200 masl (Ghimire et al 2011).
Species occurrence data
We gathered species occurrence data from several sources: our field surveys carried out in different parts of Nepal (2004–2011), distribution records in the literature, the Global Biodiversity Information Facility ( www.gbif.org), and the Flora of Nepal database ( http://padme.rbge.org.uk/floraofnepal/index.php?page=home). We checked herbarium records for georeferencing error, misidentification, and duplication based on our own field surveys and information from published reports. We only included records collected after 1980 in our analysis to reduce potential error related to geographic coordinate accuracy. We used a total of 240 records for modeling: 94 for A. spectabilis, 85 for B. utilis, and 61 for P. wallichiana (Figure 1).
We gathered 19 bioclimatic variables with a 30 arc second (∼1 km) spatial resolution from the WorldClim–Global Climate Data ( www.worldclim.org) database to model the extent of suitable habitat of treeline tree species (Hijmans et al 2005) and their spatial climatic variation within the study area (Supplemental material, Table S1: http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1). Additionally, we determined the topographic characteristics of Nepal (eastness, elevation, northness, Topographic Position Index [TPI], and Solar Illumination Index [SII]) from the ∼1 km spatial resolution digital elevation model (DEM) obtained from Diva GIS ( www.diva-gis.org/) (Weiss 2001; Oke and Thompson 2015). Eastness (sine of aspect) and northness (cosine of aspect) are the linear component of aspect. TPI is a measure of surface undulation that allows an area to be classified both as to its topographic position (eg ridge top, valley bottom, mid-slope) and landform category (eg gentle valley, plain, steep narrow canyon, open slope). SII approximates the amount of direct solar radiation that hits an area as a function of its aspect, slope, and elevation.
We removed redundant climatic or topographic variables from our analysis that were highly correlated (Pearson correlation coefficient [r] greater than 0.9) with other variables. This resulted in 11 variables: isothermality (size of the day-to-night temperature oscillation relative to the summer-to-winter oscillation), annual temperature range, mean temperature of the coldest quarter, precipitation of the driest month, precipitation seasonality, precipitation of the warmest quarter, precipitation of the coldest quarter, eastness, northness, TPI, and SII (Supplemental material, Table S2: http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1). We retained elevation in spite of its high correlation with other variables because elevation has significantly improved the predictive ability of species distribution models (SDMs) for other high-elevation plant species (Oke and Thompson 2015) and is considered an important determinant of species distributions in mountain habitats (Körner 2012).
To determine the potential distribution of treeline tree species under projected climate conditions, we obtained downscaled WorldClim data from the Consultative Group on International Agriculture Research's (CGIAR)'s Climate Change, Agriculture and Food Security data archive ( www.ccafs-climate.org/). We used data generated by the Global Coupled Atmosphere-Ocean General Circulation Model developed by the Beijing Climate Center ( http://bcc.ncc-cma.net/). This is a regional model and performed relatively well over China and for nearby regions. It is based on the Intergovernmental Panel on Climate Change's fifth assessment report, which modeled greenhouse gas emission trajectories (called Representative Concentration Pathways [RCPs]). We used 3 trajectories (RCP 2.6, RCP 6.0, and RCP 8.5) and 2 time periods (2050 and 2070.) The first, RCP 2.6, is based on a reduction in greenhouse gas concentration. According to this lowest emissions scenario, global annual greenhouse gas emissions peak between 2010 and 2020, and temperatures are projected to increase in range from 0.3 to 1.7°C. The second, RCP 6.0, is based on a stable greenhouse gas concentration. Stable emissions assume that emissions peak around 2080 and then decline, with temperature increasing by 0.8–3.1°C. The last, RCP 8.5, proposes an increase in greenhouse gas concentration. Highest emissions assume greenhouse gas emissions will continue to increase throughout the 21st century, with temperature increasing from 1.4 to 4.8°C by 2100. Table S3 (Supplemental material, http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1) indicates climatic ranges and mean values in current and forecasted scenarios.
Species distribution model
We used the maximum entropy (Maxent) software version 3.3.3k (Phillips et al 2006) to conduct species distribution modeling. Maxent uses a machine learning process that estimates the suitable habitat for a species in an area based on the association between landscape variables and known distribution records (Kumar and Stohlgren 2009). The Maxent model works well with presence-only data (Elith et al 2011) and is capable of identifying contemporary and future habitat suitability under projected climate conditions (Hijmans and Graham 2006).
We ran our Maxent analysis using the following parameter values in our simulations: random test percentage, 30%; regularization multiplier, 1; maximum number of background points, 10,000; maximum iterations, 5000 or until convergence; convergence threshold, 0.00001; and cross-validated replicated run type. We set the model to remove duplicate presence records at the spatial resolution of the topographic-climatic variables (Soria-Auza et al 2010), so that we included only 1 presence record within the ∼1 km2 grid cell for each species. We ran 15 replicates for each species and averaged the results. The program created background data by using known occurrence points.
We used a Jackknife test to measure the performance of topographic-climatic variables in the model, reporting the importance in explaining the species' occurrence and the quantity of unique information each variable provided (Baldwin 2009). We used the omission and predicted area curve, and the area under the receiver operating characteristic curve, to assess the quality of the model (Fielding and Bell 1997). The final raster product for the study area contained pixel values between 1 (high habitat suitability) and 0 (low habitat suitability). We used a standard cutoff that defined pixels with values greater than 0.5 as being suitable for the species at that location (Kumar and Stohlgren 2009; Porfirio et al 2014). To determine which topographic-climatic variables best explain the distributions of Himalayan treeline-forming species in Nepal (question 1), we compared the contribution of each climatic and topographic variable in the model.
To determine if the suitable habitat for Nepalese treeline-forming species is likely to expand or contract under predicted future climatic shifts (question 2), we compared the elevation and area of the contemporary suitable habitat to that in the projected climate scenarios. To determine potential changes in the elevational distribution of treeline species, we extracted the elevation values of the pixels from the predicted and current distribution maps generated by the model using the DEM (Shrestha and Bawa 2014). To determine if the elevation for each species shifted significantly between the contemporary climate and the 3 projected climate scenarios, we used an independent sample t-test to compare the mean elevation value of the suitable habitats under each scenario (Shrestha and Bawa 2014).
We developed SDMs for 3 treeline species on the basis of 240 presence records. Omission and predicted area curves indicated moderate to good model fit for the 3 study species (Supplemental material, Figure S1: http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1). All of the observed omission rates on training samples and test samples were close to the predicted omission rates. The models for all 3 species performed better than random, as indicated by the average test AUC (area under the curve) values for the 15 replicate runs of A. spectabilis (0.89 ± 0.05), B. utilis (0.87 ± 0.05), and P. wallichiana (0.87 ± 0.06) (Supplemental material, Figure S2: http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1). The results indicated that current suitable habitat for A. spectabilis, B. utilis, and P. wallichiana cover 5130 km2, 9822 km2, and 9764 km2, respectively (Figure 2; Table 1). The model also predicted that the most suitable habitat for A. spectabilis and B. utilis was in the eastern and central part of Nepal, and the most suitable habitat for P. wallichiana was in the central part of Nepal.
Area of suitable habitat for the 3 treeline species under the current climate and 3 future climate scenarios. Numbers in parentheses indicate percent of Nepal's total land area; numbers outside parentheses represent blocks of about 100 km2.
A Jackknife test (Supplemental material, Figure S3: http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1) indicated that for all 3 species elevation is the most important predictor of habitat suitability. The climatic variables isothermality (Bio3) and mean temperature of the coldest quarter (Bio11) were consistently important predictors for all 3 species. The percentage of each variable's contribution to species distribution in the model is presented in Table 2. Elevation, isothermality, mean temperature of the coldest quarter, precipitation of the coldest quarter (Bio19), and topographic position index (TPI) had the most influence on the distribution of A. spectabilis. Elevation, mean temperature of the coldest quarter, isothermality, eastness, and precipitation of the warmest quarter (Bio18) had the most influence on the distribution of B. utilis. The distribution of P. wallichiana depended primarily on elevation, isothermality, topographic position index, mean temperature of the coldest quarter, and temperature annual range (Bio7). Response curves showing relationships between the top 3 topographic-climatic variables and species occurrence probability are presented in Figures S4–S6 (Supplemental material, http://dx.doi.org/10.1659/MRD-JOURNAL-D-17-00071.S1).
Contribution of different topographic-climatic variables to the distribution of the 3 study species.a)
The proportion of suitable habitat for the 3 study species varied under the projected climate scenarios (Figure 2). The model predicted that suitable habitat area for A. spectabilis will increase under all 3 scenarios in 2050 and 2070. The increase in area is projected to be greater in 2070 than in 2050, except under the RCP 6.0 scenario, where a dip in suitable habitat is predicted in 2050 (Table 1). The average elevation of A. spectabilis was higher under all warming scenarios (Table 3; Figure 3). The model predicted that the suitable habitat for B. utilis will increase in both area (Table 1) and elevation (Table 3; Figure 4) under all scenarios in 2050 and 2070. The model predicted that the suitable habitat for P. wallichiana will not significantly decrease or increase except under the RCP 6.0 scenario in 2070 (Table 1), but the average elevation of P. wallichiana was higher in all warming scenarios (Table 3; Figure 5).
Elevation (masl, average ± standard deviation) of suitable habitat for the 3 study species under the current climate and 3 future climate scenarios. We compared all future elevation estimates to current values using an independent sample t-test; all comparisons were significant (P < 0.001).
The distribution of suitable habitat for treeline species in the Nepalese Himalaya is likely to shift in response to climate change. The distribution of the 3 species examined here was mostly explained by climatic variables, indicating that climatic shifts will likely influence future habitat suitability. Explicit analysis confirmed this conclusion. Our simulations indicated that the suitable habitat for the 3 treeline species will shift toward higher elevations under predicted future climate conditions and may undergo an overall range expansion. These predicted shifts have serious consequences for both forest and tundra species in this area.
The contemporary SDMs performed well in describing current patterns when examined through observational and statistical analysis. All model AUC values were close to 0.90, which is considered a moderately good fit (Kramer-Schadt et al 2013) and comparable to values obtained in similar studies in this region (Rhododendron spp. AUC = 0.78, Kumar 2012; B. utilis AUC = 0.92, Schickhoff et al 2015). Our model results matched the actual existing range of these species-based fieldwork data and existing land cover maps and ecological maps.
The primary topographic-climatic variables describing the contemporary models for our 3 study species were meaningful given existing knowledge of the ecology of treeline species. Elevation was the most important predictive variable in the model for all 3 of our species. Elevation may not be directly associated with plant physiology, but it plays an important role in controlling atmospheric pressure, solar radiation, precipitation, and cloud cover (Oke and Thompson 2015), and thus it is a strong indicator of climatic variables that influence physiology. Mean temperature of the coldest quarter was one of the most important climatic variables explaining the contemporary distribution models for all 3 study species. This finding is consistent with ecological studies that have found that the upper elevational distribution limit of many treeline species is often determined by low temperatures (Körner 2012). Additionally, mean temperature of the coldest quarter and isothermality are related to growing season length, which has been identified as a secondary factor limiting treeline distribution (Körner 2012).
We should remain conservative in interpreting these findings, as we removed variables that were highly correlated with those included in our model, and simple spatial autocorrelation is an inherent effect in any climate study (Kumar 2012). Nevertheless, the strong climatic associations we observed do suggest, in the absence of other limitations, that ongoing climate change will influence the distribution of these treeline species.
Although the shift in suitable habitat distribution varied under the projected climate scenarios, future treeline advance is likely for the 3 study species. In terms of change in the area of suitable habitat, our models predicted minimal levels of expansion for A. spectabilis and P. wallichiana relative to B. utilis, which will double or triple under RCP 6.0 and RCP 8.5 scenarios. This is likely due to differences in species-specific requirements that define their distributions. P. wallichiana prefers dry valleys and foothills and is an early successional species (Ghimire et al 2011). In contrast, both A. spectabilis and B. utilis grow on moist north-facing slopes (Ghimire and Lekhak 2007). A. spectabilis has a broader distribution than B. utilis, often also occurring on south-facing slopes. B. utilis is often found in areas that receive snowmelt from mountain peaks above the treeline (Shrestha et al 2007).
Future distribution range will depend on how these species-specific microsites will change. General circulation models project warmer days and nights in the future, which will result in more snowmelt, which benefits B. utilis more than A. spectabilis and P. wallichiana. Western Nepal is drier than eastern Nepal, and monsoon precipitation is projected to increase in western Nepal, which will make western parts more suitable for tree establishment. Our model projected that B. utilis will occupy most of western Nepal under both RCP 6.0 and RCP 8.5 scenarios. Nevertheless, we saw average treeline extent advancing under all future climatic scenarios. Many climate models—including the Long Ashton Research Station Weather Generator (Agarwal et al 2014), Hadley Centre Coupled Model Version 3, and Regional Climate Model Version 3 (NCVST 2009)—project an increase in Himalayan winter temperatures, which will prolong the growing season, making the area above the treeline more favorable for treeline species.
These results are similar to those of previous studies examining species in this area. Schickhoff et al (2015) predicted that B. utilis habitat will shift northward throughout all of the Himalayas under future climate scenarios. Zomer et al (2014) predicted that subalpine conifer forest zones will shift upward by over 400 m between 2000 and 2050. On the other hand, Kumar (2012) predicted that the habitat of Rhododendron spp. will shrink under future climate scenarios in the Sikkim region of the Indian Himalayas, so not all species will necessarily follow this same path.
These results should be considered within the limitations of using a topographic-climatic approach to predict future species distributions (Macias-Fauria and Johnson 2013; Zong et al 2014). Limitations not addressed in this analysis may prevent these species from occupying the entirety of their fundamental niche in the future. These could include dispersal limitations of each species (Boisvert-Marsh et al 2014), which may limit their ability to establish in newly available suitable habitat. In addition, our projections do not factor in human disturbance (timber and fuelwood harvesting), herbivory (cf. Cairns and Moen 2004), plant physiology, soil type, snow cover, or land cover. Human disturbance and topographical factors can limit upslope treeline advance (Leonelli et al 2009).
Shrestha and Bawa (2014) emphasized the need for high-resolution environmental data that capture microclimates, edaphic conditions, vegetation dynamics, and landscape heterogeneity in SDMs. Kollas et al (2014) also emphasized the use of high-spatial-resolution temperature data for predictive modeling of temperature-based niche envelopes, recommending the use of topographic-climatic variables with a resolution of less than 100 m. They also suggested the use in modeling of absolute minimum temperatures instead of long-term means, because minimum temperatures determine the phenology of tree species at the cold limit. Building on higher-resolution datasets in future studies will improve the accuracy of these results. High-resolution climate data are not available for Nepal. Nevertheless, this study provides results based on the best available data, isolating climatic influences from other ecological limitations, and providing essential information for future environmental management.
Our model produced a good fit of contemporary species distributions, identifying suitable habitat for 3 dominant treeline-forming species of the Nepalese Himalayas under present and potential future climates. This approach could be usefully applied to other treeline species in the region. Our results indicate that the treeline ecotone is likely to transition throughout this region, which will likely have significant impacts on the associated plant and animal species. Our models predicted that an area above the existing treeline will become suitable for tree establishment. Nevertheless, this establishment will be controlled by factors like natural or human disturbances and ecological interactions with the surrounding shrub communities. Future work examining disturbance factors and species interaction, and incorporating high-resolution satellite imagery and a DEM, will improve the accuracy of this research.
TABLE S1 Topographic-climatic variables used in the Maxent model.
TABLE S2 Correlations among topographic-climatic variables.
TABLE S3 Range and mean value of key climatic variables in current and forecasted scenarios. Mean values are shown in parentheses.
FIGURE S1 Omission and predicted area curves for the 3 study species.
FIGURE S2 Area under the receiver-operating-characteristic curve (AUC) for the 3 study species.
FIGURE S3 Jackknife test of the relative importance of predictor topographic-climatic variables for the 3 study species. Blue bars indicate the importance of individual variables; red bars represent all topo-climatic variables.
FIGURE S4 A. spectabilis response to the top 3 topographic-climatic variables in order of strongest impact: elevation, isothermality, and mean temperature of coldest quarter.
FIGURE S5 B. utilis response to the top 3 topographic-climatic variables in order of strongest impact: elevation, mean temperature of coldest quarter, and isothermality.
FIGURE S6 P. wallichiana response to the top 3 topographic-climatic variables in order of strongest impact: elevation, isothermality, and topographic position index.
All found at DOI: 10.1659/MRD-JOURNAL-D-17-00071.S1 (676 KB PDF).