www.nature.com/scientificreports
OPEN
received: 22 August 2016 accepted: 01 February 2017 Published: 08 March 2017
Potential effects of climate change on geographic distribution of the Tertiary relict tree species Davidia involucrata in China Cindy Q. Tang1, Yi-Fei Dong1, Sonia Herrando-Moraira2, Tetsuya Matsui3, Haruka Ohashi3, Long-Yuan He4, Katsuhiro Nakao5, Nobuyuki Tanaka6, Mizuki Tomita7, Xiao-Shuang Li8, Hai-Zhong Yan1, Ming-Chun Peng1, Jun Hu9, Ruo-Han Yang1, Wang-Jun Li1, Kai Yan10, Xiuli Hou11, Zhi-Ying Zhang1 & Jordi López-Pujol2 This study, using species distribution modeling (involving a new approach that allows for uncertainty), predicts the distribution of climatically suitable areas prevailing during the mid-Holocene, the Last Glacial Maximum (LGM), and at present, and estimates the potential formation of new habitats in 2070 of the endangered and rare Tertiary relict tree Davidia involucrata Baill. The results regarding the mid-Holocene and the LGM demonstrate that south-central and southwestern China have been longterm stable refugia, and that the current distribution is limited to the prehistoric refugia. Given future distribution under six possible climate scenarios, only some parts of the current range of D. involucrata in the mid-high mountains of south-central and southwestern China would be maintained, while some shift west into higher mountains would occur. Our results show that the predicted suitable area offering high probability (0.5‒1) accounts for an average of only 29.2% among the models predicted for the future (2070), making D. involucrata highly vulnerable. We assess and propose priority protected areas in light of climate change. The information provided will also be relevant in planning conservation of other paleoendemic species having ecological traits and distribution ranges comparable to those of D. involucrata. Since the maximum extension of ice sheets during LGM (Last Glacial Maximum, from 26,000 to 19,000 calendar years before the present, with cold conditions also prevailing until 14,500 yrs ago), the global climate has undergone rapid and remarkable changes toward a generally warmer current level1. To reduce the rate of species extinctions in a world dominated increasingly by human, natural protected areas are often the strategy on which conservation measures are built. However, in light of global climate change, conservation strategies and decisions regarding the location of protected areas must consider future regional climate changes and their effects on species distribution, including retractions and expansions of species ranges2–4. Determination of the extent of threatened plants’ response to climate change can be useful in formulating flexible conservation strategies for China5. 1 Institute of Ecology and Geobotany, Yunnan University, Kunming 650091, China. 2Botanic Institute of Barcelona (IBB-CSIC-ICUB), Passeig del Migdia s/n, Barcelona 08038, Spain. 3Center for International Partnerships and Research on Climate Change, Forestry and Forest Products Research Institute, Matsunosato 1, Tsukuba-shi, Ibarakiken, 305-8687, Japan. 4Kunming Institute of Forestry Exploration and Design, the State Forestry Administration of China, Kunming 650216, China. 5Kansai Research Center, Forestry and Forest Products Research Institute, Momoyama, Kyoto, 612-0855 Japan. 6Tokyo University of Agriculture, 1-1-1 Sakuragaoka, Setagaya-ku, Tokyo 1568502, Japan. 7Tokyo University of Information Sciences, 4-1 Onaridai Wakaba-ku, Chiba 265-8501, Japan. 8Yunnan Academy of Forestry, Kunming 650204, China. 9Key Laboratory of Mountain Ecological Restoration and Bioresource Utilization and Ecological Restoration Biodiversity Conservation, Key Laboratory of Sichuan Province, Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu 610041, China. 10Center for Mountain Ecosystem Studies, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming 650201, China. 11Department of Life Science and Technology, Kunming University, Kunming 650214, China. Correspondence and requests for materials should be addressed to C.Q.T. (email:
[email protected]) or M.C.P. (email:
[email protected])
Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
1
www.nature.com/scientificreports/ Identifying past, contemporary and future climate refugia of Tertiary relict plants can clarify their conservation significance. China enjoys great species diversity, with ca. 30,250 seed plant species. The key areas of south-central China and Indo-Burma (including most parts of Yunnan) are two of the 25 designated global biodiversity hotspots6. Two plant diversity centers with high levels of endemism are located in this region7. The area was never covered by extensive, unified ice-sheets, owing in part to its complex topography8 and to moderate cooling as compared to other regions of China9. Moreover, these areas did not suffer the extreme aridity that is usually associated with the glacial periods: within China, south-central and southwestern areas were those with the smallest LGM precipitation deficits, rainfall being even higher than at present during the autumn and winter seasons9. So the mountains of south-central and southwestern China have had a relatively stable long-term environment, and are now referred to as the Pleistocene glacial refugia10–12. Before and during the Ice Ages, a number of famous gymnosperms such as Ginkgo, Metasequoia and Cathaya disappeared from the rest of the Northern Hemisphere, and coincidentally, some temperate deciduous broad-leaved tree taxa died off in North America (e.g., Davidia13,14; Tetracentron15; Cercidiphyllum in review16), Europe (e.g., Tetracentron17; Cercidiphyllum in review16), the Russian Far East (e.g., Tetracentron in Kamchatka18), and Japan (e.g., Davidia19–21; Tetracentron22). Most of them now only survive in the humid subtropical areas of south-central and southwestern China. A number of these paleoendemics have highly restricted distributions, usually in the form of small, localized populations23. Refugia represent climatically stable areas and remain a high conservation priority as important sites for the long-term persistence of species. Today, they are in danger of overexploitation. Climate change has also become a great threat to these species and their forests. Davidia involucrata Baill. (Figure 1) (dove tree or handkerchief tree in English, gongtong or gezishu in Chinese), included in Davidiaceae (by some in Nyssaceae or Cornaceae), is a rare and endangered tree species with small populations restricted to the humid mountains of south-central and southwestern China, mostly confined to the mountains surrounding the Sichuan Basin. The forests containing D. involucrata as one of the dominants are mainly scattered at 1300–1900 m altitude (but at 2300–2800 m in Yunnan) in areas with cool, foggy and cloudy climatic conditions. At these altitudes, D. involucrata co-dominates with other Tertiary relict plants including Tetracentron, Cercidiphyllum, Tapiscia, Dipteronia, Pterostyrax, Carya, Liquidambar, and Decaisnea (Fig. 2). Davidia involucrata also occurs, but very rarely, between 900–1300 m, growing with some evergreen broad-leaved trees of Cyclobalanopsis, Castanopsis, Lithocarpus, Phoebe, Machilus and Neolitsea. No adequate studies have been carried out to estimate the potential impact of climate change on this Tertiary relict species, and so far no conservation strategies that consider climate change have been proposed. Species distribution modeling (SDM) combines species occurrence data with environmental variables, under the assumption that the distribution of known localities reflects survival patterns of species24. SDM has proven useful for estimating such patterns and the resultant risk of extinction. It can produce spatially explicit and comprehensive maps that are particularly valuable for identifying areas where conservation efforts and management strategies are most needed. Our main objective is to apply an SDM approach based on Tertiary relict tree species D. involucrata’s occurrence data and climatic surfaces at four different time points (LGM, mid-Holocene, present, and the year 2070), so as to estimate past distributions of climatically relevant areas, model the present potential distribution range, and predict potential distribution and vulnerability under future climate change. We hypothesize that the mountains of south-central and southwestern China have been long-term stable refugia of D. involucrata and that there would be a decline in areal extent of suitable habitats for this tree species under future climate change. In order to provide a basis for conservation management of D. involucrata, our aim is to map the past, current and future suitable habitat distributions of this paleoendemic species, and propose establishment of protected areas.
Methods
Study species and areas. Davidia involucrata is a deciduous broad-leaved tree and its fruit stone is heavy and woody, composed of fibers, each nut containing 3–6 seeds. It is a protected species in China, already included in the National List of Rare and Endangered Plant Species of 1984 (listed as “first grade” nationally protected) and later in the Catalogue of the National Protected Key Wild Plants of 1999 (also as “first grade”). Its relative rarity has also been recognized: the species was already included in the first red book of China (as a “rare” species)25, although in the red list of 2013—which follows the internationally recognized IUCN categories—the species is listed, surprisingly, as LC (“Least Concern”), that is, with no conservation concerns26. Davidia involucrata is scattered in isolated mountain slopes, or in valleys where the soil often contains much gravel, or by streams, or scree slopes harboring unique assemblages of plants27. Today scattered stands of D. involucrata range approximately from 98–110°E, 26–32°N in south-central and southwestern China, including Hunan, Hubei, Chongqing, Sichuan, southern Shaanxi, southern Gansu, Guizhou and Yunnan Provinces. The latitudes are subtropical. Vegetation types in this area are diverse, according to altitudes, mainly including subtropical evergreen broad-leaved forests, warm temperate deciduous broad-leaved forests, deciduous and evergreen broad-leaved mixed forests, warm temperate coniferous forests, coniferous and broad-leaved mixed forests, cold-temperate coniferous forests, alpine scrub and meadow; however, by far most of the vegetation consists of what is commonly referred to as the subtropical evergreen broad-leaved forest. A striking feature is that within the transitional altitudinal zone from subtropical warm and humid evergreen to the temperate deciduous broad-leaved forests, where in some cases there is also an admixture of species of gymnosperms, a number of extant Tertiary relict plant species thrive today23. The climate is dominated by the Asian monsoon system, including the East Asian summer monsoon, the Indian summer monsoon, and the East Asian winter monsoon, with dry continental winds in winter and moist oceanic winds in summer. In general, in the subtropical areas, lands east of longitude 103° are more influenced by the East Asian summer monsoon, while to the west the Indian summer monsoon dominates. Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
2
www.nature.com/scientificreports/
Figure 1. (A) An inflorescence with two showy white bracts of D. involucrata (Photograph by Shi-Liu Wang). (B) Fresh fruits with autumn foliage D. involucrata (Photograph by Cindy Q. Tang). (C) A forest of D. involucrata at ca. 1480 m in Longcanggou, Sichuan (Photograph by Jun Hu).
Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
3
www.nature.com/scientificreports/
A Heilongjiang Neimenggu Gansu
Ningxia
nxi
Xinjiang
Hebei
Sha
Qinghai
Shandong
Shaanxi Henan
Tibet
Jilin
Liaoning
Sichuan
Guizhou
Jiangsu
Anhui
Ch Hubei
Zhejiang Hunan Jiangxi
Fujian Guangxi Guangdong
Yunnan
Taiwan
B
33
Hainan
32
Wenxian 1900
Davidia-Betula
1300
Lithocarpus-Davidia-Cornus
Dabashan
1800 1400
Lithocarpus-DavidiaTetracentron
Beichuan
31
Wolong
1900 1400
Baiguolinchang
Davidia-Carya-Fagus
Changyang Xingdoushan 1800 1100 Cyclobalanopsis-Davidia 900 1500
Davidia-Juglans Davidia-Machilus-Phoebe
Houhe
Emeishan
30
Davidia-TetracentronCercidiphyllum
Davidia-Neolitsea
1570
Labahe 1800 1600
1680
1450 1300
Tetracentron-Tapiscia-Davidia
1700 1600
Fagus-Davidia Dipteronia-Davidia DavidiaTetracentron-Cornus
Davidia-CercidiphyllumTetracentron
Badagongshan Davidia-TetracentronFagus
1200
Davidia-Cyclobalanopsis
29
1600
2400
Weixi
Davidia-Acer-Tetracentron
Yiliang 1860
DavidiaPterocarya
1950
DavidiaCarya-Padus
Picea-Davidia-Populus
28
2400
Gongshan
2680
Baizhishan
Sanjiangkou
Fanjingshan
Davidia-Acer-Tetracentron
Beilidujuan Yaoshan
Davidia-CercidiphyllumPterocarya
2400
2600 2350
1500 1300
Davidia-StyraxLiquidambar-Pterostyrax
DavidiaPrunus-Litsea
27
Nayonglinchang Davidia-DecaisneaDipentodon
26
1930
98
99 102
103
104
105
106
107
108
109
110
111
Figure 2. (A) Provinces of China. Ch = Chongqing. (B) The spatial distribution pattern of forests containing D. involucrata as one of the dominants. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v. 9.3.
Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
4
www.nature.com/scientificreports/
Figure 3. Schema of species distribution modeling in Davidia involucrata, taking into account the uncertainty derived from the evaluation of model performance. (A) A “standard” model for the present time period is obtained by averaging the multiple runs (20 in our case) from the cross-validation (or from another sampling technique, e.g. subsample); (B) a map showing the variability in the predicted suitable areas among all generated models (20) from the cross-validation (“uncertainty” map); (C) A “refined” map is obtained by removing the pixels that are forecast by 10 were excluded to avoid multi-collinearity. VIF were calculated using vif function of usdm package40 in R. The best combination among 6236 candidate parameter combinations was chosen following corrected Akaike Information Criterion (AICc41), and consisted of six variables: annual mean temperature (bio1), isothermality (bio3), temperature seasonality (bio4), precipitation seasonality (bio15), precipitation of the warmest quarter (bio18), and precipitation of the coldest quarter (bio19). The “best” model was evaluated by 20-fold cross-validation with MaxEnt, with performance assessed using the area under the curve (AUC) of the receiver operating characteristic (ROC) plot. AUC scores may range between 0.5 (randomness) and 1 (exact match), with those above 0.9 indicating a good performance of the model42. The MaxEnt jackknife analysis was used to evaluate the relative importance of the six bioclimatic variables employed, based on their gain values when used in isolation. In k-fold cross-validation, occurrence data is divided randomly into k equal-size groups, and models are built using k −1 bins for calibration in each iteration, with the left-out bin used for evaluation; background data are sampled by Maxent from the entire study region43. We applied the maximum sensitivity plus specificity (MSS) logistic threshold, which is very robust with all types of data44, to obtain a map of absence/presence (with the probability of presence shown as continuous values from the
Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
7
www.nature.com/scientificreports/
Figure 5. A comparison of potential habitats under the present climate and three climatic scenarios in the mid-Holocene. (A) Potential habitats under the present climate; (B) Potential habitats under the climatic scenario mid-Holocene-CCSM; (C) Overlap areas of the mid-Holocene-CCSM and the present. (D) Potential habitats under the climatic scenario mid-Holocene-MIROC; (E) Overlap areas of the mid-Holocene-MIROC and the present. (F) Potential habitats under the climatic scenario mid-Holocene-MPI. (G) Overlap areas of the mid-Holocene-MPI and the present. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v. 10.2.
threshold to 1). The resulting map, which is the average of the multiple runs from the cross-validation, can be regarded as the “standard” model (Fig. 3A). As suggested by many authors, uncertainty should be addressed in Ecological Niche Modeling (ENM) because it can produce large biases on the niche predictions45 and references therein. Thus, although our model performed well (AUC = 0.982 ± 0.006), a “consensus ensemble approach” or “consensus approach”46 was applied in order to Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
8
www.nature.com/scientificreports/
Figure 6. A comparison of potential habitats under the present climate and three climatic scenarios in the LGM. (A) Potential habitats under the present climate; (B) Potential habitats under the climatic scenario LGMCCSM; (C) Overlap areas of the LGM-CCSM and the present. (D) Potential habitats under the climatic scenario LGM-MIROC; (E) Overlap areas of the LGM-MIROC and the present. (F) Potential habitats under the climatic scenario LGM-MPI. (G) Overlap areas of the LGM-MPI and the present. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v. 10.2.
visualize the variability in the predicted suitable areas among all generated models. A map showing such variability (or uncertainty) was obtained as follows (Fig. 3B): first all 20 continuous projection maps were converted into binary maps (applying the MSS threshold), and then were calculated how many models predict each pixel as a suitable area. Thus, a continuous map with pixels ranging from 0 to 20 was obtained. Converting these values to probabilities (i.e., 0 as 0%, 10 as 50%, or 20 as 100%), we finally obtained an estimate of each pixel’s probability of being predicted as suitable area. Every pixel for which at least 95% of the models (i.e., 19 of 20 and 20 of 20) Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
9
www.nature.com/scientificreports/
Figure 7. A comparison of potential habitats under the present climate and three climatic scenarios in future (2070). (A) Potential habitats under the present climate; (B) Potential habitats under the climatic scenario 2070-CCSM RCP 2.6; (C) Overlap areas of the 2070-CCSM RCP 2.6 and the present; (D) Potential habitats under the climatic scenario 2070-GFDL RCP 2.6; (E) Overlap areas of the 2070-GFDL RCP 2.6 and the present; (F) Potential habitats under the climatic scenario 2070-MPI RCP 2.6; (G) Overlap areas of the 2070MPI RCP 2.6 and the present. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v.10.2.
forecast species presence was regarded as a pixel of high probability of being predicted by modeling (“presence pixel”). In an attempt to reduce the uncertainty in our models, those pixels not ranking as “presence pixels” were removed from the “standard” output maps to obtain a “refined” model where only those predicted areas with 95% of confidence are shown (Fig. 3C), and this procedure was applied for the present time and for all past and future scenarios (Supplementary Fig. S1). Although the areas predicted as suitable are somewhat reduced (for Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
10
www.nature.com/scientificreports/
Figure 8. A comparison of potential habitats under the present climate and three climatic scenarios in future (2070). (A) Potential habitats under the present climate; (B) Potential habitats under the climatic scenario 2070-CCSM RCP 8.5; (C) Overlap areas of the 2070-CCSM RCP 8.5 and the present; (D) Potential habitats under the climatic scenario 2070-GFDL RCP 8.5; (E) Overlap areas of the 2070-GFDL RCP 8.5 and the present; (F) Potential habitats under the climatic scenario 2070-MPI RCP 8.5; (G) Overlap areas of the 2070MPI RCP 8.5 and the present. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v. 10.2. example, about 37% for the present time in D. involucrata; Fig. 3 and Supplementary Table S2) with this method, the robustness of the forecast is improved46. A proof of this is that the reduction of suitable area when the uncertainty is taken into account is mainly focused on those pixels with low logistic probability (threshold to 0.5, area loss = 46.54% for the present time; 0.5 to 1, area loss = 3.99 for the present time; Supplementary Table S2). All ENM predictions were visualized in ArcGIS v. 10.2 (ESRI, Redlands, CA, USA). The suitable area (in km2) for all models at each time slice for D. involucrata was also calculated in ArcGIS. To estimate suitable area gains Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
11
www.nature.com/scientificreports/
Figure 9. (A) Current presence points of D. involucrata within and outside existing protected areas (nature reserves) of China. (B) Modeled potential habitats with and without protection under the present climate in China, with areas proposed for conservation. (C) Predicted potential habitats with and without protection under a future 2070 climate in China, with areas proposed for conservation, exemplified by the scenario 2070CCSM RCP 2.6. Purple lines: national boundaries between China and India (at issue). Maps were generated using the software ArcGIS v. 10.2.
Scientific Reports | 7:43822 | DOI: 10.1038/srep43822
12
www.nature.com/scientificreports/
Figure 10. Area of potential habitats of D. involucrata with and without protection in China.
or losses (or unchanged areas) for both past and future scenarios with respect to the present, binary output maps were overlapped with the Intersect Tool of ArcGIS. To provide conservation advice on D. involucrata, the digitized map of Chinese protected areas was overlapped with the species occurrences and both binary maps of the present and future scenarios. The climate data under various scenarios for the three time periods (present, mid-Holocene, LGM) are provided in the Supplementary Tables S3‒S4.
Results
Current distribution, model performance, and potential distribution under the present climate.
The present distribution area of D. involucrata is confined to mountains with complex topographies surrounding the Sichuan Basin, the westernmost being in the Gaoligong Mountains of northwestern (NW) Yunnan, an extension of the southeastern Himalayas (Fig. 4A). Though Davidia fossils have been found dating from Pliocene to Early Pleistocene sediments in central Japan19–21, none have survived there (Fig. 4A). We projected possible distribution of climatically suitable areas (potential habitats) for D. involucrata as expressed by occurrence probability (Fig. 4B). Our model showed excellent performance (AUC = 0.982 ± 0.006). Temperature seasonality (bio4), followed by the precipitation of the warmest quarter (bio18) and the annual mean temperature (bio1) (data not shown) were the most important variables determining the potential distribution of the species. Present-day records are located in potential habitats under current climate (Fig. 4B). The area of potential habitats, after accounting for uncertainty, is 534,953 km2 using the threshold of 0.0644, but the habitats under high probability (0.5‒1) account for 34.5% of the total area of potential habitats (threshold‒1) (Table 1). Only a very few small areas in southern and central Japan appear as suitable but with low probability (