1. Introduction

The global mean temperature has increased by approximately 0.74 ℃ in the last one hundred years and this warming trend is projected to accelerate in the coming decades, increasing 2.8-5.3 ℃ by 2085 (Nogués-Bravo et al., 2007; Pacifici et al., 2015). Numerous studies have shown that in response to such warming trends the distribution of plant species tends to shift towards suitable habitat, specifically, upward in elevation and northward in latitude (Hickling et al., 2006; Chen et al., 2011). Given that the evolution of climatic niche partitioning tends to be much slower than rates of climate change (Quintero and Wiens, 2013), climate change-induced contractions in range sizes threaten many species with local extinction (Wiens, 2016). However, in response to similar climate changes, range expansion has been documented (Berry et al., 2003; Hamann and Wang, 2006). Thus, understanding how species ranges respond to climate change remains a major challenge.

Alpine species are considered especially sensitive to climate change (Lenoir et al., 2008). Due to their relatively low biotic complexity, climate change is suspected to be a dominant factor affecting range shifts in alpine species, which are often accompanied with marginal extinction at the warm edge and expansion at the cold edge (Parmesan, 2006; Pauli et al., 2007). In addition, because species distributed at high elevations are usually endemic to specific regions and have relatively narrow climate tolerances, theyoften experience a "nowhere-to-go" scenario, i.e., cold habitats are limited to the summits of mountainous regions (Nogués-Bravo et al., 2007). Models for Primula and Gentiana species predict decreased niches in response to climate change and grazing enforcement (Wehn and Johansen, 2015; Fois et al., 2016). Five comprehensive model projections have also shown that approximately 1/3 to 1/2 of the species distributed at higher elevations will lose more than 80% of their suitable habitats by 2070-2100, and cold-adapted mountain flora are at the highest risk for local extinction (Engler et al., 2011).

The Himalaya-Hengduan Mountains (HHM) region stretches over two global biodiversity hotspots, the Mountains of Southwest China and the Himalayas, both of which harbor a diversity of habitats and a high level of plant biodiversity (Mittermeier, 2004). A large proportion of alpine and sub-nival plant species are endemic to the region (Xu et al., 2014) and the area has experienced extensive glacial retreat and tree line advancement in the last century (Baker and Moseley, 2007). Two alpine genera in particular, Cyananthus and Primula, are suspected to be at risk of extinction due to climate change in the HHM (Zhou et al., 2013; Yan et al., 2015). The genus Cyananthus (Campanulaceae) comprises approximately 26 species with the highest level of diversity in the HHM (Zhou et al., 2013) and Primula (Primulaceae) contains about 500 species, among which 200 species are mainly distributed in this region (Yan et al., 2015). Most of the species in these two genera are perennial herbs, and although climate change is suspected to be a significant threat to remaining populations, there are few empirical studies that have investigated how these species may respond to climate change.

In this study, we explore how the distribution of selected Cyanathus and Primula alpine species endemic to the HHM may respond to climate change. Specifically, we extract species distribution information from herbarium specimens collected in the last century and, combined with climate data from WorldClim, we use species distribution models (SDM) to project species distributional dynamics under future climate warming conditions.

2. Materials and methods 2.1. Herbarium data

Species occurrence data were obtained from the Chinese Virtual Herbarium (CVH: http://www.cvh.ac.cn/). To assess possible changes in distribution for alpine species, we targeted Primula and Cyananthus species that had at least 100 specimen records distributed in the HHM: five species from each genus satisfied this criterion; their elevations ranged from ~2800 m to 4000 m. The original dataset contained 2, 866 specimen records, ranging from 123 specimens (Cyananthus lobatus) to 614 specimens (Primula sonchifolia). The following specimens were removed from the dataset: duplicate records with the same collector and collection number; specimens with elevation or collection records located far beyond the species distribution ranges based on the Flora of China; specimens with obvious identification errors; specimens that lacked spatial coordinates for georeferencing in Google Earth; and, to reduce the effect of spatial autocorrelation and the consequent overfitting, records within five kilometers of each other (Veloz, 2009). Our final dataset contained 708 records, which were used to conduct species distribution modeling (Table 1, Fig. 1).

2.2. Environmental data

Climate data were downloaded from WorldClim ver. 1.4 (http://www.worldclim.org/). We downloaded 19 bioclimatic variables (Table S1) for both current (Current: the average for 1960-1990) and future (2070: the average for 2060-2080) conditions at the highest available spatial resolution (30 arc-seconds; ~1 km). We chose the Global Circulation Model (GCM) ACCESS1-0, which has been shown to perform the best in the HHM region (Zhang et al., 2015; Wu et al., 2017). We used two greenhouse gas concentration trajectories for future climate conditions: an optimistic scenario whereby emissions peak around 2040 and then decline (Representative Concentration Pathway (RCP 4.5), and a pessimistic scenario whereby emissions continue to rise throughout the century (RCP 8.5).

We performed a Pearson's correlation test for the 19 climatic variables (Table S2) and the variables with correlation coefficients above 0.90 were removed to reduce the effect of multicollinearity. After this procedure, eight bioclimatic variables remained in our analyses, among which four variables were associated with temperature (bio1, bio2, bio3, bio4) and another four variables were associated with precipitation (bio12, bio14, bio15, bio18). Additional static environmental variables previously shown to contribute to model performance were used in the analysis (Stanton et al., 2012), including a land-cover layer (resolution of 10 arc-seconds) downloaded from the European Space Agency (ESA: http://due.esrin.esa.int/page_globcover.php), a soil layer (resolution of 30 arc-seconds) downloaded from the International Institute for Applied Systems Analysis (IIASA: http://webarchive.iiasa.ac.at/Research/LUC/External-World-soil-database/HTML/index.html?sb=1), and an elevation layer (resolution of 30 arc-seconds) downloaded from the Geospatial Information Authority of Japan (GIAJ: https://globalmaps.github.io/el.html). These variables were assumed to be constant under current and future scenarios. All layers were cut to restrict the projection area to the HHM, standardized to the same resolution as the climatic layers (30 arcseconds) with the same coordinate system (WGS1984), and transferred to ASCII format to be prepared for modeling operations.

2.3. MaxEnt modeling and output data analysis

We used MaxEnt to project species distribution under current and 2070 climate conditions (Phillips et al., 2006; Elith et al., 2011). Occurrence data and environmental layers for the current time period and the 2070 projection were used as input data in the model. We used the default settings for MaxEnt parameters and set the iteration time as 10 for each species, then we used the average logistic output layer (based on the ten replicates of cross-validation runs for each species) for subsequent calculations. Values for area under the receiver operating characteristic curves (AUC) were used to evaluate the performance of the models, where the value 0.5 indicates a random prediction and 1.0 a perfect model prediction (Bai et al., 2018). We reclassified the MaxEnt output file using the 10-percentile training presence logistic (10TPL) threshold value (Table 1) to define a species potential distribution region, above which species were considered "present" in the region, a method widely recognized for distinguishing suitable from unsuitable regions (Escalante et al., 2013; Kramer-Schadt et al., 2013; Radosavljevic et al., 2014; Hughes, 2017). We then calculated the longitude, latitude, elevation and range size for each pixel of potential presence and estimated the mean values for the current time period and 2070 for each species. The average value of longitude, latitude and elevation between the current time period and 2070 were compared within and among species. The shifts in total range size from the current time period to 2070 were also estimated. We then overlaid the distribution maps for the two time periods to detect potential changes in the occupied area.

All analyses were performed using MaxEnt 3.3.3k (Phillips et al., 2006), R ×64 3.3.3 (R Core Team, 2016), and ArcGIS10.2 (Environmental Systems Resource Institute (ESRI), 2014).

3. Results

All models in this study had AUC values > 0.978 (Table 1), indicating high performance. The mean 10TPL value was 0.229 (Table 1), which means that on average, we considered the cell with occurrence possibility above 0.229 to be suitable for species survival.

For all species, our models predict that, in response to climate change, distributions will shift upward in elevation, northward in latitude, and westward in longitude. Furthermore, these changes in distribution are predicted to be stronger in more pessimistic scenarios such as RCP 8.5. The mean shifts in elevation across the ten species were 532.49 m in the RCP 4.5 scenario (range: 348.11 m for Primula sikkimensis to 871.38 m for Primula poissonii) and 661.43 m in the RCP 8.5 scenario (range: 419.17 m for P. sikkimensis to 1042.38 m for P. poissonii). Mean shifts in latitude were 1.36° under the RCP 4.5 scenario (range: 0.82° for C. lobatus to 1.83° for P. sikkimensis) and 1.61° under the RCP 8.5 scenario (range: 1.03° for C. lobatus to 2.28° for P. sikkimensis) (Tables 2 and 3, Fig. 2).

In the RCP 4.5 scenario, our model predicts that from now until 2070 the range for all species will expand (range: 11.8% for Primula tangutica to 192.9% for P. sikkimensis); moreover, the newly occupied area will be larger than the area lost. Under the RCP 8.5 scenario, six of the 10 species showed greater range expansion than under the RCP 4.5 scenario (range: 75.5% for Primula secundifolia to 224.6% for P. sikkimensis), two species showed less range expansion (72.0% for Cyananthus macrocalyx and 11.4% for P. sonchifolia), while Cyananthus inflatus and P. tangutica were projected to lose area by 11.2% and 11.8%, respectively (Table 2, Fig. 2). Loss in distribution area was mainly confined to the southern portion of species distribution ranges, whereas gains in distribution area were mainly confined to the northwest on the Qinghai-Tibet Plateau.

4. Discussion 4.1. Changes in species distribution

Our models predict that the distribution of all Cyanathus and Primula species in our study will show a consistent response to a warming climate; specifically, the distributions of these species will shift upward in elevation and northward in latitude so that species will occupy living conditions similar to that of their previous habitats. Our predictions are consistent with previous studies that modeled distributional changes in response to climate change for other species in the region. For example, northward and upward shifts in distribution have been detected in wild soybean (He et al., 2016) and yews (Poudel et al., 2014). Furthermore, the westward shifts in distribution projected by our models are not surprising. Because high mountain habitats converge in the western portion of the HHM, when species distributions shift upwards in elevation, they will likely shift westward in longitude. This same trend has been projected for other alpine species in the HHM (Liang et al., 2018; He et al., unpublished).

4.2. Range shifts under future scenarios of climate change

All the species in our study are projected to expand their ranges under the more optimistic RCP 4.5 scenario but some of the species are projected to show greater range expansion under the more pessimistic RCP 8.5 scenario. These projections are in stark contrast to those of previous studies that have projected decreases in species ranges due to climate change. For example, the distribution of beech trees is expected to decrease by 2100 due to increased climate variability (Saltre et al., 2015), and by 2070 a decline in suitable climate space is projected for the two threatened dipterocarp trees (Deb et al., 2017) and the Tertiary relict tree species Davidia involucrate (Tang et al., 2017). Under the global warming scenarios projected in our study, however, suitable habitats may not be limited for most of the alpine species, which appear to be able to expand their distribution under our model projections. Similarly, a high proportion of plant species (63.6%) in the HHM have shown range expansion under future climate warming scenarios, but only a small proportion of plant species (11.9%) have been projected to experience range contraction (Liang et al., 2018). According to our model, the range of two species (C. inflatus and P. tangutica), which expand in the RCP 4.5 scenario, are projected to contract under the RCP 8.5 scenario. Under this same scenario, the ranges of another two species (C. macrocalyx and P. sonchifolia) are expected to expand less. Unsurprisingly, the species (P. tangutica) that occupies the highest latitude in our study is projected to encounter the greatest limitations to range expansion, likely due to the lack of suitable habitat in the north. While other three species (C. inflatus, C. macrocalyx and P. sonchifolia) that occupy relatively low latitudes and low elevations are projected to lose the largest areas in the southern margin of their distribution. Collectively, our projections indicate that although there is suitable habitat for future range expansion, if more extreme climatic changes persist, some of the alpine species in our study may decrease in range size, possibly resulting in local extinction.

We also found that the expansion area was in most cases larger than the contraction area, indicating that most of the species in our study will actually increase their area of occupied sites in response to climate change. Other models have projected that in response to a warming climate the land surface area suitable for plant occupation in some mountain habitats would increase rather than decrease (Elsen and Tingley, 2015). For example, Liang et al. (2018) concluded that, in the Hengduan Mountains, species would expand their ranges because the complex topography in high elevation areas provide sufficient land surfaces, which would promote not only upward shifts in elevation, but also expansion into the adjacent Qinghai-Tibet Plateau, which can accommodate northward shifting species and mitigate range contractions. Similarly, most Rhodiola species in the HHM are expected to expand their range sizes through upward shifts in elevation, as climate warming is expected to reduce permafrost and increase potential habitats for upshifting species (You et al., 2018). In this study, northwestern Yunnan, western Sichuan, and southeastern Tibet appear to be stable areas with suitable habitats that can accommodate alpine herbaceous species in response to climate change. These regions may be particularly useful for future in situ conservation efforts (Fois et al., 2017; Heywood, 2018), whereas immediate actions to minimize the effects of grazing and tourism may help to provide future habitats for species affected by climate change.

4.3. Future considerations

In this study, we assumed that potential habitats could be fully occupied by species with full dispersal capability when only considered the impacts of climate change. Our projections, however, did not consider several biotic factors that may impact our results, such as limited dispersal capability, pollination syndrome, and competition among species. For example, a loss of distribution area of up to 61.10% (full dispersal assumption) and 68.91% (no dispersal) was observed in SDM projections of crop wild relatives (Aguirre-Gutiérrez et al., 2017); species range sizes may also be controlled by mating systems, where self-pollinated species typically have larger geographic ranges than those of close relatives that are out-crossing (Grossenbacher et al., 2015). Although these sources of variation may affect the outcomes of our projections on a local scale, we anticipate that the general trends shown here may not change as climate is likely a dominant factor leading species distributional changes at a broader regional scale (Pearson and Dawson, 2003).

The model projections in this study span ten species across two genera and reveal almost the same distributional changes across all species. The consistency of these projections provides confidence in the validity of our models. However, although our modeling process was previously shown to perform well in the HHM region, our projections rely on this single GCM. To expand the confidence in model projections, future studies should include more GCMs to evaluate the magnitude of variation that may be associated with different climate models.

The species selected in our study are mainly those with an abundant number of herbarium specimens in order to have enough valid data for modeling projections. As such, the species included in our study are those with a relatively wide distribution and abundance. Species with a more localized distribution and less abundance are largely ignored in this analysis due to the limited number of herbarium specimens available for data collection. Such species, however, may be the most susceptible to climate change due to their narrow niche occupation and habitat requirements. This limitation highlights the need for future studies to include fieldbased investigations on the distribution of locally endemic species and their fate under future climate change scenarios.

5. Conclusion

Our study highlights the usefulness of herbarium specimen data for evaluating species distribution and predicting species distributional changes under the effects of climate change. The results are consistent with previous studies that show species will shift upward in elevation, northward in latitude and westward in longitude under a warming climate scenario, however, our species are also projected to expand their ranges as the loss of previously occupied habitat is less than that of newly acquired areas. Great opportunities for the exploration of species distributional dynamics under climate change remain to be explored, including using digitalized herbarium specimens, citizen scientists, and targeted field surveys. These new approaches to examining responses to climate change may provide insights for targeted species protection. The utilization of more comprehensive models that consider both the biotic and abiotic factors contributing to species persistence will provide valuable data for making conservation decisions under future climate scenarios.

This study was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31010000), the National Key Basic Research Program of China (2014CB954100), and the Program of Science and Technology Talents Training of Yunnan Province (2017HA014). We appreciate the Chinese Virtual Herbarium for providing specimen data of Cyananthus and Primula species.

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2019.01.004.