Potential changes in the distribution of Delphinium bolosii and related taxa of the series Fissa from the Iberian Peninsula under future climate change scenarios

A particular threat posed by climate change for biodiversity conservation, one which has scarcely been studied, is the overlap of the potential distribution areas in phylogenetically closely related species. In this study, Species Distribution Modelling (SDM) was used to investigate the potential changes in the distribution of Delphinium bolosii and D. fissum subsp. sordidum under future climatic scenarios. These two closely related and endangered endemic species from the Iberian Peninsula do not have complete reproductive barriers between them. The two models selected different predictors with a similar effect in the biological cycle. Both taxa need low winter temperatures to break seed dormancy and sufficient rainfall to complete the flowering and fruiting stages. The current potential distribution areas of both taxa do not currently overlap. However, the results showed that potential changes may take place in the species’ distribution range under future climate scenarios. The models predict a reduction of the potential distribution area of D. bolosii while, conversely, the potential distribution area of D. fissum subsp. sordidum increased. In both cases, the predicted contraction in range is very high, and loss of habitat suitability in some current localities is worrying. Notwithstanding, the models do not predict overlaps of potential areas under climate change scenarios. Our findings can be used to define areas and populations of high priority for conservation or to take action against the impacts of climate change on these endangered species.


Introduction
Recent studies suggest that climatic change may well be one of the most important challenges that humanity will face in forthcoming decades to avoid precisely what models forecast, in the worst-case scenarios -namely the sixth mass extinction in the history of the Earth (Barnosky et al. 2011). Climate change effects on biodiversity involve many components, such as: (1) decreased genetic diversity and/or populations; (2) alterations to interspecific relationships (e.g. plant-animal systems); and (3) changes in plant communities (Bellard et al. 2012).
Rare species are more susceptible to environmental changes (e.g. climate warming) than common species (Lomba et al. 2010). The International Union for Conservation of Nature (IUCN) recommends developing new quantitative methods to determine species ranges, potential suitable habitats and potential shifts in distribution ranges when faced with climate change, especially for rare and threatened species (Fordham et al. 2012). By means of species distribution models (SDMs), a technique that is being increasingly employed in the 21 st century and uses environmental variables together with species occurrences (Guisan and Zimmermann 2000), it is possible to obtain this information; e.g. ecological drivers of species' ecological niches and geographic distributions and predicting how future climatic change will affect biodiversity Thuiller et al. 2011). SDMs are crucial for anticipating the ecological impacts of climate change in more susceptible areas, and for devising adequate conservation planning and decision making (Riordan and Rundel 2014).
The interpretation and usefulness of SDMs have been discussed in several works (e.g. Sinclair et al. 2010;Austin and Van Niel 2011). It has been found that SDMs almost never consider interspecific interactions, species' genetic and plastic ability and temporal and physiological responses, despite some exceptions (such as Lavergne et al. 2010). Likewise, species are considered static elements, but their dynamism has been documented (Schröder and Seppelt 2006). However, Araújo and Luoto (2007) confirmed that models based only on environmental variables on a regional scale (e.g. the Iberian Peninsula) can adequately quantify the impact of climatic change on species distributions. For species with generalist requirements, such as heterogeneity of habitats and plant communities, broad spectrum of pollinators and non-specific seed dispersal adaptations, SDMs can offer good results to address several ecological and biodiversity conservation issues (Sinclair et al. 2010). This is the case with the endangered endemic Iberian species of the Delphinium L. series Fissa B. Pawl.
The origin of Delphinium series Fissa (Delphinieae, Ranunculaceae) lies in Central Asian mountains and extended westward throughout the Mediterranean Basin in the Late Miocene, namely during the Messinian Salinity Crisis (6-5.3 Ma) (Blanché 1991). At present, Delphinium series Fissa is represented in the western Mediterranean region by three taxa: D. fissum Waldst. and Kit. subsp. fissum is distributed from northern Italy and the French Maritime Alps to the Pyrenees, whereas two other taxa, Delphinium bolosii C. Blanché and Molero and D. fissum subsp. sordidum (Cuatrec.) Amich, E. Rico and J. Sánchez, are differentiated in the Iberian Peninsula , 2020. According to Bosch (1999), the reproductive barriers during this differentiation process are still incomplete, and the formation of fertile hybrids is possible through interbreeding. These two taxa have undergone an allopatric speciation process to currently give rise to non-overlapping distribution areas. One interesting application of SDMs is to study changes in the potential distribution areas of phylogenetically closely related species as the overlap of the potential areas induced by climate change might lead to hybridisation and competition processes, which could threaten species' persistence. However, these biological impacts of climate change have scarcely been studied (Krosby et al. 2015).
Recently, Rus et al. (2018) generated potential distribution models according to current climate conditions for D. fissum subsp. sordidum, and identified the environmental factors to explain its distribution, which are related to the taxon's biological and ecological characteristics. The selected predictors and its percentage contribution were: mean temperature of wettest quarter (26.2%), temperature annual range (25.7%), precipitation of warmest quarter (22.1%), slope (21%), precipitation of coldest quarter (2.7%), and hillshade (2.3%).
The present study aims to: (1) generate the current model for D. bolosii with the same methodology used for D. fissum subsp. sordidum, and compare the predictors of both models; (2) determine and evaluate the potential impacts of climate change on the distribution range of the series Fissa species in the Iberian Peninsula, as an example of closely related species without complete reproductive barriers; and (3) provide information for the medium and long-term conservation of these taxa under climate change scenarios.

Study area and test species
The study area corresponds to the Iberian Peninsula, which is located to the southwest of the European continent ( Fig. 1). This territory covers an area of 582918 km 2 and is characterised by its asymmetrical distribution of main mountainous systems and a strong climatic gradient. Most of the Iberian Peninsula has a Mediterranean climate, characterised by a summer with at least two very dry and hot months. Moreover, it has high levels of biodiversity due to its great bioclimatic and biogeographic heterogeneity determined by geological and historical factors (Rivas-Martínez et al. 2002). Consequently, the Iberian Peninsula is one of the floristically richest areas of Europe and the Mediterranean Basin (Aedo et al. 2013).
According to the phylogenetic results obtained by cpDNA (Ramírez-Rodríguez et al. 2019), field observations and reviewed herbarium material (Ramírez-Rodríguez et al. 2020), the populations located in the Pre-Pyrenean and Catalan Mediterranean System correspond to D. bolosii (Blanché and Molero 1983;López-Pujol et al. 2015), as well as the cited populations on the eastern slope of the Iberian System (Mateo and Pisco 1993, sub D. fissum subsp. sordidum;Pitarch 2002, sub D. mansanetianum). Delphinium bolosii is a perennial rhizomatous larkspur endemic of the north-eastern Iberian Peninsula (Fig. 1). This species lives in Quercus rotundifolia Lam. forests on schists on the edge of communities with Buxus sempervirens L., Jasminum fruticans L., and Acer monspessulanum L., and in open clearings in Quercus humilis Mill. forests ( Fig. 2A, B). It is catalogued as EN [B1ab(iii,iv,v)+2ab(iii,iv,v); C1] on the Red List of Spanish Vascular Flora .
Delphinium fissum subsp. sordidum is a rhizomatous hemicryptophyte endemic of the central-western and southern Iberian Peninsula ( Fig. 1) that lives in chestnut forests, oak forests, holm-oak forests and open areas between shrubs (Fig. 2C, D). The subspecies is included in the Red List of Spanish Vascular Flora  under the category EN [B2ab(v)c(iv); C2b]. For further details, see Rus et al. (2018) and Ramírez-Rodríguez and Amich (2019).

Occurrence data and environmental variables
Fourteen occurrences for D. bolosii were obtained, but those that occurred in the same 1 × 1 grid cell were removed to reduce spatial correlation (Benito et al. 2013). This left 12 valid occurrences at the defined grid resolution. Although the number of occurrences is low, they represent the whole distribution range of D. bolosii after decades of study of the species (Blanché and Molero 1983;Mateo and Pisco 1993;Bosch et al. 1998, . As potential predictors, 19 bioclimatic variables and one topographic variable (elevation) were considered and downloaded from the WorldClim database (Hijmans et al. 2005, available at http://www.worldclim.org). Three additional topographic variables (aspect, hillshade and slope) were obtained from the elevation variable by ArcMap 10.3.1 (ESRI, Redlands, California, USA). A list of all the used variables is available in Suppl. material 1: Table S1. These variables were the same as those used by Rus et al. (2018) to generate the potential model of D. fissum subsp. sordidum. No edaphic, phytocenotic, land use and landscape features variables were used due to the study taxa develop in a broad variety of soils and habitats.
For the future predictions, the CCSM4 model was used, which was downloaded from the Fifth Assessment Report of the International Panel on Climate Change (IPCC AR5 WG1 2013) for the years 2050 and 2070 with two IPCC Representative Concentration Pathways: "stabilisation" (RCP 4.5) and "high increase" (RCP 8.5). They were selected as they represent contrasting scenarios to predict potential climate change effects. Thus, the emissions in RCP 4.5 peak around 2040 before declining. In RCP 8.5, emissions continue to rise throughout the 21 st century (Meinshausen et al. 2011). CCSM4 is an efficient global climate projection that predicts the influence of future climatic changes on the distribution of plant species in the Mediterranean Basin (Abdelaal et al. 2019).

Spatial modelling
The model for D. bolosii was created using Maxent v. 3.3.3 (Phillips et al. 2006;Merow et al. 2013), following the same methodology used for D. fissum subsp. sordidum by Rus et al. (2018). For details about data preparation, Maxent configuration, model fitting, model evaluation and jack-knife procedure, see Rus et al. (op. cit.). This methodology is adequate for modelling species with a low number of occurrences (Pearson et al. 2007). To test the fit of model, the area under the curve (AUC) of the receiveroperating characteristic (ROC) was used (Fielding and Bell 1997). Although this metric has been criticised in some recent works (Jiménez-Valverde 2012), it is still the most widely used metrics, and is a good tool to evaluate models for the same species and within the same geographical scope.
To determine the potential distribution areas, binary outputs of presence/absence were generated by setting thresholds in the logistic models. Firstly, the 10 th percentile training presence logistic threshold (10P) was used to transform habitat suitability, as estimated by the models, into a binary prediction. This method is well recognised for distinguishing suitable regions from unsuitable ones (Hughes 2017). Then the current distribution maps were overlaid on those obtained for all the four future climate scenarios selected to identify any potential changes in range: contraction, expansion, refuge, non-suitable areas (Hatten et al. 2016). Secondly, more restrictive thresholds were set to differentiate three classes of potential habitats (Yang et al. 2013): low potential (10P-0.6), medium potential (0.6-0.8) and high potential (0.8-1). After taking into account these thresholds, the extent of the potential distribution areas (both current and future) were calculated in all the potential zones (see Fig. 1) as the studied species showed a markedly disjunct distribution. For the known localities, the current values of habitat suitability were compared with those obtained under the climate change scenarios.

Current model for Delphinium bolosii
The AUC value was 0.962 ± 0.032 for the D. bolosii model, indicating that the model performed well at predicting the current distribution of the species. The current potential distribution area was 35814.9 km 2 for threshold 10P = 0.4585, of which 14127.1 km 2 (39%) showed a medium-high suitable habitat probability (HS ≥ 0.6). This potential distribution was found mainly in two areas of the north-eastern Iberian Peninsula: the Pre-Pyrenean, along with Catalan Mediterranean System, and the south-eastern Iberian System (Fig. 3).
For the D. bolosii model, the mean temperature of the driest quarter had the strongest relative influence on habitat suitability (82.3%). Other predictors for this model and their contribution were: precipitation of the coldest quarter (5.2%), temperature seasonality (4.7%), precipitation seasonality (4.2%), slope (3.5%) and Isothermality (0.2%). The effects of these environmental variables on habitat suitability for D. bolosii are shown in Fig. 4. For example, the probability of presence was maximum when the mean temperature of the driest quarter went close to 0 °C. The selected variables and its percentage contribution to the D. bolosii and D. fissum subsp. sordidum models are compared in Table 1.

Potential effects of climate change on the series Fissa
Our results showed that potential changes could take place on the distribution range of the series Fissa species, i.e., when comparing the extent and quality of the present and future suitable habitats for D. bolosii (Table 2) and D. fissum subsp. sordidum (Table 3).
For D. bolosii, the total potential distribution area reduced for all the studied climate change scenarios, except for one (2070 RCP 4.5) ( Table 2), mainly in the outer areas of its potential distribution (Fig. 5). The predicted contraction in range for this taxon was 55.5% in the harsher scenario (2070 RCP 8.5) ( Table 4). In contrast, the potential distribution area of D. fissum subsp. sordidum increased in all the climate change scenarios (especially for the harsher scenario). This expansion of potential areas occurred mainly towards inland areas of the northern Sub-Plateau, the eastern Central System, the western Iberian System and the north-eastern area of the Baetic Systems (Fig. 6). Despite this expansion, the models forecast relatively high loss or contraction percentages (47.3% for 2070 RCP 4.5) (Table 4), especially to the east of the northern Sub-Plateau and the Sierra Morena range (Fig. 6). However for both taxa, the potential   area with a medium-high suitability habitat increased. In both cases, habitat quality also significantly increased, i.e. the proportion with a medium-high suitability habitat in relation to the total potential habitat (Tables 2, 3).   Regarding the habitat suitability of both taxa in currently known populations, models forecast some variations under climate change scenarios (Table 5). For Delphinium bolosii, habitat suitability would reduce, generally in the populations located in the Catalan Mediterranean System and the eastern Iberian System, whereas habitat suitability would remain or increase in the Pre-Pyrenean populations. For  D. fissum subsp. sordidum, the main changes would be increased habitat suitability for some populations of the Central System in the RCP 8.5 scenario, and major loss of suitability for almost all the populations in the northern Sub-Plateau. In the Baetic Systems, habitat suitability would remain very high.

Discussion
The models obtained for the series Fissa species reinforce the hypothesis of the differentiation of one eastern taxon (D. bolosii) and another western taxon (D. fissum subsp. sordidum) in the Iberian Peninsula following the migratory patterns proposed by Bocquet et al. (1978). The Iberian System may act as a geographic barrier during this differentiation process. Apparently, the D. bolosii model selected different predictors that the model provided by Rus et al. (2018) for D. fissum subsp. sordidum. However, when we contemplated the differences between the climate features of the north-eastern Iberian Peninsula and the rest of the study area, a remarkable similarity was observed in the effect of these variables on the biological cycle of both taxa. Such similarities, which are discussed later, would be in accordance with the niche conservatism theory (i.e. retention of fundamental niche characteristics over evolutionary time) in phylogenetically closely related taxa (Peterson et al. 1999).
Firstly, the mean temperature of the driest quarter had the strongest influence on the habitat suitability of D. bolosii. In contrast, the variable that contributed the most in the model of D. fissum subsp. sordidum was the mean temperature of the wettest quarter (Rus et al. 2018). The driest quarter in the north-eastern Iberian Peninsula (corresponding to potential distribution of D. bolosii) was winter, and the wettest quarter in the western Iberian Peninsula (corresponding to potential distribution of Delphinium fissum subsp. sordidum) was also winter (Rodríguez-Ballesteros 2015). Thus habitat suitability was maximum for both taxa when the average winter temperature was 0-4 °C. According to Blanché et al. (2014), low temperatures are required to break seed latency in D. bolosii.
Secondly, both models' selected variables were related to the continentality of climate: with D. bolosii it was temperature seasonality, but was temperature annual range for D. fissum subsp. sordidum (Rus et al. 2018). A marked contrast between low winter temperatures and high summer temperatures increased habitat suitability for the two taxa. Both are hemicryptophytes with a rhizome that remains in a latent state during the cold season. Thus subsequently during a period of suitable conditions, the rhizome can develop a rosette of palmatisect leaves and a dense inflorescence with a height that equals 1 m or more.
Both models selected the precipitation of the coldest quarter (i.e. winter season), although the optimal values for D. bolosii (approximately 100 mm) were lower than those indicated for D. fissum subsp. sordidum (approximately 200 mm). This difference can be explained because when D. bolosii restarts its biological cycle at the end of winter, it depends less on winter water reserves than D. fissum subsp. sordidum as spring is generally rainier in the north-eastern than in the western Iberian Peninsula.
Regarding the precipitation of the warmest quarter (i.e. summer season), D. fissum subsp. sordidum requires approximately 50 mm to complete blooming and fruiting stages (Rus et al. 2018). In contrast, the D. bolosii model did not select this variable, despite having a similar biotype and biological cycle. This may be due to summer being relatively rainy from storms in the north-eastern Iberian Peninsula. In such territories, spring and autumn are the rainiest seasons, followed by summer, and finally by winter (Rodríguez-Ballesteros 2015). The habitat suitability for D. bolosii is greater when precipitation seasonality is low. In this context, it should be pointed out that the distribution of rainfall over the year is more important than the annual global rainfall in the distribution of many species (Del Río et al. 2018).
Our models showed potential impacts of climate change on the series Fissa species in the Iberian Peninsula, a result that is consistent with previous climate change analyses conducted with mountain species (Lenoir et al. 2008;Engler et al. 2011) andthreatened taxa (Mendoza-Fernández et al. 2021). Studies about the effects of climate change on biodiversity have predicted range contractions for many species, but range expansion has also been documented (see e.g. Berry et al. 2003;Hamann and Wang 2006). For D. fissum subsp. sordidum, our results indicated that the areas with an oceanic climate in the central-western Iberian Peninsula and mountain ranges with lower elevations, such as Sierra Morena, are more sensitive to changes and, thus, face marked habitat suitability reductions. Conversely, the models forecast range expansion in the areas with a continental climate in the central Iberian Peninsula, as well as the mountain ranges with a bigger scope and higher elevation. For D. bolosii, the outside areas of mountain ranges showed more marked habitat suitability losses given their lower elevation. The habitat suitability loss around the Pyrenees for the most pessimistic scenario (RCP 8.5) was remarkable. Engler et al. (2011) indicated that the flora of the Pyrenees appears especially sensitive to climate change due to increased temperature, but mainly due to reduced precipitation. However, the inner areas of the Eastern Iberian and the Pre-Pyrenean Systems showed good persistence and an increase in potential areas. These results can be explained by the complex meteorological processes indicated by López- Moreno et al. (2008), in which several variables take part, such as wind speed and direction, surface temperature and relative humidity.
Both taxa have limited dispersal ability since seed dispersal occurs by boleochory. These limitations can prevent species from successfully tracking climate to the potential areas of future overlaps Krosby et al. 2015). However, sporadic events of long-distance seed dispersal by herbivores seem to be possible in D. fissum subsp. sordidum (Melendo et al. unpublished data). Anyway, no overlapping was found between either the potential distribution areas in current climatic conditions or the potential areas of the projected models. The possible range expansion of D. fissum subsp. sordidum towards the Iberian System, especially in the most pessimistic scenario, coincided with the retraction of the potential areas of D. bolosii towards the Eastern Iberian System.
From the conservation point of view, loss of habitat suitability seems worrying for some known populations of these two endangered endemic taxa. In particular, the Delphinium fissum subsp. sordidum populations located on the northern Sub-plateau, and some Delphinium bolosii populations located in the Catalan Mediterranean and Eastern Iberian Systems. According to Wiens (2016), climate change-induced contraction in range sizes poses a local extinction threat to many species. However, other factors must be considered. For example, some highly threatened populations have only a few individuals and a high inbreeding rate (Orellana et al. 2007;Bosch et al. 2019; Ramírez-Rodríguez and Amich 2019), while others have many individuals Ramírez-Rodríguez and Amich op. cit.). Some populations are also affected by other threats, such as too many herbivores (Ramírez-Rodríguez et al. 2016;Ramírez-Rodríguez and Amich op. cit.).
During the modelling process of the studied species, biotic variables like vegetation cover, species interactions or dispersal ability were not used. Biotic interactions were implicitly considered. Thus the realized niche of the studied species was modelled, which partly incorporates these interactions. In this way, models were simplified, but provided a static representation of the biotic interactions. This implies a source of uncertainty for models if interactions between species were modified, which could occur with climate change (Suttle et al. 2007). For instance, Lepidopterans are an important group of pollinators for the studied species (Bosch et al. 1998;, which might shift northward as a future response to climate change (Parmesan et al. 1999). Although continuous efforts are being made to explicitly integrate biotic interactions into distribution models, the pitfall of predicting how such interactions will evolve under future environmental conditions still needs to be overcome.
Another source of uncertainty in modelling species distribution stemmed from our inability to predict how species will be able to genetically adapt or express their phenotypic plasticity when faced with changing environmental conditions (Theurillat and Guisan 2001). In this context, long-life species and those species with limited dispersal ability better conserve the realized niche . Indeed, the studied species have these features: a rhizome that can remain for decades and short-distance seed dispersal. Given such a slow migration rate, it is unlikely that the studied species could reach all its recent suitable habitats ). Accordingly, the projections that forecast a large gain in the distribution range might be overestimated.
The changes in vegetation cover associated with climate change might also be another source of uncertainty. However, we think that this might not be relevant because the studied species are able to live in wide-ranging plant communities. However, the selected spatial scale may prove more significant because, in some cases, it might not reflect the more suitable microhabitats for the studied species ).

Conclusions
In the Iberian Peninsula, the series Fissa of the genus Delphinium has diversified into two endemic and endangered taxa, D. bolosii and D. fissum subsp. sordidum, which do not have complete reproductive barriers and need similar climatic requirements. Despite their ecological similarities, climate change would cause different effects in the distribution area of each taxon: while for D.bolosii the total potential area would decrease, it would increase for D. fissum subsp. sordidum. In both cases, the potential distribution area would shift towards areas of higher continentality at present. The aforementioned orographic barriers may play an important role in the maintenance of non-overlapping potential areas. However, some of the populations would face a high risk of local extinction; therefore, monitoring efforts would be conducted for these populations as well as a joint conservation plan which includes in situ and ex situ conservation measures for both taxa.