Projecting Suitability and Climate Vulnerability of Bhutanitis thaidina (Blanchard) (Lepidoptera: Papilionidae) with Conservation Implications

Bhutanitis thaidina is an endemic, rare, and protected swallowtail in China. Deforestation, habitat fragmentation, illegal commercialised capture, and exploitation of larval food plants are believed to be the four major causes of population decline of B. thaidina in the recent decade. However, little attention was paid to the impact of climate change. This study used ecological niche factor analysis and species distribution model to analyse the current suitable areas for B. thaidina with BioClim variables as well as its future suitable areas under four future climate scenarios (represented by four Representative Concentration Pathways: RCP2.6, RCP4.5, RCP6.0, and RCP8.5). Statistical analysis was carried out to compare the possible area and altitude changes to the distribution of B. thaidina under changing climate. Our analyses showed that the suitable areas for B. thaidina are fragmented under the current climate, with four suitable centres in northwestern Yunnan, northeastern Yunnan and northwestern Guizhou, the western margin of Sichuan Basin, and Qinling mountains. Apart from further habitat fragmentation under climate change, slight range expansion (average 6.0–8.9%) was detected under the RCP2.6 and RCP4.5 scenarios, while more range contraction (average 1.3–26.9%) was detected under the RCP6.0 and RCP8.5 scenarios, with the two southern suitable centres suffering most. Also, a tendency of contraction (2,500–3,500 m) and upslope shift (~600 m) in suitable altitude range were detected. The findings of this study supported the climate-vulnerable hypothesis of B. thaidina, especially under future climate like the RCP6.0 and RCP8.5 scenarios, in terms of contraction in suitable areas and altitude ranges. Conservation priority should be given to northwestern Yunnan, northeastern Yunnan, and northwestern Guizhou to alleviate the stress of massive habitat loss and extinction. Refugial areas should be established in all four suitable centres to maintain genetic diversity of B. thaidina in China.

Habitat losses associated with human activities are undoubtedly imminent threats to the survival of certain populations of these two Bhutanitis species in China. In recent years, a few conservation studies were carried out on B. thaidina in attempt to alleviate the situation from biological and ecological aspects 4,14,15 , while little could be done with B. mansfieldi, a bionomic and distribution data-poor species. Apart from human activities, climate change is another factor which is attributed to many cases of extinction of rare and endemic species globally 16,17 . However, unfortunately, little attention has been paid to such slow but prolonged effect of climate change on the future of these butterflies.
Species distribution models (SDM) contains a range of effective analytical tools for simulating and visualising suitable areas (potential distribution range) of organisms, and has been widely applied to species of conservation interests as well as policy making over the past decade [18][19][20][21][22][23][24][25][26][27][28][29] . Among these methods, ecological niche factor analysis (ENFA) and maximum entropy (MaxEnt) modelling are the most frequently applied SDMs which project the suitable area of a species using the presence-only data without depending on bionomical parameters of the focal species, or being biased by pseudoabsence data 30,31 .
In an attempt to fill the gap in conservation of Bhutanitis, the present study chose B. thaidina, a data-rich species as our model, analysed the current distribution and the future distribution shift under different climatic change scenarios 32 using SDMs of ENFA and MaxEnt. The results will provide us an overview of its suitable areas in China and facilitate our understanding of how the suitable areas would shift in the process of climate change. The findings of the present study are beneficial to conservation management in current time as well as to formulate countermeasures to alleviate population decline of this rare butterfly in the future.
Nineteen BioClim 38 variables were used to represent the current climate features (averaged over 1970-2000), the 19 BioClim variables and the altitude mask with 30 arc seconds resolution were obtained from the WorldClim database (www.worldclim.org). All data was further cropped by the political boundary of the People's Republic of China, and will be referred as 'environmental factors' hereafter.
The CMIP5 climate projections under the IPCC-AR5 (the 5 th Assessment Report of the Intergovernmental Panel on Climate Change) frame were used to represent the future climate 32 . Four representative concentration pathways (RCPs), RCP2.6, RCP4.5, RCP6.0, and RCP8.5 were selected to simulate possible climate changes 32 . Data with 30 arc seconds resolution was also obtained from the WorldClim database and cropped by the political boundary of the People's Republic of China.
Species distribution points and environmental factors were transformed into two formats, with the IDRISI format for ENFA analysis 31,39 and the ASCII format for MaxEnt analysis 40 . of the 20 environmental factors were tested using a UPGMA dendrogram in Biomapper 4.0, and factors with correlation coefficients above 0.95 were removed from the dataset. When removing autocorrelated factors, those representing short-period extremes (e.g., minimum temperature of the coldest month, maximum precipitation of the wettest month) were removed, while those representing longer periods (e.g., mean temperature of the coldest quarter, precipitation of the driest quarter) were kept, as such type of environmental factors often play an important role in species distribution.
In an attempt to analyse the distribution prevalence of B. thaidina, values of previously screened environmental factors in the distribution area of B. thaidina and the entirety of China were extracted using DIVA-GIS 5.7 (www.dive-gis.org) 41 . The distribution frequencies were calculated in Biomapper 4.0. The importance of the environmental factors was measured using the jackknife method in MaxEnt 3.4.1 with 1,000 iterations 42,43 .

Species distribution model (SDM).
Environmental factors with ENFA scores over 0.2 were selected and assigned to the MaxEnt 3.4.1 40 to project suitable areas (current and future climate variables were analysed separately) based on the presence data points, among which 25% were extracted for random testing. The logistic output method was selected to estimate the distribution (or presence) probability of B. thaidina considering certain assumptions of species' prevalence and sampling effort 44 . The resultant map was saved as ASCII format and then redrawn using Surfer 10.0 (Golden Software Inc., Golden, CO, USA). Model robustness was evaluated using the receiver operation curve (ROC) and the area under the ROC curve (AUC) 45,46 , where the AUC value [AUC ∈ (0, 1)] approaching 1.0 is usually considered acceptable, whereas it should be rejected when approaching the random turquoise line of 0.5 47 .
Statistical analyses. The number of grid cells (further transformed into area using 1 grid cell = 1 km 2 ) as well as their elevation property were extracted in ArcGIS 10 (ESRI, USA) from projection maps under the current climate and the four future climate scenarios in both the 2050 s and the 2070 s. Comparative bar charts for suitable areas and curve line charts for suitable altitude range were made to 2050 s vs. current, 2070 s vs. current, and 2070 s vs. 2050 s, mainly focusing on suitability ranks from 0.5 to over 0.8.
Frequency distribution of B. thaidina against the nine key influential environmental factors for the entirety of China showed evident preference for each factor. For temperature factors, B. thaidina occurs in areas where annual mean temperature (Bio1) ranges between 1-17 °C, mean temperature of the warmest quarter (Bio10) ranges between 10-23 °C, and mean temperature of the coldest quarter (Bio11) ranges between −9-10 °C (Fig. S1). For precipitation factors, B. thaidina occurs in areas where annual precipitation (Bio12) ranges between 630-1,400 mm and precipitation of the wettest quarter (Bio16) ranges between 250-750 mm (Fig. S1). For temperature variabilities, B. thaidina mainly occurs in areas where mean diurnal temperature range (Bio2) varies between 7-12 °C, temperature annual range (Bio7) varies between 22-35 °C, relatively higher isothermality (Bio3) and low temperature seasonality (Bio4) (Fig. S1). On the large scale, the current suitable areas for B. thaidina are still confined to west China, as mirrored by its actual distribution localities (Fig. 2). Four areas with higher suitability were identified. (1) Northwest Yunnan. This area occupies the Hengduan Mountains in Yunnan and southwest Sichuan, including the mountains separated by the upper Irrawaddy, Salween, Mekong, and Yangtze watersheds. The eastern edge of this area At smaller and local scales, the suitable areas for B. thaidina are highly fragmented, even within the four isolated patches mentioned above. Despite ridges of high mountains and deep valleys of large rivers cutting these patches into separate pieces, the suitable areas for B. thaidina are further isolated by complex terrains in a small range (Fig. 3). In the 2050 s, the overall distribution pattern of the suitable areas showed obvious but non-radical changes. The change under RCP2.6 scenario is very limited, making the distribution pattern very similar to that under the current climate, but the 0.7 grade suitable areas expanded in northwest Yunnan, occupying the 0.6 grade suitable areas under the current climate; while the 0.7 grade suitable areas retreated in southwest Sichuan bordering with northeast Yunnan (Fig. 4A). Under the RCP4.5 scenario, the distribution pattern remained the same with that under the RCP2.6 scenario, except for an elevation of suitability grade (0.8-0.9) in the western edge of the Sichuan basin (Fig. 4B). Under the RCP6.0 and RCP8.5 scenarios, distribution pattern changed more obviously with the low-medium suitability grades (0.5-0.6) retreating in the southern edge and lower-altitude areas of the distribution range, but a higher suitability grade (0.8-0.9) appearing in northwest Yunnan and the western edge of the Sichuan basin (Fig. 4C,D). Gain of high suitability grade (0.7) was also detected in south Qinling under the RCP8.5 scenario (Fig. 4D).

Current suitable areas. The
In the 2070 s, the overall distribution pattern of the suitable areas under the RCP2.6 and the RCP4.5 scenarios almost remained the same as that in the 2050 s (Fig. 5A,B). However, dramatic changes to the suitability  www.nature.com/scientificreports www.nature.com/scientificreports/ distribution under the RCP6.0 and the RCP8.5 scenarios were detected. Not only the low-medium suitability grades (0.5-0.6) largely retreated in the southern margin and the lower-altitude areas, the medium-high suitability grades (0.6-0.8) also dramatically retreated in the southern portion of the distribution range, especially in northwest Yunnan, northeast Yunnan, and northwest Guizhou (Fig. 5C,D). Similar to the 2050 s, gain of higher suitability grade (0.8-0.9) was detected in south Qinling under the RCP8.5 scenario (Fig. 5D).
Analysis of frequency distribution change of suitable altitudes showed that, under RCP2.6 and RCP4.5 scenarios, the frequency distribution of suitable altitudes did not shift obviously but were more contracted between 2,500-3,200 m in the 0.5-0.6 and 0.6-0.7 ranks, while they shifted to 2,600-3,500 m in the 0.7-0.8 rank (Fig. 7). The frequency of suitable altitudes in the 2050 s is almost equal to that in the current, while the frequency of suitable altitudes in the 2070 s was significantly higher than that in the 2050 s (Fig. 7). Under RCP6.0 scenario, the frequency distribution patterns were similar to those under the preceding scenarios, but the peak frequency of the suitable altitudes was lower, with that in the 2070 s even significantly lower than the 2050 s in the 0.6-0.7 and the 0.7-0.8 ranks (Fig. 7). It is noticeable that under the RCP8.5 scenario, the peak frequency of the suitable altitudes was further significantly lowered, especially in the 2070 s (Fig. 7).

Discussion
Distribution shift and climate vulnerability. Based on the EGV frequency distribution characters, B. thaidina mainly occurs in temperate climate zones with less precipitation, relatively higher diurnal temperature range, and lower temperature seasonality (Fig. S1). Such habitat is represented by montane broadleaf forest and subalpine evergreen needle leaf forest, which may extend from 2,000 m to below the treeline in west China, with its lower and upper limits varying from the south to the north 48,49 . With climate change, the habitat belt will be forced to 'move' upslope 50 . However, shift rates of the lower and upper limits are not expected to be the same, and such asymmetric shift rates will eventually result in a decline in suitable altitude belt 51 .
Although some research indicates that species will move or expand their ranges upslope and poleward with climate change 52-56 , our analysis implied that B. thaidina would more likely suffer from rapid habitat compression or be driven to extinction during this process [57][58][59] . On the one hand, the distribution shift rate of species can hardly keep up with the pace of habitat shift and compression 51 ; and on the other hand, ascent of the upper part of its suitable range will be limited by unfavourable climate or soil condition (e.g., increased precipitation, www.nature.com/scientificreports www.nature.com/scientificreports/ erosion linked to permafrost degeneration) 60 . More practically, the ability of B. thaidina to move from one habitat to another is largely limited by lack of connecting mountain ridges with suitable habitat, since its suitable area is already highly fragmented even under the current climate. Our speculation on this issue was similarly demonstrated in the tropical forests 61 .
The adult is the only stage at which all butterflies can achieve distant movement, while other stages from egg through pupa can only stay in the same locality. Such a life history further reduces their ability to escape from bad climate. Bhutanitis, a univoltine group with a maximum adult stage of only 2.5 months per year 14,62 , can hardly defend themselves against any climate induced incidents, including extreme precipitation, long-lasting drought, or even forest fire 63,64 . The extinction of B. lidderdalii on Doi Chiang Dao of northern Thailand in 1983, caused by a severe forest fire in the dry season as a result of the 1982-1983 El Niño, is a most recent case (A. M. Cotton, pers. comm.).
As a result, the distribution shift for B. thaidina in light of climate change would compress its suitable habitat and further reduce its refugial areas (Figs 3-7). Globally, the climate vulnerability of B. thaidina is much higher under the RCP6.0 and RCP8.5 scenarios compared to that under the RCP2.6 and RCP4.5 scenarios. Regionally, climate vulnerability of the two southern distribution centres is higher than the two northern ones (Figs 4 and 5) (discussed in detail below).
Biodiversity significance. The current suitable areas revealed a patchy and highly fragmented distribution pattern for B. thaidina (Fig. 3), while future projection showed a compressed and further fragmented distribution pattern (Figs 4 and 5). Nonetheless, as the present study only applied climatic factors in SDM simulation, the actual distribution pattern of B. thaidina could be more fragmented on an unsuitable matrix when availability of host resources and vegetation are taken into consideration. In population genetics, highly fragmented distribution would result in a reduction of gene flow and genetic diversity [65][66][67][68][69][70][71] .
B. thaidina is a morphologically variable species with four subspecies recognised to date: ssp. thaidina in west Sichuan, ssp. hoenei Bryk in northwest Yunnan, ssp. melli Bryk in Qinling and Taibai Shan (probably also in Shennongjia), and ssp. dongchuanensis Lee in northeast Yunnan and northwest Guizhou 3,9,11,33,35,72 (Fig. 2), mirroring our identification of four suitability centres (Figs 3-5). The biological and ecological issue underpinning the taxonomic complexity is that each subspecies possesses a distinct genetic profile. The genetic and www.nature.com/scientificreports www.nature.com/scientificreports/ morphological profiles altogether constitute the entire biodiversity integrity of B. thaidina in China, and any degeneration or loss of such profiles will directly lead to a loss of biodiversity of this endemic species.
The distribution pattern of suitable areas for B. thaidina is highly fragmented, thus making each subspecies a metapopulation comprised of multiple scattered and isolated smaller local populations. A recent population genetic analysis suggested very low genetic diversity among all populations of B. thaidina in China 73 , implying vulnerability to degeneration or extinction in the dynamic wild 74 . In the process of climate change, the future distribution pattern of B. thaidina will be further fragmented and isolated, which would inevitably bring more restriction to the gene flow between the four distribution centres as well as within each one. Our future projections showed significant suitability loss in the distribution range of ssp. dongchuanensis in northeast Yunnan and northwest Guizhou, followed by ssp. hoenei in northwest Yunnan (Figs 4 and 5), making these two subspecies more prone to extinction under the RCP6.0 and RCP8.5 scenarios, compared to the other two subspecies.
Conservation implications. The present study showed that the fragmented suitable areas for B. thaidina will undergo further fragmentation and reduction in the process of climate change (Figs 3-5). Hence, maintaining current existing suitable areas is vital to the conservation of this rare species.
To conserve B. thaidina with genetic integrity, conservation strategies must firstly take all four suitable centres into equal consideration, as each suitable centre represents a distinct subspecies of B. thaidina in China. Next, combining the degree of rarity and vulnerability, priority should be given to areas with ssp. dongchuanensis and ssp. hoenei, which will be the most threatened in the future (Figs 4 and 5); followed by ssp. melli, which is only found in a narrow area in Qinling (Fig. 3); and ssp. thaidina being of least concern.
Availability of larval food plants is also crucial in conservation of B. thaidina. This species shows a strong host association in nature, and subspecies use different Aristolochia species, e.g., in northwest Yunnan and most parts of west Sichuan, B. thaidina mainly uses A. moupinensis and A. delavayi 3,14 , while using A. mandshuriensis in the Qinling Mountain area 4 . These food plant species have being exploited for traditional Chinese herbal medicines until recent years 75 . Such long lasting exploitation has already depleted wild resources of Aristolochia in some places 2 . Apart from human exploitation, deforestation of virgin forests on the median altitude to subalpine mountains is another important threat to wild Aristolochia resources, as most Aristolochia species are shade plants. Deforestation will also destroy the suitable habitats of Aristolochia.
Establishing refugial areas for B. thaidina could be an optimal in situ protection method. When selecting sites for refugial areas, vegetation surveys must be performed in advance to ensure the best vegetation type being included, e.g., Quercus stands associated with multiple local Aristolochia species 76 . However, the optimal planning of refugial areas for B. thaidina must rely on future in-depth research involving bionomics, dispersal capability, food plant adaptability, habitat matrix composition and connectivity, as well as a thorough evaluation of population genetic diversity.
When in situ protection faces the challenge of high and progressive habitat fragmentation due to climate change found in the present study, coupled with the limited genetic diversity described earlier 73 , other measures must also be considered in the future to increase the genetic diversity and evolution flexibility of B. thaidina to respond to rapid environmental changes in a certain area. Possible measures include introducing ex situ populations containing new genetic profiles from other distribution areas, or even releasing laboratory genetic recombinants 80 . Published: xx xx xxxx