Distributional dynamics of a vulnerable species in response to past and future climate change: a window for conservation prospects

Introduction

Of the four billion species that have evolved on Earth over the past 3.5 billion years, 99% are considered to have disappeared (Novacek, 2001), most notably in the ‘Big Five’ mass extinctions (Raup & Sepkoski, 1982; Jablonski & Chaloner, 1994; Bambach, 2006). Considerable species losses over the past few centuries and millennia as a result of anthropogenic climate change have sounded the alarm about a possible sixth mass extinction (Barnosky et al., 2011). Effective conservation planning is necessary to avoid potential biodiversity loss, and an understanding of the distribution dynamics of organisms in response to climate change underlies the development of effective conservation strategies (Razgour et al., 2013).

Climate change is recognized as a key driver of species’ range dynamics, both locally and globally, over time (Huntley et al., 1995; Pearson & Dawson, 2003; Yates et al., 2010; Hamer et al., 2015). Quaternary climatic oscillations, particularly the most recent late Quaternary, characterized by markedly recurring glacial-interglacial cycles have played a crucial role in shaping the contemporary geographical distribution of plant species (Comes & Kadereit, 1998; Dynesius & Jansson, 2000; Hewitt, 2000; Sandel et al., 2011). The Last Interglacial (LIG, ∼120–140 ka) and Last Glacial Maximum (LGM, ∼21 ka) periods mark contrary extremes during the late Quaternary (Dawson, 1992), with the latter especially representing one of Earth’s most extreme periods of environmental variability (Clark et al., 2009). Indeed, the climate has warmed from the LGM to the present, and the temperature variations during this interval cover almost the entire temperature range of the Quaternary (Imbrie, McIntyre & Mix, 1989; Ruddiman, 2008). The dramatic climatic cooling during the LGM drove many species to glacial refugia (Nogués-Bravo et al., 2010), though populations began to recolonize during postglacial climate warming (Davis & Shaw, 2001; Normand et al., 2011).

As human activity intensifies, global temperatures are expected to rise by 1.1–6.4 °C during the 21st century (Intergovernmental Panel on Climate Change, 2007). Due to this unprecedentedly rapid rate of warming, climate change is predicted to be the greatest force in reshaping the geographical distribution of species in the 21st century (Leadley et al., 2010). Moreover, given the rapidity of climate change over the coming decades, whether populations can shift rapidly enough successfully tracking climate change is the central concern of many ecology studies (Davis & Shaw, 2001). Therefore, although many organisms have survived multiple climate cycles during their evolutionary histories (Meyers & Bull, 2002), some species are unable to disperse or adapt fast enough to track the rapidly changing climate, leading to increased extinction risk (Warren et al., 2001; Menéndez et al., 2006). In addition, landscape modifications resulting from the intensification of human activity may aggravate the negative effects of climate change by impeding species migration. Consequently, there is great concern about the challenges posed to extant species by the ongoing unprecedented change (Thuiller, 2007).

The response of species to climate change can be synthesized as evolutionary adaptation, dispersal or extinction (Parmesan, 2006; Aitken et al., 2008; Dawson et al., 2011), processes that are related to the velocity of the climate change and the species’ capacity to adapt and migrate (Jump & Penuelas, 2005; Sandel et al., 2011). Populations that are unable to keep up with climate change or to adapt to new climate conditions, especially those with narrow climatic tolerances (Thuiller, Lavorel & Araújo, 2005), face a very high risk of extinction (Hofreiter & Stewart, 2009; Dawson et al., 2011; Molinos et al., 2016). Endemic species are unique to a defined geographic range. These restricted-range species may be highly vulnerable to rapid climate change because of their narrow climatic tolerances (Malcolm et al., 2006; Ohlemüller et al., 2008), indicating that endemism and extinction risk are closely related (Petit, Hu & Dick, 2008). Therefore, surveying the response of endemic species to climate change is particularly important.

Pseudolarix amabilis is a representative of monotypic genus of the family Pinaceae; it is endemic to China, inhabiting a highly restricted area of the lower Yangtze River at an elevation range of 180–1,000 m (Yang & Christian, 2013). This deciduous tree’s branchlets are dimorphic: long branchlets (leading shoots) with helically borne leaves and short branchlets (brachioblasts) with fascicularly arranged leaves. The bract-scale complexes of the seed cones shed at maturity; two winged seeds, located at the base of the seed scales, mature in the 1st year (Fu, Li & Robert, 1999). P. amabilis is a wind-pollinated (Zanni & Ravazzi, 2007) and wind-dispersed species (Fordham & Spraker, 1997). The species grows on a variety of soils derived from acidic rock and is distributed in mixed-mesophytic and evergreen sclerophyllous broad-leaved forests (Wang, 1961; Farjon, 1990).

Currently, this species is ranked as at least ‘Vulnerable B2ab (iii, v) ver. 3.1’ and possibly the more threatened category ‘Endangered’ in the International Union for Conservation of Nature (IUCN) Red List of Threatened Species (Yang & Christian, 2013). Its population is severely fragmented, and the quality of the habitat and the population size continue to decline (Yang & Christian, 2013). To prevent deterioration, establishment of protected areas has been advocated by the IUCN (Yang & Christian, 2013). Compared to its current distribution, the extensive distribution of this species across the Northern Hemisphere was much wider during geological time (from the Cretaceous to the P1io-P1eistocene) (LePage & Basinger, 1995; Fig. S1). The sharp contraction suggests that the relict species P. amabilis has kept up with and survived past changes in climate. However, whether the ‘living fossil’ P. amabilis can cope with the challenge presented by the unprecedentedly rapid climate warming of the 21st century is a matter of concern. In addition, mast seeding—the synchronous intermittent production of large seed crops by populations (Kelly, 1994)—may aggravate the vulnerability of P. amabilis to climate change.

Reconstructing the range dynamics of species under past climate fluctuations is helpful for revealing their response mechanism to climate change. Clarification of a species’ response mechanism in combination with future distribution prediction can assist in developing effective conservation strategies (Petit, Hu & Dick, 2008; Désamoré et al., 2012). Ecological niche models (ENMs) link species occurrence records with environmental variables on multiple spatial and temporal scales to model past, present, and future distribution (Peterson et al., 2011). Such climate envelope models are important for revealing the past response of organisms to climate change, particularly for periods for which there is no fossil record. Overall, this tool has been regarded as an effective approach for reliably assessing the potential effects of future climate change on the fate of organisms (e.g., Huntley et al., 2004; Rodríguez-Sánchez & Arroyo, 2008; Elith & Leathwick, 2009).

Niche conservatism is a key assumption underlying the application of ENMs. In the case of P. amabilis, the evolutionary history indicates evolutionary stasis for this species since at least the Eocene (LePage & Basinger, 1995). As P. amabilis is an extraordinary example of evolutionary stasis, its niche is considered to be evolutionarily conservative because evolutionary stasis is perceived as related to niche conservatism (Eldredge et al., 2005; Stigall, 2012). The inference of niche conservatism for P. amabilis is supported by the general consistency of its temperature requirement throughout geological time (Bai & Li, 2017). In addition, this concern can be somewhat alleviated by the fact that a very general pattern of niche conservatism among species has been rather broadly confirmed in a recent review (Peterson, 2011).

In this study, we use an ENM to model the range dynamics of P. amabilis across four temporal frameworks representing four extreme moments of climatic variability during the Quaternary: the LIG, the LGM, the present and the future (2080). Our overall aim is to hind-cast the response of P. amabilis to chronological climate change and ultimately to provide some meaningful information about future conservation strategies for this vulnerable species based on the hind-casted response mechanism to climate change and the future distribution prediction.

Materials & Methods

Natural occurrence data for P. amabilis were compiled from the Chinese Virtual Herbarium (http://www.cvh.ac.cn/), the Global Biodiversity Information Facility (http://www.gbif.org/), and references (Ge et al., 1998; Hao et al., 2000; Duan et al., 2012; Huang et al., 2016).

Bioclimatic layers at a resolution of 2.5 arc minutes for the present (representative of 1960–1990), the LIG (∼120–140 ka) and the LGM (∼21 ka) climate were downloaded from the WorldClim dataset (available at http://www.worldclim.org/). Present climate data were generated through interpolation of average monthly climate data from global weather stations (Hijmans et al., 2005). LGM climate data were based on general circulation model (GCM) simulations from the Community Climate System Model (CCSM; Kiehl & Gent, 2004), and LIG climate data were based on models from Otto-Bliesner et al. (2006).

Future climate data were provided by the CGIAR Research Program on Climate Change, Agriculture and Food Security (CCAFS) (available at http://www.ccafs-climate.org). These climate data were generated by GCMs under a set of emissions scenarios. By the end of 2012, the closest emissions scenarios resulting from the Intergovernmental Panel on Climate Change (IPCC) process to the observed emission trends were the Special Report on Emissions Scenarios (SRES) A1B used in the IPCC Fourth Assessment Report and the Representative Concentration Pathways (RCPs) 8.5 used in the IPCC Fifth Assessment Report (Peters et al., 2012). The SRES A1B scenario describes a future of very rapid economic growth, a balance between the use of fossil fuels and non-fossil fuels and moderate human population growth (Intergovernmental Panel on Climate Change, 2007); the RCP 8.5 scenario depicts a world characterized by an atmospheric CO₂ concentration that continues to increase at current rates (Intergovernmental Panel on Climate Change, 2013). For this study, we chose four types of future (2080) climate layers generated from simulations using a set of GCMs under the above two emission scenarios: UK Meteorological Office (UKMO) Hadley Centre Coupled Model, version 3 (HadCM3) (Gordon et al., 2000; Pope et al., 2000), under SRES A1B; Met Office Hadley Centre (MOHC) Hadley Centre Global Environmental Model, version 2 (HadGEM2-ES (Earth System)) (Jones et al., 2011), under RCP 8.5; Canadian Centre for Climate Modelling and Analysis (CCCma) third-generation Coupled Global Climate Model with T63 spectral resolution (CGCM3.1-T63) (Flato, 2005), under SRES A1B; and Meteorological Research Institute (MRI) Coupled Global Climate Model, version 3 (CGCM3) (Yukimoto et al., 2012), under RCP 8.5.

The maximum entropy (MaxEnt) model, based upon the maximum entropy principle and using ‘presence-only’ species data (Phillips, Anderson & Schapire, 2006), has excellent predictive performance for threatened or range-restricted species (Elith et al., 2006; Hernandez et al., 2006; Hijmans & Graham, 2006; Pearson et al., 2007). We employed MaxEnt 3.4.1 for modelling. The responses of multivariate nonlinear models based on highly correlated climate variables can result in model overfitting (Morlini, 2006), which will also occur when such models are applied to data outside the training conditions (Graham, 2003; Morlini, 2006). To prevent multicollinearity and increase transferability effectiveness, we removed all environmental predictors with high pairwise correlation (Pearson’s correlation > 0.85, Table S1). Among highly correlated variables, we selected those with more direct physiologically roles in limiting the survival and reproduction of P. amabilis. We also deleted climate variables with a relatively low percentage contribution to the model performance. The contribution of each climate variable was assessed by the Jackknife procedure in MaxEnt (Fig. S2). These operations created a model of better performance with a balance between an underfitted model with few parameters and an overfitted model with too many correlated explanatory variables (Burnham & Anderson, 2002). The final set of selected bioclimatic variables included mean annual temperature (MAT, Bio1), temperature seasonality (TS; standard deviation of monthly mean temperature value, Bio4), min temperature of coldest month (MTCM, Bio6), mean temperature of driest quarter (MTDQ, Bio9), annual precipitation (AP, Bio12), precipitation of driest month (PDM, Bio14), precipitation of wettest quarter (PWeQ, Bio16) and precipitation of warmest quarter (PWaQ, Bio18).

ENMs assume that species are in equilibrium with their environments (Guisan & Zimmermann, 2000), i.e., species occur in all climatically suitable regions while being absent from all unsuitable ones (Hutchinson, 1957). In the process of investigating species occurrence data, however, some areas in a landscape were sampled more intensively than others; this sampling bias will lead to a lack of equilibrium of species distributions with climate. To reduce the effect of sampling bias, the approach of spatial filtering of a species occurrence dataset was applied in this study (Boria et al., 2014). We spatially filtered the distribution data of P. amabilis to obtain the maximum number of localities 10 km apart. As a result, 49 occurrence localities (Table S2) were rarefied to 47 points (Table S3). To minimize the effect of sampling bias, we also used a bias file representing a Gaussian kernel density of the species occurrence localities sampled within a 60-km search radius. The bias file upweights presence-only data points with fewer neighbours in a geographic space (Elith et al., 2011). The species occurrence data and climate layers were both projected to the Asia north equidistant conic projection in ArcMap.

We chose area under the receiver operating characteristic curves (AUCs) to assess the predictive performance of the ENMs (Fielding & Bell, 1997). The AUC statistic is bounded from 0.5 to 1.0, in which 0.5 indicates a random prediction (useless model) and 1.0 a perfect model prediction of presence versus absence; the closer AUC is to 1, the better the predictive accuracy of the model. Each model was run 10 times, with 75% of the species occurrence data selected for model training and 25% for model testing. To obtain binary predictions of the climate suitability of ENMs, MaxEnt’s logistic probability of occurrence output was converted to a binary mode (presence-absence output) using the maximum training sensitivity plus specificity logistic (MTSS) threshold and the 10 percentile training presence logistic (10% TP) threshold. By maximizing the proportions of actual positives and negatives that are correctly identified, prediction based on the MTSS threshold represents the most accurate forecast of presence/absence (Liu et al., 2005; Jiménez-Valverde & Lobo, 2007). The prediction based on the 10% TP threshold represents the core of species ranges by excluding the 10% of training localities with lowest prediction in the modelling (Morueta-Holme, Fløjgaard & Svenning, 2010; Anderson & Gonzalez, 2011).

To explore the distributional change between the two ENMs during each climate transition (e.g., from the LIG to the LGM, from the LGM to the current period, and from the present to 2080), we reduced the distribution to a single central point (known as a centroid) and created a vector depicting the magnitude and direction of the predicted change. Furthermore, to reveal potential migration in detail, we applied migration vector analysis. First, we calculated the geographic centroids of the ranges for every 60 × 60 km² for each period and then determined the centroid in the second time interval nearest to the centroid in the first time interval (i.e., LGM to LIG, present to LGM, and future to present); finally, we evaluated the potential migration from one period to the next.

To determine the variable with greatest influence for the model prediction in each grid cell, we performed limiting factor analysis. For any given point, the limiting factor is the variable with a value change at that point that results in the greatest change in the predicted probability of species occurrence (Elith, Kearney & Phillips, 2010; Elith et al., 2011). We focused on the limiting climatic factors affecting the contraction and expansion of climate suitability. The algorithm for the limiting factor analysis was programmed into MaxEnt, and the ‘density.tools.LimitingFactor’ command was used.

Results

The observed AUC values of 0.9802 (training data) and 0.9744 (testing data) indicate a good discrimination of the models between absence and presence cells, and the relatively high AUC values suggest that the species distributions were well predicted by climate. The MTSS and 10% TP thresholds were calculated to be 0.3048 and 0.3842, respectively.

It is noteworthy that the predicted area based on the MTSS threshold is very similar to that based on the 10% TP threshold (Figs. 1A–1D). The high similarity implies a consistence in the trend of distributional change over time between the predictions based on the two types of thresholds. Therefore, investigation of the predicted areas was based on binary ENMs generated by MTSS and 10% TP thresholds; however, assessment of the dynamic change of distribution over time was based only on the chosen MTSS threshold. Besides, compared to the present distribution of P. amabilis, its future distributions predicted using four types of climate layers all indicate an expansion northward and southern contraction, though the areas of the expansion and contraction differ (Fig. S3). Given the consistent trend of change and the focus of this study to explore the distributional change and potential migration routes of species, for conciseness, we choose one of the four future prediction for elaboration below: the prediction based on the HadCM3 simulation under SRES A1B.

The current potential distribution of P. amabilis involves three main disjunct districts (Fig. 1C): southeast China, where the actual population exists; and the southern frontier regions of Japan and the Korean peninsula, where extant populations are absent. The ENM projections for the other three periods (LIG, LGM and 2080) depict potential distributions with altered locations. The modelled potential climate suitability during the LIG shows four main disjunct regions, including southeast China, the southern frontier regions of the Himalayas and the southern frontier regions of Japan and the Korean Peninsula (Fig. 1A). During the LGM, fragmented distribution in southeast China is revealed (Fig. 1B), whereas the central-eastern regions of China, most parts of the Korean Peninsula and southern Japan are included in the future (Fig. 1D).

Contraction and expansion occur during each transition, and the resulting range varies significantly through time (Fig. 2A, Table 1). The ratio of the area between adjacent moments indicates a dramatic shrinkage from the LIG to the LGM, followed by significant expansion from the LGM to the current time and moderate expansion from the current period to the future (Fig. 2B, Table 1). The pattern of climate suitability change is dominated by contraction from the LIG to the LGM (Fig. 3A, Table 2), and by expansion since the LGM (Figs. 3C and 3E, Table 2).

With regard to distributional change, overall migration route analysis shows migration southeast from the LIG to the LGM (Fig. 3A), followed by migration northeast from the LGM to the present (Fig. 3C) and migration northward from the present to 2080 (Fig. 3E). From the LIG to the LGM, the expansion occurred mainly in the Sichuan Basin, with a source population from the most southeastern regions of China migrating northwest (Fig. 3B). From the LGM to the present, the expansion occurred in southeast China, with a relatively short migration distance (Fig. 3D). Northward expansion is the dominant tendency from the present to 2080 (Fig. 3F).

Limiting factors analysis indicates that the main climate variables influencing the model prediction for the change in climate suitability are temperature seasonality (TS), min temperature of coldest month (MTCM), mean temperature of driest quarter (MTDQ), annual precipitation (AP) and precipitation of driest month (PDM) (Figs. 4A–4C, Table 3, Fig. S4). The changes in the five variables set them within or outside the range of the physiological tolerances of P. amabilis (Fig. S5) and result in distribution expansion or contraction.

Discussion

Conclusions

In summary, in common with the general responses of species to climate change but with individual pattern, P. amabilis responds to glacial-interglacial cycles with high sensitivity, supporting the view that restricted-range species are sensitive to climate change (Sandel et al., 2011). The combination of investigating the response mechanism of P. amabilis to past climate change and predicting future climate suitability is beneficial for devising an effective conservation strategy. Our findings highlight the importance of combining a historical perspective with future predictions to develop a global conservation planning strategy for organisms in a changing world.

Uncertainties related to our model should be kept in mind when interpreting the results of this study. One important uncertainty derives from the fact that non-climatic factors, such as soil conditions, could not be integrated into our modelling because of the lack of sufficient data so far. In addition, many other ecological and evolutionary processes, such as biological interactions and interactions between the functional traits of an organism and its habitat, will also affect the distribution of species (Kearney & Porter, 2009; Fordham et al., 2012). These constraints are being addressed by some mechanistic modelling approaches. Compared to the MaxEnt-style correlative/statistical model, which statistically links spatial data to species distribution records, mechanistic models incorporate mechanistic links between the functional traits of organisms and their environments (e.g., Renton, Shackelford & Standish, 2012; Tomlinson et al., 2017). These two types of models have both strengths and weaknesses (Kearney & Porter, 2009; Kearney, Wintle & Porter, 2010). For poorly studied taxa with a paucity of knowledge about the physiological constraints on their survival and reproduction, such as P. amabilis, MaxEnt-style correlative/statistical models are a better choice. We will never be able to reconstruct the past and predict the future with accuracy, but we need a strategy for utilizing existing knowledge to reveal the likely effects of climate on species survival. Optimistically, we hold the opinion that as more data become available, ENMs will generate more realistic simulations and provide a solid basis on which to draw a more practical conservation strategy.

Supplemental Information