Central Asian wild tulip conservation requires a regional approach, especially in the face of climate change

Tulips (Tulipa spp.) are one of the most widely appreciated plants worldwide, nevertheless species taxonomy and biogeography are often poorly understood. Most wild tulips inhabit the mountains of Central Asia, a recognised biodiversity hotspot, and a centre of tulip diversity. Despite the presence of several country-level endemic Tulipa species, most taxa span the borders of several nations. With no globally Red Listed tulip taxa from this region national level conservation assessments are an important resource. Nonetheless, threats posed to tulips are still inadequately understood, especially climate change, and given the trans-national nature of most species, distributional information is restricted and often misleading. Here we collate 330 species records from the Global Biodiversity Information Facility with 85 newly collected records, to undertake species distribution modelling (MaxEnt) for ten native Central Asian species. This work showed that regional level models provide a much more comprehensive understanding of species’ extinction risks, proportions of habitat in different countries, and limitations in protected area coverage. Furthermore, our climate modelling, the first of its kind for tulips, suggests that climate change will have a significant negative impact on the range size of all species; including those that are currently widespread. We therefore add climate change to the list of threats affecting tulip populations in Central Asia, which already includes livestock overgrazing, urbanisation, wild collection, and mining. Overall, our work shows that although national information is important, a regional approach is crucial not just for tulip conservation efforts, but likely for Central Asian plant conservation in general.


Introduction
Conservation, including conservation assessments such as the global IUCN Red List (IUCN 2020), are underpinned by an array of data. Primarily, they rely on knowledge of a species' distribution and threats to its survival (IUCN 2012). For many species, global IUCN Red Listing has not been undertaken, and this is often due to a lack of information to support these reports (Rodrigues et al. 2006). Local experts can provide a valuable insight into threats to species (Keppel et al. 2015), but commonly only country-level efforts have been undertaken to assess the status of flora, and many National Red Lists are therefore more comprehensive than global ones (Tojibaev and Beshko 2015). This is exacerbated by plants being globally and nationally less commonly assessed than animals (Nic Lughadha et al. 2020), in part due to plant-blindness (Balding and Williams 2016), whilst crucially, many are known to be declining towards extinction and require urgent conservation action (Nic Lughadha et al. 2020). Therefore, there is a need for more information to support plant conservation assessments, especially across the full geographic range of a species. Species distribution modelling presents a method through which to improve understanding of the suitable habitat area of a species (Phillips et al. 2017), developing knowledge of its extinction risk, and promoting population monitoring, management, and related policy, especially across multiple countries (Pearson et al. 2014;Wilson et al. 2019).
Tulips are a group of perennial geophytic monocots (Everett 2013). They are widely recognised for their horticultural varieties which support a billion Euro industry (Christenhusz et al. 2013), but also for their wild species which form a genus that is estimated to comprise between 76 and 90 species (Everett 2013). This array of wild taxa has historically underpinned the large horticultural trade, something that has significantly complicated the taxonomy of this clade (Christenhusz et al. 2013). Wild species remain a significant genetic resource for horticultural breeders, have considerable cultural value, and play an important role in ecosystems, especially for pollinators and insects (Kashin et al. 2016;Su et al. 2020). Wild tulips grow in the temperate regions of the northern hemisphere. Their distribution covers much of Eurasia, extending from western China across to western Europe, whilst a single species occurs on the Mediterranean coast of Africa. Nonetheless most species distributions centre around Central Asia (Botschantzeva 1982), and specifically the Mountains of Central Asia biodiversity hotspot (Critical Ecosystem Partnership Fund 2016).
Currently, there are only five tulip species published on the global IUCN Red List although none of these are native to the Central Asian centre of diversity (IUCN 2020). For Central Asian tulips an assortment of national level assessments have been published (SAEPF et al. 2006;Khassanov and Prastov 2009;Baitulin 2014;Tojibaev and Beshko 2015;Nowak et al. 2020), but these focus on country-wide distributions, not accounting for the fact that many species' distributions cross national borders. Poor representation is a common problem for Central Asian species, with data often lacking (Yesson et al. 2007;Paton et al. 2020), and international evaluation rare due to geopolitical tensions between neighbouring countries (Nowak et al. 2020). Regardless, national conservation assessments undertaken in this region indicate that a range of tulip species are threatened (SAEPF et al. 2006;Khassanov and Prastov 2009;Baitulin 2014;Tojibaev and Beshko 2015;Nowak et al. 2020). Threats previously recorded for these species include wild collection and trade, livestock overgrazing, and climate change (Nowak et al. 2020), although there remains limited literature and research focused on understanding their impacts. Nonetheless research suggests that the life history, ecology, 1 3 and cultural value of tulips makes them exceptionally vulnerable to disturbance and consequently to extinction (Tojibaev and Beshko 2015;Nowak et al. 2020).
Most tulips grow in the lower and middle elevations of mountain belts (Botschantzeva 1982). These alpine regions are thought to be extremely sensitive to climate change (Rangwala and Miller 2012) and therefore tulips may be especially vulnerable to this increasing threat (Nowak et al. 2020). Nonetheless the lower semi-desert and steppe areas of Central Asia, where fewer tulip species grow, are also reportedly fragile to a changing climate with many areas predicted to becoming increasingly arid in future years (Lioubimtseva and Henebry 2009;Chen et al. 2019). The impacts of climate change in these areas could be exacerbated by tulips' geophytic growth habit, relying on bulb-driven rapid spring growth to survive summer drought conditions common in the temperate latitudes where these species grow (Botschantzeva 1982). However, tulip bulbs also rely on a cold winter period as a trigger for initial growth, a process known as vernalisation. Furthermore, tulips require dry and freely draining soil conditions, with dampness often leading to rot and disease (Wilford 2013). Thus, tulip distributions are tightly linked to seasonal triggers, and both temperature and rainfall patterns, meaning changes in these may lead to declines in population numbers and even local extinctions. Moreover, climate change may cause shifts in suitable habitat and therefore there may be an increasing need for species to rapidly adapt or relocate to survive.
Initial flowering of tulips will not occur until there is a large energy store in the bulb, which can take several growing seasons. In addition, if damage occurs to the bulbs or leaves during the short growing season this can greatly weaken the plant, limit growth and reproduction, and sometimes even lead to death (Wilford 2013). Tulips' relatively long reproductive cycle and vulnerability to damage means that colonisation and repopulation of areas is slow, exacerbated by their limited ability to disperse pollen and seeds (Kashin et al. 2016). Given this and the widespread nature of livestock grazing across the grasslands, pasturelands, shrublands, steppes, and semi-deserts in which this plant commonly grows and the opportunistic collection of flowers by communities, tulips may be at increased risk from disturbance and may be unable to migrate to compensate for climate change (SAEPF et al. 2006;Tojibaev and Beshko 2015;Nowak et al. 2020). Additionally, many species are thought to have small, restricted distributions, especially many described endemics that are often only known from a specific hillside or gorge (Millaku and Elezaj 2015;de Groot and Tojibaev 2017), a trait widely associated with a heightened risk of extinction (Pearson et al. 2014).
Overall, Central Asian tulips are likely broadly threatened, but the extent and shape of this threat is significantly underreported, especially with regards to the impact of climate change. Strikingly, most species span the borders of the mountains of Central Asia, an ecosystem thought to be particularly vulnerable to the complex impacts of a changing climate (Xenarios et al. 2019;Nowak et al. 2020), and it is of increasing importance to expand upon national efforts to provide a regional perspective on the distribution and threats to wild tulips. This is especially crucial given that new species are frequently described in the region which often have extremely small distributions and may therefore have an immediate risk of extinction (de Groot and Tojibaev 2017;de Groot and Zonneveld 2020). Strikingly the genus Tulipa is only one of many geophytic clades that has a large diversity of species in the region including Amaryllidaceae, the broader Liliaceae of which Tulipa is only a small section, Iridaceae, and Asphodelaceae (Tojibaev et al. 2018). There are also a number of plant communities unique to the region including the walnut-fruit forest (Wilson et al. 2019) and the mountain grasslands (Borchardt et al. 2011). It is therefore likely that the issues surrounding tulips and 1 3 the corresponding threats are not exclusive to this plant group and are broadly indicative of the state of flora in this remote corner of the world.
Here we use species distribution modelling with MaxEnt to examine the range of ten Central Asian tulip species and the predicted changes in habitat suitability linked to climate change. Crucially, this work provides the first regional approach to tulip conservation in Central Asia, showcasing how such a practice can provide a more robust evidence base for conservation decision making in this region.

Study site
Fieldwork for this study was undertaken in the Republic of Kyrgyzstan in the spring of 2019 and 2020. We performed several field surveys: covering south-western Kyrgyzstan specifically the Batken and Osh regions; western Kyrgyzstan specifically the Jalal-Abad and Talas regions; and northern Kyrgyzstan, explicitly the Chuy and Issyk-Kul regions (see Online Resource 1). This fieldwork recorded 85 new location points for the ten species of focus (Table 1). The broader area of Central Asia, specifically the countries of Afghanistan, China, Kazakhstan, Kyrgyzstan, Mongolia, Tajikistan, Turkmenistan, and Uzbekistan, were included in species distribution modelling (Fig. 1).

Species distribution modelling
MaxEnt v3.4.0, based on the maximum-entropy approach, is a widely used modelling technique, especially in conservation (Trisurat et al. 2013;Liang et al. 2017;Wilson et al. 2019). This software is open source, can model past, present and future species distributions given suitable environmental layers, and relies only on presence data (Phillips et al. 2006). Specifically, the model uses location data and habitat vectors to predict the probability of presence of a species across a selected area. MaxEnt has been used frequently in conservation to model species distributions primarily because it is highly accurate with small sample sizes (Elith et al. 2011;Qin et al. 2017), characteristic of Threatened species. This software also models distributions under future climate change scenarios, allowing this threat to be assessed (Qin et al. 2017; Hof and Allen 2019). MaxEnt's inferences are correlative, with the software using a regression framework to produce predictions of occurrence. This method of climate change modelling often does not provide the same resolution as both mechanistic and trait-based approaches. Yet, these other approaches rely on detailed information of taxon-specific parameters, population sizes, interspecific relationships, and well-defined species distributions (Pacifici et al. 2015). In most cases this type of data is just not available for tulips, or not available in sufficient detail. This is especially the case in Central Asia, which is a relatively data deficient region (Pearson et al. 2014). So, although correlative approaches do not provide the same resolution, they require fewer initial data and therefore can be exceptionally useful for not only modelling present distributions but also future habitat in data poor areas.
In this study we focused on ten Tulipa species representing a range of habitats, distributions, and threat levels. We primarily selected species that had over 20 datapoints available for modelling to ensure there was enough GPS points present to establish a significant relationship with environmental variables. For comparison, we also selected to model T.  jacquesii (Zonneveld 2015) as a representative of the recently described endemic species in the region that is relatively data deficient (Table 1). Location data for each species was downloaded from the Global Biodiversity Information Facility (GBIF; GBIF. org 2020) database through the GeoCAT tool (Bachman et al. 2011) and combined with data gathered in our field surveys. We selected 23 environmental variable layers to be used as inputs for the MaxEnt programme. These consisted of 19 bioclimatic layers from the WorldClim2 database with a resolution of 30 s or ~ 1km 2 at the equator (Fick and Hijmans 1 3 2017), high resolution altitude data from the Shutter Radar Topography Mission (SRTM), slope and aspect layers generated from the altitude data using the QGIS 2.14 terrain analysis tool (QGIS Development Team 2009), and the land cover data GlobCover2009 (Arino et al. 2012). These layers have previously been effective in determining the distribution of Threatened species where data are limited, including in Central Asia (Kumar and Stohlgren 2009;Wilson et al. 2019). All layers were checked for multicollinearity using the ENM-Tools R package (Warren et al. 2019). A single variable from a group of highly correlated variables was chosen for modelling using a threshold (r > 0.85) commonly applied in Max-Ent work (Syfert et al. 2013;Wilson et al. 2019). We selected the variable from within this group which showed the average highest correlation to all other variables in the group; within this we avoided selecting altitude because this data would be uninformative in any climate modelling as linking a species presence to altitude would not allow for migration of a species.
Multiple location points of a species in the same grid cell were deemed duplicates and were treated as a single point. This reduced the available data for training and testing the model (Table 1). Given the limited data of some species, different modelling features were applied. For species with under ten points we used only linear features, for species with < 25 points we used linear and quadratic features, whilst for those with ≥ 25 datapoints linear, quadratic, and hinge features were used in modelling (Elith et al. 2011). K-fold cross-validation (K = 5) was used to generate an average present-day model for each species, as this has been empirically shown to neither be affected by excessively high bias nor very high variance (James et al. 2013) and reported as better than simple thresholding (Merow et al. 2013). The model's accuracy was determined by the test data area under the receiver operating characteristic curve (AUC) value with a value of > 0.9 representing a very good model, a value between 0.7 and 0.9 showing a good model, and anything with < 0.7 deemed uninformative (Swets 1988;Baldwin 2009). We assessed a range of regularization multipliers for several species deeming the default (1.0) to be the most effective modelling parameter based on both TestAUC (see Online Resource 2) and current understanding of taxa distributions (Everett 2013), and so this was used in all models. We also generated a dataset of locational data for Liliaceae across Central Asian using GBIF, which we used to perform background manipulation to account for sampling bias (Syfert et al. 2013;Kramer-Schadt et al. 2013). This generally helps to ameliorate overfitting, however, in our case it distorted distributions somewhat, both broadly lowering the TestAUC value (see Online Resource 3) and making predictions of suitable habitat in areas where the species is known not to occur (Everett 2013). We display only present-day modelling efforts and final climate models produced without a bias file. All outputs were presented using a five tier classification system and using the protected area data from Protected Planet (UNEP-WCMC and IUCN 2020).
Climate models were produced for nine of the ten species; T. jacquesii was excluded from this analysis due to its extremely limited location data. Two climate models were selected for comparison, Community Climate System Model version 4 (CCSM4 GCM) and Model for Interdisciplinary Research on Climate (Miroc ESM), to allow us to assess the reliability of our results across models. We selected these as they have both been used in previous studies to assess changes in habitat suitability linked to climate change (Rej and Joyner 2018; Hof and Allen 2019). We also selected two climate change scenarios to model, RCP 2.6 which represents a best-case scenario (BC) where emissions peak in 2020, and RCP 8.5 which is a worst-case scenario and similar to business as usual (BU). We selected to use data from both the years 2050 and 2070 so as to investigate the change in species distribution across multiple future periods. Given the number of models that needed to be generated we trained climate models using 75% of location data and tested them using the remaining 25% avoiding time consuming five-fold cross-validation (Wilson et al. 2019). All climate models were assessed using the same TestAUC classifications as the present-day models and presented in the same format.
Areas of habitat suitability were calculated for native and non-native sections of the models produced. The native area of a species was estimated based on the model output, location data, and literature (Bachman et al. 2011;Everett 2013;POWO 2019). To capture the native area, we created polygon layers on QGIS covering all location data and connected areas deemed suitable in habitat by the MaxEnt outputs and within or closely adjoining the previous estimated range of the species. We maintained a lenient approach to ensure we did not exclude important parts of the natural distribution from calculations and so likely captured some areas outside of the true distribution of the species. Non-native areas were therefore those outside the range of this polygon. Using these polygons, we extracted the number of cells in each of the five suitable habitat categories for each species. In addition, we used the estimated native areas to extract altitudinal information and calculate protected area coverage. To do this we selected cells from our models with a 0.5 or greater habitat suitability within the predefined natural distribution of the species, and specifically for protected area calculations those that overlapped with recorded protected areas from Protected Planet (UNEP-WCMC and IUCN 2020). We selected this threshold value as it has been previously used in research to represent a presence-absence separation (Carrasco et al. 2020) and would capture areas of habitat deemed to be highly suitable or very highly suitable, which we think represent the most important areas of the species range. To statistically compare altitudinal values from different years, we conducted an analysis of variance (ANOVA), which if significant (p < 0.05) was followed by a Tukey's HSD test. We used the area function of the raster 3.1.5 package (Hijmans 2014) to calculate mean cell size in km 2 across the modelled region. We then used this value, 0.648 km 2 , both for general area calculations as well as alongside the calculated number of cells meeting the threshold criteria to estimate habitat within protected areas. For both altitudinal and area computation raster layers were manipulated on QGIS v2.18.25 (QGIS Development Team 2009) and all calculations were carried out on R v3.4.0 (R Core Team 2020).

Results
After assessing for multicollinearity, we selected to use 14 layers for present-day modelling, whilst landcover was excluded from any climate modelling as future land cover predictions were not available leaving 13 layers for this aspect ( Table 2). All present-day models had TestAUC values larger than 0.9 and were therefore classified as very good (Table 3; Fig. 2). Yet, we report that the distribution of T. jacquesii be used cautiously as only four location points were available to train the model. The most important environmental predictors for every species varied greatly, however, each were linked to precipitation or temperature patterns in the winter or summer months or the seasonality of climate in general (Table 4). Models broadly highlighted areas of mountain ranges in the Tien Shan as suitable habitat, although a few species were more closely linked to semi-desert areas. Some species had relatively extensive distributions across the region, occurring in several mountain ranges, yet there was also several taxa which were much more spatially restricted and occurred in single mountain ranges, valleys, or steppe areas. Our models also show that tulip species are often not spatially separated, with considerable overlap between many 1 3 species' ranges. All species modelled, barring T. jacquesii as its model was deemed unreliable, were predicted to occur in at least one protected area. Nonetheless five species had less than five percent of their range in protected areas and both T. korolokowii and T. ferganica had two percent or less captured in the protected area network (Fig. 3).
For climate modelling, all models had TestAUC values greater than 0.9 and were deemed 'very good' providing informative results (Table 3). We present only the results from the CCSM4 GCM modelling (Fig. 4) as the results of Miroc ESM were consistently similar to CCSM4 GCM models not only across climate scenarios, but also years, and taxa (see Online Resource 4). There were two species reported that appeared to retain or expand in suitable habitat area, T. biflorformis and T. kaufmanniana, however, much of the future suitable habitat area for these species occurred in the Pamir mountains of Tajikistan, an area hundreds of kilometres from their native range. Given this information all models consistently showed a considerable reduction in the size of suitable habitat for tulip species in their natural distribution (Table 5; Fig. 5; see Online Resource 5). The suitable habitat area in the native species range of T. bifloriformis, T. dasystemon, T. greigii, and T. kaufmanniana declines in a stepwise manner from present-day to 2050 to 2070, with on average 78% of high and very high suitable habitat areas lost by 2050 and 83% lost by 2070. Whilst suitable habitat for the semi-desert dwelling species, T. korolkowii and T. kolpakowskiana and the alpine species of T. turkestanica, T. heterophylla, and T. ferganica is predicted to disappear completely. Overall, BC and BU scenarios were broadly similar however generally BU scenarios showed marginally less suitable habitat than the BC scenarios (Fig. 4).
Under future climate scenarios protected area coverage in the species native range decreased for all species except T. bifloriformis. For seven out of the nine modelled species coverage dropped to below one percent by 2050 and for six no suitable habitat was protected ( Fig. 3; see Online Resource 6). Our analyses also revealed that there was a significant difference between the mean altitude (metres) between the present day, 2050 and 2070 for T. bifloriformis [F(2, 296,712) = 224,119, p < 2.2 × 10 -16 ], T.   (Fig. 6; see Online Resource 7). More specifically, the altitudinal range significantly narrowed in future years for all species, whilst the suitable habitat areas for T. bifloriformis, T. dasystemon, and T. kaufmanniana were predicted to shift to higher altitudes. Surprisingly and unlike these other species the suitable habitat for T. greigii in 2050 was predicted to occur on average at lower altitudes than the present day, specifically at the base of previously more broadly suitable mountains. Nonetheless, the area of suitable habitat in 2070 was then predicted to shift back to near the present day mean altitude. All comparisons between years for these four species were deemed significant through a post-hoc TukeyHSD test (Table 6).

Discussion
Central Asia is the primary diversity hotspot for wild tulips and many species in this region have an elevated risk of extinction. Several new species have recently been described with highly restricted ranges. However, distributional understanding of tulips across this area is often poor and threats are inadequately reported. This study is the first to take a regional level approach to model current distributions of Tulipa, including one newly described endemic species, and assess how habitat suitability may change under different climate scenarios. Our models highlight a range of important results both for present day modelling as well as under future climate scenarios that allows us to draw a number of conclusions, primarily about tulips, with implications for the wider plant community and conservation. We recognise the limitation of modelling approaches to current and future species distributions, nonetheless our models provide an important resource, especially to aid future Red Listing efforts for the Tulipa genus, and to guide appropriate conservation interventions.
First, our models showcase the tight link between most Central Asian tulip's distributions and the mountain ranges of this region in both current and future climate scenarios. Our work underscores the importance of mountains in the niche occupancy of tulips in line with previous studies (Botschantzeva 1982). Far fewer species inhabit the lowland steppes and semi-desert areas of the region (Everett 2013), and so although these may still be of importance for conservation efforts, targeting mountainous areas is more urgent given the limited resources available to conservation practitioners (Bottrill et al. 2008). Moreover, our models highlight that many species distributions overlap across these alpine areas. This includes a number that are superficially similar showcasing the taxonomic difficulties of this genus, which are so often reported (Zonneveld 2009;Christenhusz et al. 2013). Given this, we urge researchers to be cautious when using tulip location data, especially where not supported by herbarium specimens, as well as our models and to critically assess these based on the current known natural distributions. This identification problem is exacerbated by the inconsistent use of the taxonomy of tulips (Christenhusz et al. 2013). We currently recommend using the species concepts of Christenhusz et al. (2013) so as to ensure consistent use of names across the scientific community, until further taxonomic work can update species concepts, which we note is urgently needed.
Many of the transnational species modelled in this project are reported Threatened across parts of their range. Yet, our results highlight that, frequently, the countries in which they are reportedly most Threatened often harbour only a small proportion of the overall distribution, and potential distributions under climate change. This is especially apparent for T. bifloriformis, T. korolkowii and T. kaufmanniana which are recorded in Tajikistan as Endangered, Endangered, and Critically Endangered respectively, yet they are recorded in an exceptionally small area of northern Sughd region, which may represent the extremity of their range. Nonetheless, this is also the case for other species such as T. greigii and T. dasystemon. This trend highlights that relying on national assessments for an understanding of the extinction risk of the whole species may be misleading and that global   assessments provide a much more informative and reliable approach. Nevertheless, remote regions of a species' distribution must still be considered in conservation planning as they potentially represent important sites of local adaptation and therefore genetic novelty (Flanagan et al. 2018). We also note through our work that several species are not evaluated in countries where they are reported as native, for example T. dasystemon in Kyrgyzstan. This is often because the species is widespread and national documents only focus on Threatened species (SAEPF et al. 2006). However, many of these taxa are considered Threatened elsewhere in their range e.g., T. dasystemon is Vulnerable in Uzbekistan. National level assessments may therefore present species as Threatened and in need of urgent action when across their broader range they could be considered relatively secure. Our work reinforces that although national level information remains an important resource, it needs to be critically assessed and considered in a broader context for use in directing conservation actions for non-endemic species. Furthermore, our climate models show that suitable habitat in future scenarios will remain trans-national and so international efforts will be crucial for tackling the impacts of climate change in Central Asia.
Due to its recent description, the endemic T. jacquesii, unlike the other species modelled here, had very limited location data available. We decided to model this species even given its limited GPS data to present an understanding of the challenges associated with a newly described taxon, as in the past decade a number of new tulip species have been described representing a considerable degree of newly discovered diversity (Tojibaev et al. 2014;de Groot and Tojibaev 2017;de Groot and Zonneveld 2020). Currently, like T. jacquesii, these taxa generally lack location information and conservation assessments. Our modelling of T. jacquesii supports previous research suggesting that predicting distributions with the extremely low number of points is highly constrained (Pearson et al. 2007);  in our case the predicted range was much larger than expected. Our work therefore importantly highlights the need for efforts to explore distributions for recently described tulip species to enable more accurate modelling and assessment of true distribution size. This forms part of a broader need for more information about these species to facilitate reliable assessment of their conservation status. Our T. jacquesii model provides a resource to aid in the search for new populations of this endemic species (Fois et al. 2018), albeit lacking significant resolution. Protected areas remain essential to conservation efforts globally (Naidoo et al. 2019), as reflected in policy in the Aichi Biodiversity Targets (Venter et al. 2014), and present a useful tool for safeguarding tulip populations. Using our models, we explored the overlap between predicted distributions and coverage of protected areas. Overall, this work emphasises the poor coverage of the protected area network of Central Asia in capturing tulip diversity. In general, most species have only been reported in one protected area (GBIF. org 2020) and our models support the view that only small parts of most species' distributions are captured in this network. Nonetheless our models do highlight that most species likely occur in several protected areas, but not always with confirmatory location data. Further efforts are needed to document the presence of species in many protected areas across the region. For example, our models of the species T. korolkowii and T. ferganica underline the restricted representation of these taxa in protected areas. Current knowledge also suggests that T. jacquesii does not occur in any protected areas, yet our model lacks the resolution to confidently assess this. Given the importance of protected areas for plant conservation (Chape et al. 2008;Souza and Prevedello 2020), the limited coverage provided for Threatened tulips needs to be addressed. Here, our models together with previous work (Botschantzeva 1982;Everett 2013) show that large areas of suitable habitat for these species are situated away from settlements in remote mountainous areas where increased protected area coverage may be feasible (Venter et al. 2018). These remote areas form part of the broader Mountains of Central Asia biodiversity hotspot (Critical Ecosystem Partnership Fund 2016) and so protection of these habitats may improve the survival chances of an array of Threatened species (Nowak et al. 2020). Even so, it is important to recognise that protection of lowland areas, especially semi-desert areas, will clearly also be essential for conserving species such as T. korolkowii, that are currently overly exposed to extinction due to the significant underrepresentation of their habitat in the Central Asian protected area network. Importantly, protected areas will not offer a silver bullet, with populations known to have declined in some strictly protected reserves (Krasovskaya and Levichev 1986), therefore a combination of conservation actions will need to be put in place alongside a strengthened protected area network. Regular monitoring of populations and stronger enforcement of environmental laws will also be critical components of successful implementation.
We undertook climate change modelling of this region to offer the first ever perspective on how this threat may impact future tulip habitat suitability. Across all models, seasonality, and precipitation or temperature patterns in the winter or summer months were deemed important predictors of distribution and this emphasises the importance of seasonal triggers in the life cycle of tulips (Botschantzeva 1982). Broadly, our models show that areas of habitat suitability will decline in all species including even the widely distributed and relatively common species, such as T. dasystemon, which exhibit a significant loss of suitable habitat in their native range. The severity of these declines is captured most clearly for T. bifloriformis which showed the lowest recorded loss of native habitat under future climate scenarios, yet even in this species only 40% of the present-day area of 'high' and 'very high' suitable habitat areas was predicted to remain in 2050. Clearly, climate change poses a significant threat to all tulip diversity in this region, mirroring the situation of many plants worldwide (Parmesan and Hanley 2015). We note that BA models were broadly worse than BC and so climate change mitigation may play a role in tulip conservation, but our models emphasise the severe plight of tulips even under best case climate scenarios. Although there is uncertainty surrounding our models, they reveal that distinct tulip habitats vary in vulnerability to climatic shifts. For example, our work shows that all semi-desert and steppe dwelling species are predicted to see a complete loss of suitable habitats by 2050, whereas only some alpine species show this. This is likely due to changes in rainfall patterns across these areas with aridity predicted to increase (Lioubimtseva and Henebry 2009). Given this information there is an urgent need to better protect populations of these semi-arid species now to allow genetic diversity to develop that may enable better resilience to climate change impacts in the future (Jump and Peñuelas 2005). Suitable habitat for some alpine species undergoes observable shifts to higher altitudes and declines in a stepwise manner as time progresses. This shift in altitude has been previously observed in different plant groups (Lenoir et al. 2008), but we provide supporting evidence that some tulip species may also show similar migratory trends. Our models suggest that this will also increase fragmentation of alpine refuges, leading to reduced gene flow between populations and an increased risk of extinction (Halloy and Mark 2003) escalating the need for more targeted conservation actions.
Protected areas and their expansion would likely play a significant role in the conservation of some species under future climate scenarios. Mountainous areas, including the mountains of Central Asia, are predicted to be extremely sensitive to the impacts of climate change (Rangwala and Miller 2012), yet across Central Asia there are a number of protected areas that encompass high altitudinal habitat, which notably are already connected to landscapes where tulips grow. Broadly our modelling shows that protected area coverage of species will decrease under future climate scenarios yet, they also suggest that several protected areas already encompass suitable habitat into which tulips may eventually migrate and so could be of increasing importance to populations as climate patterns begin to change. As an exception to this rule, T. bifloriformis appears to have more suitable habitat in protected areas in future climate scenarios. Even so, given the poor dispersal range of tulips (Kashin et al. 2016), a trait deemed important for survival in alpine areas (Rumpf et al. 2019), migration may be slow and could prevent species from reaching suitable habitat areas before dying out. There is already significant evidence that extinction debts and colonisation credits will be widespread in future climate change scenarios in mountainous regions (Rumpf et al. 2019). So, although several alpine tulips are predicted to have suitable habitat at higher altitudes, including within protected area, their survival may still rely on human intervention. Interestingly, suitable areas way outside of several species natural ranges were highlighted in our modelling. This was especially apparent for T. kaufmanniana and T. bifloriformis, whose native range currently encompasses the mountains around the Fergana valley, but where large parts of Tajikistan's more southerly Pamir mountains became suitable in 2050 and 2070. We therefore suggest that future species translocation initiatives (Berger-Tal et al. 2020) may be necessary, although considerable further work is needed to determine the effectiveness and ecological safety of such an action.
We note here that our models do not account for a range of factors. Genetic variation and species adaptability to climate change has not been incorporated, but can be important for persistence in areas deemed unsuitable (Graae et al. 2018;Razgour et al. 2019). Moreover, previous research has shown tulips can actively populate and survive in highly disturbed landscapes including agricultural land (Krasovskaya and Levichev 1986;Pratov et al. 2006) and therefore may survive better in a changing landscape than 1 3 our models predict. Although changes in climate may decouple seasonal triggers such as flowering time (Wadgymar et al. 2018) which could be exceptionally damaging for tulips and similar plants that rely heavily on these for the timing of their short growing season. Furthermore, alpine habitats encompass a range of microclimatic niches which broad scale modelling overlook as potential refuge areas (Scherrer and Körner 2011). Some areas deemed unsuitable may therefore in fact present adequate microclimatic conditions for the survival of local populations. The structural composition of communities, which is especially important to tulips due to their requirement for direct sunlight for growth, may mean that areas within predicted suitable habitat cannot in fact support populations (Vittoz et al. 2009). We therefore acknowledge these limitations and accept that some taxa may be more resilient than suggested by our models. Nonetheless, we suggest climate change will be an important threat to tulip populations, and highlight that there are other factors we have not examined, such as the shifting of invasive species into mountainous areas, that could exacerbate impacts further (Petitpierre et al. 2016).
Climate change is not the only threat posed to wild tulip species. Poorly managed livestock can cause significant damage to ecosystems (Wilson et al. 2019) and livestock overgrazing continues to degrade habitat across much of Central Asia (Tojibaev and Beshko 2015;Nowak et al. 2020). Given that livestock populations are thought to be on the increase across Central Asia, overgrazing appears to pose a growing threat to tulips. Furthermore, although many settlements are in rural areas (Djanibekov et al. 2016) urbanisation also poses a threat to tulips. Many of the major cities in Central Asia are situated close to mountainous tulip habitat, including Bishkek, Almaty, Dushanbe, and Samarkand. Given the rapid development of these cities and the corresponding loss of habitat, urbanisation needs to be urgently considered as part of any tulip conservation activities in the region; a similar but more localised threat is presented by mining activities. Finally, the horticultural history of the genus and the demand for tulips worldwide has meant that wild collection and trade has been reported as a threat and is believed to have led to previous extinctions and populations declines (Maunder et al. 2001;Menteli et al. 2019). Central Asian tulips have been an important part of tulip horticulture throughout the existence of this trade (Christenhusz et al. 2013) and now many Threatened tulip species are protected by law (SAEPF et al. 2006). Yet, opportunistic collection continues, and this may exacerbate the impact of other threats including climate change.
Overall, here we have shown that climate change will pose a significant threat to wild tulips, whilst current distributions of most species are tightly linked to the mountains of the broader Central Asia region and are poorly captured in protected areas. This leaves many populations already declining, spanning borders that scientific research and conservation collaboration has not yet bridged, and increasingly exposed to an array of threats and their interactions. Whilst our work has focused on the genus Tulipa, and specifically Central Asian species, we recognise that many plant groups require similar focused attention and so, although we advocate for urgent efforts to protect wild tulips from growing threats, we also suggest that efforts are made to carefully assess and use available data, including national level assessments, to improve conservation of plants across broader Central Asia. Yet, most importantly in this paper we have shown that a regional approach is essential for an accurate understanding of a species' risk of extinction, especially with respect to the growing impacts of climate change. Given this, now is the time for the broader conservation community to work together to ensure a more aligned regional approach in Central Asia.