Next Article in Journal
Default-Mode Network Connectivity Changes Correlate with Attention Deficits in ALL Long-Term Survivors Treated with Radio- and/or Chemotherapy
Previous Article in Journal
Evaluation of the Cardiometabolic Disorders after Spinal Cord Injury in Mice
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

‘Fly to a Safer North’: Distributional Shifts of the Orchid Ophrys insectifera L. Due to Climate Change

by
Martha Charitonidou
1,*,
Konstantinos Kougioumoutzis
2,
Maria Chara Karypidou
3 and
John M. Halley
1
1
Laboratory of Ecology, Department of Biological Applications & Technology, University of Ioannina, 45110 Ioannina, Greece
2
Laboratory of Botany, Department of Biology, University of Patras, 26504 Patras, Greece
3
Department of Meteorology and Climatology, School of Geology, Faculty of Sciences, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece
*
Author to whom correspondence should be addressed.
Biology 2022, 11(4), 497; https://doi.org/10.3390/biology11040497
Submission received: 28 February 2022 / Revised: 15 March 2022 / Accepted: 22 March 2022 / Published: 24 March 2022

Abstract

:

Simple Summary

Climate change is one of the major threats to plant diversity and is expected to force species distributions into latitudinal or altitudinal shifts. The complex biology of orchids, and their many interactions with other organisms, increases their vulnerability in a changing climate. This study focuses on how climatic alterations will affect the distribution of the fly orchid (Ophrys insectifera L.), one of the most well-known and distinctive Ophrys species in Europe, using models that predict the species range changes in the future, based on environmental factors. The orchid’s environmentally suitable area is projected to shift northwards but downhill in the future, experiencing a moderate overall range contraction. More specifically in near- and long-term future, it is expected to be lost in South Europe, especially from the Balkans, while it will gain areas in North Europe, with the UK, Scandinavia, and the Baltic countries being among the winners. These results, although conservative since they are based only on abiotic variables, provide useful insights on the fly orchid’s response to future climatic change, and can serve as a basis for further studies on a finer scale.

Abstract

Numerous orchid species around the world have already been affected by the ongoing climate change, displaying phenological alterations and considerable changes to their distributions. The fly orchid (Ophrys insectifera L.) is a well-known and distinctive Ophrys species in Europe, with a broad distribution across the continent. This study explores the effects of climate change on the range of O. insectifera, using a species distribution models (SDMs) framework that encompasses different climatic models and scenarios for the near- and long-term future. The species’ environmentally suitable area is projected to shift northwards (as expected) but downhill (contrary to usual expectations) in the future. In addition, an overall range contraction is predicted under all investigated combinations of climatic models and scenarios. While this is moderate overall, it includes some regions of severe loss and other areas with major gains. Specifically, O. insectifera is projected to experience major area loss in its southern reaches (the Balkans, Italy and Spain), while it will expand its northern limits to North Europe, with the UK, Scandinavia, and the Baltic countries exhibiting the largest gains.

1. Introduction

During the last few decades, abrupt climatic changes have induced considerable shifts in the spatiotemporal climatic regimes across the globe [1], with strong effects on all levels of biodiversity already being reported [2,3,4,5,6,7]. Several studies have shown that many organisms will or have already experienced significant range shifts due to climate change, with a poleward or uphill direction. Although the majority of species observed to shift their distributions are abiding to this pattern, there are some that migrate southwards or downhill, showing a wide variety of range alterations [4,5,8]. On the other hand, organisms that fail to migrate or adapt may be more vulnerable to extinction [9]. The number of such species is predicted to increase on a global scale [10,11]. However, local extirpations are already widespread for several species [12]. Adding to this, dispersal due to the changing climatic conditions will not help them avoid extinction, even when combined with niche shifts [6].
Climate change has been listed as a potential threat for plants [13]. The main influence comes from changes in climatic factors like temperature and precipitation that play an important role for the plants’ life cycle [2,3], leading to altered phenologies, mismatches on their interactions with other organisms (e.g., pollinators), and distribution area changes [11,13,14,15]. Global warming seems to force plant species to shift their ranges in latitudinal and/or elevation gradients, in search of more favorable climatic conditions [4,16], while in certain cases, it can even act as an amplifier of their vulnerability to extinction [13,17,18,19]. In both hemispheres, plant species tend to migrate towards the poles, while in altitudinal gradient, they shift mostly uphill (e.g., [4,13,16,20,21]), with exceptions of species following opposite directions existing for both patterns (i.e., towards the equator and downhill, respectively (see [4,5]).
Like all plants, orchids are affected by climatic changes, and their response to those alterations has already been studied. Cases of phenological shifts have been reported [7,22,23,24,25] discussing not only changes in orchid flowering patterns, but also the impact of global warming on orchid-pollinator interactions. In addition, numerous studies have been published about the effects of climate change on the distribution of orchids (e.g., [26,27,28,29,30]), reporting different levels of area change, from low to severe range contractions or even expansions. As for the range shifts, orchids seem to follow the general poleward and uphill trends as a general rule (e.g., [31,32], although there may be some contradictory patterns at a local or regional level [33]. Both temperature and precipitation play an important role in driving distribution changes in orchids (e.g., [33,34,35,36,37]). Additionally, orchid species will be affected by climatic fluctuations in these quantities, and climatic models foresee major increases in the variability of both these factors [38,39].
One of the concerns regarding species’ distributions, is that species may be driven either into regions that are unfavorable or very limited in their extent (e.g., [10]). Orchidaceae are one of the families that are under-assessed regarding their vulnerability to extinction [40]. They could face especially significant losses to their distributions and associated population declines [41,42,43], because of the complexity of their life cycle and their symbiotic dependence on other organisms that could act as an amplifier of extinction probability. Consequently, orchids could face greater risks from global warming than in other families [44].
The fly orchid (Ophrys insectifera L.) is one of the most distinctive Ophrys, serving as the type species of the genus. It is characterized by a wide distributional range throughout the continent [45,46], and it can be found in a variety of habitats over a broad altitudinal extent [45,47]. However, in most cases, it occurs in areas with increased soil moisture (damp soils), and always in full sun or semi-shaded sites (see [47,48,49]). Thus, despite its widespread occurrence in Europe, there is cause for concern regarding the conservation of this species in the future. The scope of this paper is to investigate whether such concerns have a reasonable basis. For this purpose, the study focuses on the effects of climate change on the distribution of O. insectifera in the upcoming decades, using a species distribution models (SDMs) framework under different combinations of global circulation models (GCMs), share socioeconomic pathways (SSPs), and future time periods.

2. Materials and Methods

2.1. Study Species

The fly orchid (Ophrys insectifera L.) was first observed and described in the 1750s by Linnaeus, during his trip to the Baltic islands [50], and serves as the type species of the genus Ophrys. It is one of the most distinctive bee orchids, with its lax inflorescence bearing the characteristic fly-resembling flowers, that, in contrast to their appearance, are pollinated by male digger wasps of the genus Argogorytes (A. mystaceus and A. fargeii) [51]. It is a tuberous perennial orchid, with a height varying from (12-)15 to 50 cm, which can be found in bloom during May–July (depending on the region), in damp, full sun/semi-shade sites, in a variety of habitats (from woodlands and forest edges to fens and grasslands), and in a broad altitudinal range (0–1700 m a.s.l.) [45,47]. O. insectifera is a species native to Europe, and in contrast to the mainly Mediterranean range of other bee orchids, it has an extended distribution throughout the continent [45]. Extending from Ireland and Northern Spain to Ukraine (West–East Axis), and from Scandinavia and the Baltic countries to Italy and Northern Greece (North–South Axis) [46], its distribution center lies in Central Europe. However, marginal occurrences can be found in Russia and Northern Norway, thus characterizing O. insectifera as a rather temperate species and the northernmost bee orchid [47] (Figure 1a). Despite its wide range throughout the European continent, its populations can be locally dense or even locally rare, with only scattered occurrences [52].

2.2. Species Occurrence Dataset and Coordinate Thinning Procedures

Occurrence data for Ophrys insectifera were downloaded from the Global Biodiversity Information Facility (GBIF) database (a total of 63,574 occurrences) [54]. In this dataset, available occurrences for Greece (20 occurrences; ‘Orchids of Greece’ project database, Assist. Prof. Spyros Tsiftsis pers. comm.), Serbia (seven occurrences; [55]), Bulgaria (six occurrences; [56]), and Romania (four occurrences; Mihai Bobocea pers. comm.) were added. The original occurrence dataset for O. insectifera derives from the available data from GBIF. However, as shown in Figure 1a, the distribution of the species contains other countries as well (e.g., in Eastern Europe and the Balkans), for which there are no data in the GBIF repository. This issue could question the reliability of the predictions for the species’ range. Nevertheless, it was tackled to some extent with the aforementioned occurrence data additions, and by setting the species distribution area to be equal to the extent of occurrence (EOO) as proposed by guidelines from the International Union for Conservation of Nature (IUCN) [57]. This extent of occurrence was used as a proxy for the species distribution [58] and was calculated using the R-package ‘ConR’ 1.3.0 [59] with the alpha-hull method, adding a buffer around the calculated extent to cover all occurrence data.
All non-georeferenced occurrences were removed from the dataset, as well as any pre-1970 data in order to match the temporal resolution of WorldClim v2.1 baseline data [60]. All points with coordinate uncertainty > 9.25 km were also removed, so as to be in line with the resolution of the selected environmental variables (see Section 2.3). Moreover, the functions ‘clean_coordinates’ of the ‘CoordinateCleaner’ 2.0.18 R package [61] and ‘elimCellDups’ function from the ‘enmSdm’ 0.5.3.3 R package [62] were used to further clean the occurrence dataset and eliminate any duplicate records, respectively.
The remaining occurrence data were spatially thinned using the R package ‘spThin’ 0.1.0 [63] following Robertson et al. [64], resulting in a dataset of 3914 records (hereafter GeoThin). However, according to Varela et al. [65], environmental thinning may lead to improved model performance. To that end, spatially thinned occurrences were further thinned based on the representative and uncorrelated environmental variables occurring in the study area (see Section 2.3), following Varela et al. [65] (and the code provided at: https://github.com/SaraVarela/envSample; last accessed on 26 February 2022). After environmental thinning, the second occurrence dataset consisted of 1802 records (hereafter, EnvThin).

2.3. Environmental Data

Current and future climatic data (minimum temperature, maximum temperature, average temperature, precipitation, and 19 bioclimatic variables) were downloaded from the WorldClim v2.1 database [60] at a 5 arc minutes resolution. Regarding future projections, data from the Coupled Model Intercomparison Project—Phase 6 (CMIP6) were obtained for two time slices, 2070 (timeframe 2061–2080, representing ‘near-term future’) and 2090 (timeframe 2081–2100, representing ‘long-term future’), for three GCMs (BCC-CSM2-MR [66], MIROC-ES2L [67], MRI-ESM2-0 [68]), and four different SSP scenarios [69] (SSP1-2.6, SSP2-4.5, SSP3-7.0, SSP5-8.5). Extra-terrestrial solar radiation and 16 additional bioclimatic variables were constructed for all time-slices using functions from the R package ‘envirem’ 2.2 [70]. Elevation data were downloaded via the CGIAR-CSI data-portal [71], and five topographical variables (slope, aspect, heat load index, topographic position index and terrain ruggedness index) were then estimated, with the R packages ‘raster’ 3.3.13 [72] and ‘spatialEco’ 1.2-0 [73]. Soil variables were extracted from the SoilGrids 2.0 database [74] via Web Coverage Services (WCS) in QGIS v.3.18.0 ‘Zürich’ [75] at a 5 arc minutes resolution. The function ‘vifcor’ from the R package ‘usdm’ 1.1.18 [76] was used to assess multicollinearity. From the initial set of 52 chosen predicting variables, 20 did not show any collinearity problems (Spearman rank correlation < 0.7 and VIF < 10; [77]) and were, thus, included in the analyses (Table S1).

2.4. Species Distribution Models

The realized climatic niche of Ophrys insectifera was modeled by combining the occurrences’ dataset with current environmental predictors in an ensemble modelling scheme, to reduce model algorithm uncertainty [78,79]. This process was followed for both thinning procedures (GeoThin and EnvThin occurrence datasets).
Fine-tuned SDMs were fitted based on four different algorithms: random forests (RF), boosted regression trees (BRT), Bayesian additive regression trees (BART) and maximum entropy (MaxEnt), via the R packages ‘SDMtune’ 1.1.4 [80], ‘embarcadero’ 1.2.0.1003 [81] and ‘ENMeval’ 0.3.1 [82], respectively. Before the procedure of model fitting, spatial cross-validation was applied in order to reduce spatial autocorrelation [83], by spatially partitioning the dataset in four blocks, using the function ‘get.block’ from the R package ‘ENMeval’ 0.3.1 [82]. Regarding the algorithms RF and BRT, models’ hyperparameters were fine-tuned using the functions ‘reduceVar’ and ‘optimizeModel’ of the ‘SDMtune’ R package [80]. There, any variables with low permutation importance (nperm = 100; <5%) were sequentially removed until the optimum Jackknife test TSS was reached, and afterwards, the hyperparameter combination of the best performing model based on test TSS was identified. For BART, fine-tuning process was applied using the functions ‘bart.step’ and ‘retune’ as described in Carlson et al. [84]. For MaxEnt, the combinations between linear, quadratic, and hinge feature classes were explored, since they lead to more comprehensive and better-performing models [85]. The regularization multipliers ranged from 1 to 10 with a step value of 0.1, resulting in 600 candidate models. Optimal MaxEnt model settings were identified based on threshold-dependent (i.e., omission rate) evaluation metrics, in order to prevent overfitting and improve model transferability, as model predictions based on information criteria lead to oversimplified models with low predictive performance [86].
Pseudo-absences (PAs) for O. insectifera were generated following the recommendations of Barbet-Massin et al. [87] and Liu et al. [88]: 30,000 background points were randomly sampled within the study area (defined in Section 2.2), since poor background sampling may lead to a truncated environmental response [89].
The prediction capability of each one of the models was evaluated using a selection of discrimination and calibration metrics. In order to avoid any misleading result by relying to a single metric [90,91,92], four discrimination (AUC, AUC-PR, Sørensen’s index, TSS) and three calibration (Brier score, Cohen’s Kappa, Continuous Boyce Index (CBI)) metrics were selected for the evaluation [93,94,95,96,97,98]. The aforementioned metrics were calculated using functions from the ‘CalibratR’ 0.1.2, ‘DescTools’ 0.99.40, ‘ecospat’ 3.2, ‘enmSdm’ 0.5.3.2, ‘Metrics’ 0.1.4, ‘MLmetrics’ 1.1.1 and ‘modEvA’ 2.0, R packages [62,99,100,101,102,103,104].
Variable importance for each model was estimated via the functions ‘varImp’ (nperm = 1000), ‘varimp’ and ‘var.importance’ from the ‘SDMtune’, ‘embarcadero’ and ‘ENMeval’ R packages, respectively.
The potential distribution of O. insectifera was projected under current and future climate conditions, for all algorithms separately (regarding MaxEnt, the ‘clogclog’ output was used, which is analogous to the occurrence probability predictions from the remaining algorithms [85]), as well as via an ensemble model framework [79], based on calibrated models with TSS ≥ 0.6 (to avoid poorly calibrated ones). The contribution of each model to the ensemble forecast was weighted according to its TSS score. Five ensemble methods were selected (median, mean, weighted mean, committee average and PCA-based), to tackle performance uncertainty [105], and from these the best ensemble model was selected based on TSS and Sorensen’s index [92].
The resulting habitat suitability maps were transformed to binary based on the metric that maximizes the sum of sensitivity and specificity [96,106,107] and the one that maximizes Sorensen’s index as suggested by Leroy et al. [92]. Afterwards, those were compared to the binary maps obtained for each GCM and SSP. As a conservative approach, the suitability of any cells that had non-zero values in the clamping mask was set to NA [108]. Regarding models produced by ‘embarcadero’, the suitability of any cells that had an uncertainty value equal to or higher than the 90% of the cells was set to NA.
Finally, function ‘BIOMOD_RangeSize’ of the R package “biomod2” 3.3.7 [109] was used to assess the projected range change of the species, for the individual and ensemble models for both thinning procedures. O. insectifera, like all orchids, produces numerous dust-like seeds, with extremely low weight, that theoretically can “travel” for long distances [110]. To that end, for this analysis an unlimited dispersal capacity was assumed for the species range.

2.5. Distribution Shift in Latitudinal and Altitudinal Gradient

In order to test if the distribution of Ophrys insectifera shifts latitudinally in the future, the median centroids of distributions of current and future time slices were calculated based on the binarized habitat suitability from the ensembles of both thinning procedures, and for all GCMs/SSPs, using the function ‘st_centroid’ of the R package ‘sf’ 1.0.5 [111]. A distance matrix between current and future centroids was calculated for all cases (near- and long-term future for EnvThin and GeoThin ensemble) from the homonymous toolbox in QGIS v.3.18.0 ‘Zürich’ [75] software. Moreover, to examine the altitudinal shift of the species, elevation data from the CGIAR-CSI data-portal were used in order to extract the mean altitude for the current and future projections, for all combinations of thinning procedures, GCMs and SSPs. Finally, for both types of shifts, non-parametric tests were used, in order to check for any statistically significant difference (functions ‘kruskal.test’ and ‘pairwise.wilcox.test’ of the R-package ‘stats’).

2.6. Response of the Most Important Variable in Locations of Interest

In order to investigate the effect of the most important variable on the distribution of Ophrys insectifera, two locations were selected from the species’ range that represent areas of consistent loss and gain and are congruent under all studied combinations of models and scenarios. Three sets of simulations were utilized for each one of the selected GCMs: the historical simulation (1986–2014), SSP1-2.6 (best-case scenario), and SSP5-8.5 (worst-case scenario) for the period 2015–2100. Simulations for MIROC-ES2L and MRI-ESM2-0 were retrieved from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/; accessed on 31 January 2022) using the “cdsapi” python tool (https://pypi.org/project/cdsapi/; accessed on 31 January 2022), while BCC-CSM2-MR simulations were retrieved from the Earth System Grid Federation (https://esgf-data.dkrz.de/projects/esgf-dkrz/; accessed on 31 January 2022) using bash scripts. Data were analyzed in their native horizontal resolution, and the required variable is calculated from the original data, following the definition by [58], for every year during the period 1986–2100 for each GCM and SSP. Time series were extracted without smoothing, based on the grid point that was the nearest neighbor to each one of the two selected locations. The processing of the raw files was performed using Climate Data Operators (https://code.mpimet.mpg.de/projects/cdo/; accessed on 31 January 2022), while data were plotted using the R package ‘ggplot2’ [112].

3. Results

All models—apart from MaxEnt—performed sufficiently well (TSS: 0.747 ± 0.232; AUC: 0.895 ± 0.120; PRAUC: 0.782 ± 0.317; Cohen’s Kappa: 0.598 ± 0.318; CBI: 0.994 ± 0.009; Sorensen’s index: 0.794 ± 0.011; Brier Score: 0.059 ± 0.049; see Table S2). Ensemble models performed equally well (median AUC: 0.994 ± 0.003; Brier score: 0.000 ± 0.000; CBI: 0.978 ± 0.043; Sorensen’s index: 0.774 ± 0.052; TSS: 0.947 ± 0.028; see Table S3).
Among all the response variables, temperature seasonality had the highest contribution for the majority of combinations of thinning procedures and algorithms, followed by the Thornthwaite Aridity Index. Exceptions to this were the case of EnvThin-BART, where Thornthwaite Aridity Index was the most important variable followed by Temperature seasonality, and the ones of EnvThin-RF, GeoThin-RF and GeoThin-BRT, where elevation and Precipitation Seasonality were the second most important variables (Table S4).
The resulting habitat suitability maps that were similar for the ensemble of both thinning procedures (Figure 2 and Figure S2) were converted into binary maps, and then compared to the binary maps obtained for each GCM, SSP scenario, time-slice, and thinning procedure. Since the outcome of the ensemble modeling for the current and future potential distribution of Ophrys insectifera was largely similar across all combinations of thinning procedures and future scenarios, results for the EnvThin ensemble were selected to be shown in-text. In addition, for future predictions, only SSP1-2.6 and SSP5-8.5 are presented in-text, as best-case and worst-case scenario, respectively.
Ophrys insectifera is projected to experience moderate overall contraction of its environmentally suitable area in the near-term future (mean contraction in the worst-case scenario of the presented EnvThin ensemble: 15.59%—2070; mean current occurrence lost for the worst-case scenario of the presented EnvThin ensemble: 27.65%—2070; Table 1). However, the species is expected to face large changes across its distribution. A considerable area loss is observed (mean expected loss: 38.48% for the worst-case future scenario of EnvThin ensemble in 2070; Table 1), mainly at the southern parts of the species’ distribution, as well as at its west and east edges (Figure 3). O. insectifera is also anticipated to show an increase of environmentally suitable areas in the near future (mean expected gain: 22.9% for the worst-case scenario of EnvThin ensemble in 2070; Table 1). This gain is observed mainly at the leading edge of the distribution, with some sporadic suitable areas’ increase at regions on both sides and close to the distribution centre (Figure 3). Results for long-term future (2090) show similar trends for the range changes of the species. Contraction of environmentally suitable areas is expected to be higher than in 2070 (mean contraction for the EnvThin ensemble in 2090: 27.51%; mean current occurrence lost for the worst-case scenario of the presented EnvThin ensemble: 40.12%—2090; Table 1). Loss of currently occupied areas is expected to be larger than in 2070 (mean expected loss: 48.87% for the worst-case scenario of EnvThin ensemble in 2090; Table 1), and observed at the trailing edge of the distribution in the south (Figure S2). Area gain is slightly lower than in 2070 (21.35%; Table 1) and located in the same areas at the northern parts of the distribution (Figure S2). The results are consistent throughout all combinations of GCMs, SSPs and future time slices, for both EnvThin and GeoThin ensemble models (see Figures S3–S8 and Tables S5 and S6).
The northwards shift of the species can be observed more easily by looking at the centroid movement. The centroids of all projected future environmentally suitable areas of Ophrys insectifera for the selected ensemble appear to be always lying to the north and mainly northwest of the current centroid, with very few exceptions, e.g., BC85 centroid in the presented ensemble, that is placed to the northeast of the current (Figure 4 and Figure S10). The distance between the current and future centroids depended on the combination of GCM and SSP, and in the presented EnvThin ensemble for 2070 varied from 59.6 to 323.5 km (MI26 and BC85 respectively; Table S9). For 2090, the centroids for all combinations of GCMs and SSPs show similar patterns (Figures S9 and S11), with the distance between the current and future centroids being even larger (85.7 to 759.2 km, for MI26 and BC85 respectively; Table S9). Results are consistent among all combinations of GCMs, SSPs and future time slices for both EnvThin and GeoThin ensemble models (Figures S12–S15 and Table S10).
The results for the altitudinal shift test showed a moderate downhill movement for all different combinations of GCMs and SSPs of future projections for both 2070 and 2090 (Table 2 and Table S7; 2070 Kruskal-Wallis test x2 = 863.38, df = 6, p < 0.001; 2090 Kruskal-Wallis test x2 = 3504.2, df = 6, p < 0.001). The species can currently be found at a mean altitude of 542.6 m a.s.l. for the presented ensemble, while in future it is expected to be found in areas with a mean altitude lower by ca. 100–140 m for the worst-case scenario of 2070 and 2090 respectively. Similar results have been observed for all combinations of GCMs, SSPs and future time slices for the GeoThin ensemble (Table S8).
Figure 5 shows the time series for temperature variability (yielded as the most important variable in SDMs) for the two selected locations, Greece and the UK, at the southern and northern reach of Ophrys insectifera, respectively. As mentioned in Section 2.6, these two locations were selected as examples for comparison purposes, since they consistently exhibit high rates of area loss (Greece) and gain (UK) in all future projections. Temperature seasonality in Greece displays a consistent increase with a high intra-annual variability for all GCMs, which is maximized for the worst-case scenario. On the other hand, GCMs for the UK exhibit a considerable consensus compared to Greece. Over this location, temperature seasonality also shows an increasing trend, however the magnitudes of change are significantly smaller than in Greece.

4. Discussion

On the basis of results obtained in this study, the environmentally suitable area of Ophrys insectifera is expected to shift northwards. Its future centroid will shift by 59.6 to 323.5 km (EnvThin ensemble’s MIROC SSP1-2.6 and BCC SSP5-8.5, respectively; Table S9) from its current location, which lies close to the tripoint of France, Germany and Switzerland (Figure 4). This northward shift is robust under all investigated climatic models, scenarios, and future time frames (Figures S9–S15). By contrast, there is no consistent direction, westward or eastward in the forecasts for the environmentally suitable areas. This migration pattern of the fly orchid is in line with the prevailing trend anticipated for plants in general (see [4,13,16]), and for other orchid species in Europe, including Cypripedium calceolus [32], Epipactis helleborine [35], Himantoglossum hircinum [31,113], Orchis militaris, O. simia, O. anthropophora and O. purpurea [26].
Also expected is a moderate downhill movement (under all investigated climatic models, scenarios and future time frames; Tables S7 and S8). In the presented EnvThin ensemble model O. insectifera shows a mean altitudinal difference of 115 m by 2070 and of 265 m by 2090 (Table 2). Although this downhill shift goes against the common uphill pattern, it has been observed before. Lenoir et al. [5] found that 30.9% of forest plant species in Europe are expected to move to lower elevation in the future, while Chen et al. [4] showed that species exhibit high variation of shifts across groups and regions, with ca. 25% of the studied taxa moving downhill, in contrast to the mainstream uphill trend. This may be due to a complex interaction between the altitudinal and latitudinal temperature gradients, that strongly affect the niche patterns in areas of extended elevation shift.
There have been long-term observations of orchid species in Europe, which report decline in numbers, as well as a systematic contraction of the range in certain European orchids [114,115]. This contraction is expected to intensify due to climatic change. For example, Epipactis helleborine is expected to face a decrease of 25–40% of its suitable habitat areas until 2080 [35], Traunsteinera globosa is projected to show a future range contraction of 18–32% [29], and Orchis militaris, that is anticipated to decrease by up to 61.3% in the future worst-case scenario [26]. O. insectifera follows the same trend as other orchids with temperate extent, showing a moderate contraction pattern, for all studied models and scenarios, which will reach 27.5% in the worst-case scenario in the long-term future.
According to Evans et al. [26], as climate continues to change, areas in Northern Europe will become more suitable for orchid species, in contrast to southern European countries that are expected to experience more intense climatic alterations (see also [1]). The aforementioned future changes in the suitable areas of Ophrys insectifera are mainly driven by temperature variability, and aridity (see Table S4). Temperature is known to influence orchid distributions and abundance (e.g., [35,37,116]). Although orchids can be tolerant to temperature seasonal changes, large fluctuations can have negative effects on the life cycle and the species’ population dynamics [117]. Soil moisture is also a factor that strongly influences orchid distributions, through its effect on the below-ground stages of the orchids’ life cycle [118]. However, soil moisture, being a factor associated with microclimate [117], cannot be easily included in a large-scale analysis, and thus is absent from available databases. Nevertheless, aridity can serve as a proxy for soil moisture, and thus the distribution changes of the species exhibit considerable dependence on this factor.
Despite the overall moderate change predicted in this study, O. insectifera will exhibit strong patterns of environmentally suitable area loss and gain in certain locations. More specifically, the species will show high percentages of loss at the southern reach of its current extent (Figure 3 and Figures S2–S8). This includes Spain, Italy, and countries of the Balkan Peninsula such as Bulgaria and Greece, from which it will disappear entirely. On the other hand, the northern reach of the species distribution will see large gains of suitable areas. Among the regions projected to show gain in the north of Europe, the Baltic countries, Scandinavia and especially the UK appear as gain hotspots (Figure 3 and Figures S2–S8). O. insectifera is the only Ophrys with a European temperate distribution that is found in a wide altitudinal range across its extent, but always in sites with increased water presence [45,48], in contrast to other bee orchids, that are adapted to drier conditions occurring in Mediterranean habitats [47]. Larger fluctuations and high increase of temperature seasonality observed in the south of Europe are in contrast to the less intense ones to the north (Figure 5). This, combined with the vast increase of aridity, can explain the projected changes of the species in future decades. This result is also in line with the findings of Thuiller et al. [18], that found high percentages of plant species’ loss in the southern parts of Europe (up to 62.1% in Mediterranean mountainous areas), and high percentages of species’ turnover in Northern Europe (Central Atlantic, Continental and Boreal biogeographic regions).
The anticipated gain of area to the north reach of O. insectifera seen in this study is in contrast to trends observed so far. Declines of the species have been observed in the Netherlands, Flanders, Estonia and the UK [114,119,120,121], with an exception coming from Denmark, where the species has been showing an increase the last 30 years [122]. However, such declines from the past to present, probably reflect previous land use change (see [123]). Although the key driver of declines in the future is likely to be climatic, land use change will also play an important role.
On the other hand, the observed loss in southern Europe under all investigated cases can also raise concerns about the conservation of the species’ diversity. More specifically, at the southwest border of O. insectifera distribution, there have been described two endemic subspecies, O.insectifera subsp. subinsectifera in the Spanish Pyrenees, and O. insectifera subsp. aymoninii in the French Causses region [46,47]. Although they differentiate in a phenological aspect, their genetic differentiation is rather weak [124]. Additionally, at the southeastern limits of the species’ distribution, Greek and Bulgarian populations of O. insectifera were genetically distant from others, and were characterized by a unique haplotype, underlining a cryptic diversity within the species [124]. According to this study’s results, areas that currently host this diversity at the trailing edge of the species are projected to be lost in future time slices, under all combinations of models and scenarios of climate change, pointing to an urgent need for conservation of genetic resources.
Despite the favorable performance metrics and the congruence of results between all combinations of models and scenarios, the estimated predictions for O. insectifera could be characterized as conservative, and potentially underestimating the actual trends. Among the limitations of SDMs, is the difficulty in modelling the effect of interaction between species, population dynamics and evolution. Orchids are known to depend on symbiotic relationships with other organisms above and below ground for their survival [125,126]. O. insectifera is a specialist that relies on two specific digger wasp species (Argogorytes spp.) for its pollination [48]. Below-ground, it also relies on interactions with fungi of Tulasnellaceae (at the mature stage) and Ceratobasidiaceae (at the protocorm stage) [127]. These interactions add a level of complexity when dealing with climatic changes; in order for a species to shift its distribution without further adaptations, so should its close partners. An alternative outcome in such a situation is the evolution of new or substantially modified interactions. Evolution driven by climatic change has been stated before for the genus Ophrys, that during the Pleistocene exhibited high diversification rates and shifted from wasp- to bee-pollinated [128]. However, it should be noted that this past evolution due to climate change occurred over a much longer time scale compared to the current climatic alterations, where many species will fail to cope with its comparatively rapid pace [19]. As for their population dynamics, orchids have shown great variability of intra- and interspecific growth rates [126], and at the same time show a degree of resilience to environmental changes. Nevertheless, extremes of environmental variability can result to more dramatic population declines, which, combined with orchids’ massive fecundity, can lead to population extirpations [129] and further range changes of orchid taxa. However, such effects are beyond the intended scope of this paper.

5. Conclusions

This study explored the effects of climate change on the distribution of the fly orchid (Ophrys insectifera L.) in the near- and long-term future, by means of a species distribution models (SDMs) framework, under different combinations of global circulation models (GCMs), share socioeconomic pathways (SSPs), and three different future time periods. According to our findings, O. insectifera is expected to face moderate contractions of its environmentally suitable area in the future, with a projected loss at the warm edge of its distribution under all possible cases examined. However, a possible expansion of its environmentally suitable area is observed on its leading edge, while the species is anticipated to move northwards and downhill to cope with the changing environmental conditions.
These results provide a basis for further investigation of the effects of the changing climate on this particular orchid. Additional studies should focus on combining outcome of SDMs with dynamic population models, especially for populations that are anticipated to face major climatic changes, as well as those at the leading edges of the species’ distribution, where gain is expected. Of special importance are issues involving land-use change, since such factors are believed to be behind widely reported declines in orchid abundance in Europe, which contrasts with the increases expected in northern Europe due to climatic shifts. The effect of land use change is also listed among the factors that strongly affect the distributions of orchid species (see [123]), hence, further work should focus on how changes in land use can affect the distribution of Ophrys insectifera in future decades.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/biology11040497/s1, Table S1: Predictors that did not show any collinearity issues and were included in the analysis; Table S2: Discrimination and calibration metrics for each one of the separate models (RF: random forest; BRT: boosted regression trees; BART: Bayesian additive regression trees; MaxEnt: maximum entropy), and for both occurrence data thinning procedures; Table S3: Evaluation metrics for the five selected methods for each ensemble model, and both occurrence data thinning procedures; Table S4: Five first most important variables for each of the separate models and for both occurrence data thinning procedures; Table S5: Percentages of area loss, gain, overall area change and current occurrences lost for future projections of Ophrys insectifera of the EnvThin ensemble for all selected GCMs and two SSPs (SSP2-4.5, SSP3-7.0); Table S6: Percentages of area loss, gain, overall area change and current occurrences lost for future projections of Ophrys insectifera of the GeoThin ensemble for all selected GCMs and all SSPs; Table S7: Mean altitude for present and future projections of Ophrys insectifera distribution for the EnvThin ensemble model for all selected GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Table S8: Mean altitude for present and future projections of Ophrys insectifera distribution for the GeoThin ensemble model for all selected GCMs and SSPs; Table S9: Distance matrices between current and future distribution centroids for the EnvThin ensemble, for all GCMs and SSPs; Table S10: Distance matrices between current and future distribution centroids for the GeoThin ensemble, for all GCMs and SSPs; Figure S1: Current habitat suitability map for Ophrys insectifera, for the ensemble model using the geographical thinning procedure; Figure S2: Current and future (2090) potential distribution maps for Ophrys insectifera EnvThin ensemble for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5) as ‘best’- and ‘worst’-case scenario; Figure S3: Current and future (2070) potential distribution maps for Ophrys insectifera EnvThin ensemble model for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S4: Current and future (2090) potential distribution maps for Ophrys insectifera EnvThin ensemble model for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S5: Current and future (2070) potential distribution maps for Ophrys insectifera GeoThin ensemble model for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5) as ‘best’ and ‘worst’ case scenario; Figure S6: Current and future (2090) potential distribution maps for Ophrys insectifera GeoThin ensemble model for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5) as ‘best’ and ‘worst’ case scenario; Figure S7: Current and future (2070) potential distribution maps for Ophrys insectifera GeoThin ensemble model for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S8: Current and future (2090) potential distribution maps for Ophrys insectifera GeoThin ensemble model for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S9: Median centroids for current and future (2090) projected distributions of Ophrys insectifera, for the EnvThin ensemble for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5); Figure S10: Median centroids for current and future (2070) projected distributions of Ophrys insectifera, for the EnvThin ensemble for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S11: Median centroids for current and future (2090) projected distributions of Ophrys insectifera, for the EnvThin ensemble for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S12: Median centroids for current and future (2070) projected distributions of Ophrys insectifera, for the GeoThin ensemble for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5); Figure S13: Median centroids for current and future (2090) projected distributions of Ophrys insectifera, for the GeoThin ensemble for all GCMs and two SSPs (SSP1-2.6 and SSP5-8.5); Figure S14: Median centroids for current and future (2070) projected distributions of Ophrys insectifera, for the GeoThin ensemble for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0); Figure S15: Median centroids for current and future (2090) projected distributions of Ophrys insectifera, for the GeoThin ensemble for all GCMs and two SSPs (SSP2-4.5 and SSP3-7.0).

Author Contributions

Conceptualization, M.C. and K.K.; data curation, M.C., K.K. and M.C.K.; formal analysis, M.C., K.K. and M.C.K.; funding acquisition, M.C.; investigation, M.C.; methodology, K.K., M.C. and M.C.K.; resources, M.C., J.M.H. and M.C.K.; supervision, J.M.H.; validation, M.C., K.K. and J.M.H.; visualization, M.C.; writing—original draft, M.C., K.K., J.M.H. and M.C.K.; writing—review and editing, M.C., J.M.H., K.K. and M.C. All authors have read and agreed to the published version of the manuscript.

Funding

This study is part of first author’s doctoral research, funded by the General Secretariat for Research & Technology (GSRT) and the Hellenic Foundation for Research & Innovation (HFRI) (Scholarship Code: 2359).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank Spyros Tsiftsis for kindly providing Ophrys insectifera data for Greece from the ‘Orchids of Greece’ project database, and Mihai Bobocea for O. insectifera data for Romania from his personal observations. The authors would also like to thank the Scientific Support Centre of the Aristotle University of Thessaloniki (Greece) for providing computational/storage infrastructure and technical support for the analysis of temperature seasonality.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or in the decision to publish the results.

References

  1. IPCC. Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M.I., et al., Eds.; Cambridge University Press: Cambridge, UK, 2022; In Press. [Google Scholar]
  2. Walther, G.-R.; Post, E.; Convey, P.; Menzel, A.; Parmesan, C.; Beebee, T.J.C.; Fromentin, J.-M.; Hoegh-Guldberg, O.; Bairlein, F. Ecological responses to recent climate change. Nature 2002, 416, 389–395. [Google Scholar] [CrossRef] [PubMed]
  3. Parmesan, C. Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol. Syst. 2006, 37, 637–669. [Google Scholar] [CrossRef] [Green Version]
  4. Chen, I.C.; Hill, J.K.; Ohlemüller, R.; Roy, D.B.; Thomas, C.D. Rapid range shifts of species associated with high levels of climate warming. Science 2011, 333, 1024–1026. [Google Scholar] [CrossRef] [PubMed]
  5. Lenoir, J.; Gégout, J.C.; Marquet, P.A.; De Ruffray, P.; Brisse, H. A significant upward shift in plant species optimum elevation during the 20th century. Science 2008, 320, 1768–1771. [Google Scholar] [CrossRef] [PubMed]
  6. Román-Palacios, C.; Wiens, J.J. Recent responses to climate change reveal the drivers of species extinction and survival. Proc. Natl. Acad. Sci. USA 2020, 117, 4211–4217. [Google Scholar] [CrossRef] [PubMed]
  7. Hoffmann, A.A.; Rymer, P.D.; Byrne, M.; Ruthrof, K.X.; Whinam, J.; McGeoch, M.; Bergstrom, D.M.; Guerin, G.R.; Sparrow, B.; Joseph, L.; et al. Impacts of recent climate change on terrestrial flora and fauna: Some emerging Australian examples. Austral. Ecol. 2019, 44, 3–27. [Google Scholar] [CrossRef] [Green Version]
  8. Freeman, B.G.; Lee-Yaw, J.A.; Sunday, J.M.; Hargreaves, A.L. Expanding, shifting and shrinking: The impact of global warming on species’ elevational distributions. Glob. Ecol. Biogeogr. 2018, 27, 1268–1276. [Google Scholar] [CrossRef]
  9. Corlett, R.T.; Westcott, D.A. Will plant movements keep up with climate change? Trends Ecol. Evol. 2013, 28, 482–488. [Google Scholar] [CrossRef] [PubMed]
  10. Thomas, C.D.; Cameron, A.; Green, R.E.; Bakkenes, M.; Beaumont, L.J.; Collingham, Y.C.; Erasmus, B.F.N.; de Siqueira, M.F.; Grainger, A.; Hannah, L.; et al. Extinction risk from climate change. Nature 2004, 427, 145–148. [Google Scholar] [CrossRef]
  11. Urban, M.C. Accelerating extinction risk from climate change. Science 2015, 348, 571–573. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  12. Wiens, J.J. Climate-Related Local Extinctions Are Already Widespread among Plant and Animal Species. PLoS Biol. 2016, 14, e2001104. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Parmesan, C.; Hanley, M.E. Plants and climate change: Complexities and surprises. Ann. Bot. 2015, 116, 849–864. [Google Scholar] [CrossRef] [PubMed]
  14. Easterling, D.R.; Meehl, G.A.; Parmesan, C.; Changnon, S.A.; Karl, T.R.; Mearns, L.O. Climate Extremes: Observations, Modeling, and Impacts. Science 2000, 289, 2068–2074. [Google Scholar] [CrossRef] [Green Version]
  15. Menzel, A.; Sparks, T.H.; Estrella, N.; Koch, E.; Aaasa, A.; Ahas, R.; Alm-Kübler, K.; Bissolli, P.; Braslavská, O.; Briede, A.; et al. European phenological response to climate change matches the warming pattern. Glob. Chang. Biol. 2006, 12, 1969–1976. [Google Scholar] [CrossRef]
  16. Parmesan, C.; Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 2003, 421, 37–42. [Google Scholar] [CrossRef] [PubMed]
  17. Nic Lughadha, E.; Bachman, S.P.; Leão, T.C.C.; Forest, F.; Halley, J.M.; Moat, J.; Acedo, C.; Bacon, K.L.; Brewer, R.F.A.; Gâteblé, G.; et al. Extinction risk and threats to plants and fungi. Plants People Planet 2020, 2, 389–408. [Google Scholar] [CrossRef]
  18. Thuiller, W.; Lavorel, S.; Araújo, M.B.; Sykes, M.T.; Prentice, I.C.; Thuiller, W.; Lavorel, S.; Arau, M.B. Climate change threats to plant diversity in Europe. Proc. Natl. Acad. Sci. USA 2005, 102, 8245–8250. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Duchenne, F.; Martin, G.; Porcher, E. European plants lagging behind climate change pay a climatic debt in the North, but are favoured in the South. Ecol. Lett. 2021, 24, 1178–1186. [Google Scholar] [CrossRef] [PubMed]
  20. Bellard, C.; Bertelsmeier, C.; Leadley, P.; Thuiller, W.; Courchamp, F. Impacts of climate change on the future of biodiversity. Ecol. Lett. 2012, 15, 365–377. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  21. Calinger, K.M.; Queenborough, S.; Curtis, P.S. Herbarium specimens reveal the footprint of climate change on flowering trends across north-central North America. Ecol. Lett. 2013, 16, 1037–1044. [Google Scholar] [CrossRef] [Green Version]
  22. Robbirt, K.M.; Davy, A.J.; Hutchings, M.J.; Roberts, D.L. Validation of biological collections as a source of phenological data for use in climate change studies: A case study with the orchid Ophrys sphegodes. J. Ecol. 2011, 99, 235–241. [Google Scholar] [CrossRef]
  23. Robbirt, K.M.; Roberts, D.L.; Hutchings, M.J.; Davy, A.J. Potential disruption of pollination in a sexually deceptive orchid by climatic change. Curr. Biol. 2014, 24, 2845–2849. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Hutchings, M.J.; Robbirt, K.M.; Roberts, D.L.; Davy, A.J. Vulnerability of a specialized pollination mechanism to climate change revealed by a 356-year analysis. Bot. J. Linn. Soc. 2018, 186, 498–509. [Google Scholar] [CrossRef]
  25. Pfeifer, M.; Heinrich, W.; Jetschke, G. Climate, size and flowering history determine flowering pattern of an orchid. Bot. J. Linn. Soc. 2006, 151, 511–526. [Google Scholar] [CrossRef] [Green Version]
  26. Evans, A.; Janssens, S.; Jacquemyn, H. Impact of climate change on the distribution of four closely related Orchis (Orchidaceae) species. Diversity 2020, 12, 312. [Google Scholar] [CrossRef]
  27. Geppert, C.; Perazza, G.; Wilson, R.J.; Bertolli, A.; Prosser, F.; Melchiori, G.; Marini, L. Consistent population declines but idiosyncratic range shifts in Alpine orchids under global change. Nat. Commun. 2020, 11, 1–11. [Google Scholar] [CrossRef]
  28. Konowalik, K.; Kolanowska, M. Climatic niche shift and possible future spread of the invasive South African Orchid Disa bracteata in Australia and adjacent areas. PeerJ 2018, 6, e6107. [Google Scholar] [CrossRef] [Green Version]
  29. Kolanowska, M. The future of a montane orchid species and the impact of climate change on the distribution of its pollinators and magnet species. Glob. Ecol. Conserv. 2021, 32, e01939. [Google Scholar] [CrossRef]
  30. Foster, E.A.; Ackerman, J.D. Future changes in the distribution of two non-indigenous orchids and their acquired enemy in Puerto Rico. Biol. Invasions 2021, 23, 3545–3563. [Google Scholar] [CrossRef]
  31. van der Meer, S.; Jacquemyn, H.; Carey, P.D.; Jongejans, E. Recent range expansion of a terrestrial orchid corresponds with climate-driven variation in its population dynamics. Oecologia 2016, 181, 435–448. [Google Scholar] [CrossRef]
  32. Kolanowska, M.; Jakubska-Busse, A. Is the lady’s-slipper orchid (Cypripedium calceolus) likely to shortly become extinct in Europe?—Insights based on ecological niche modelling. PLoS ONE 2020, 15, e0228420. [Google Scholar] [CrossRef]
  33. Charitonidou, M.; Kougioumoutzis, K.; Halley, J.M. An Orchid in Retrograde: Climate-Driven Range Shift Patterns of Ophrys helenae in Greece. Plants 2021, 10, 470. [Google Scholar] [CrossRef] [PubMed]
  34. Tsiftsis, S.; Djordjević, V. Modelling sexually deceptive orchid species distributions under future climates: The importance of plant–pollinator interactions. Sci. Rep. 2020, 10, 1–12. [Google Scholar] [CrossRef]
  35. Kolanowska, M. Niche Conservatism and the Future Potential Range of Epipactis helleborine (Orchidaceae). PLoS ONE 2013, 8, e77352. [Google Scholar] [CrossRef] [Green Version]
  36. Kolanowska, M.; Kras, M.; Lipińska, M.; Mystkowska, K.; Szlachetko, D.L.; Naczk, A.M. Global warming not so harmful for all plants-response of holomycotrophic orchid species for the future climate change. Sci. Rep. 2017, 7, 1–13. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Ongaro, S.; Martellos, S.; Bacaro, G.; De Agostini, A.; Cogoni, A.; Cortis, P. Distributional pattern of sardinian orchids under a climate change scenario. Community Ecol. 2018, 19, 223–232. [Google Scholar] [CrossRef]
  38. Cos, J.; Doblas-Reyes, F.; Jury, M.; Marcos, R.; Bretonnière, P.; Samsó, M. The Mediterranean climate change hotspot in the CMIP5 and CMIP6 projections. Earth Syst. Dyn. 2022, 13, 321–340. [Google Scholar] [CrossRef]
  39. Tebaldi, C.; Debeire, K.; Eyring, V.; Fischer, E.; Fyfe, J.; Friedlingstein, P.; Knutti, R.; Lowe, J.; Neill, B.O.; Sanderson, B. Climate model projections from the Scenario Model Intercomparison Project (ScenarioMIP) of CMIP6. Earth Syst. Dyn. 2021, 12, 253–293. [Google Scholar] [CrossRef]
  40. Bachman, S.P.; Nic Lughadha, E.M.; Rivers, M.C. Quantifying progress toward a conservation assessment for all plants. Conserv. Biol. 2018, 32, 516–524. [Google Scholar] [CrossRef] [Green Version]
  41. Seaton, P.T.; Hu, H.; Perner, H.; Pritchard, H.W. Ex Situ Conservation of Orchids in a Warming World. Bot. Rev. 2010, 76, 193–203. [Google Scholar] [CrossRef]
  42. Kull, T.; Selgis, U.; Peciña, M.V.; Metsare, M.; Ilves, A.; Tali, K.; Sepp, K.; Kull, K.; Shefferson, R.P. Factors influencing IUCN threat levels to orchids across Europe on the basis of national red lists. Ecol. Evol. 2016, 6, 6245–6265. [Google Scholar] [CrossRef] [PubMed]
  43. Fay, M.F.; Pailler, T.; Dixon, K.W. Orchid conservation: Making the links. Ann. Bot. 2015, 116, 377–379. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Fay, M.F. Orchid conservation: How can we meet the challenges in the twenty-first century? Bot. Stud. 2018, 59, 1–6. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Delforge, P. Orchids of Europe, North Africa and the Middle East, 3rd ed.; A&C Black: London, UK, 2006; ISBN 071367525X. [Google Scholar]
  46. Kühn, R.; Pedersen, H.; Cribb, P.J. Field Guide to the Orchids of Europe and the Mediterranean; Royal Botanic Gardens, Kew: Richmond, UK, 2019; ISBN 978-1842466698. [Google Scholar]
  47. Pedersen, H.; Faurholdt, N. Ophrys: The bee orchids of Europe; Kew Publishing: Richmond, UK, 2007; ISBN 9781842461525. [Google Scholar]
  48. Stroh, P.A. Ophrys insectifera L. Fly Orchid Species Account; Botanical Society of Britain and Ireland: Durham, UK, 2015. [Google Scholar]
  49. Tsiftsis, S.; Antonopoulos, Z. Atlas of the Greek Orchids; Mediterraneo Editions: Rethymnon, Greece, 2017; Volume 1–2, ISBN 9789606848933. [Google Scholar]
  50. Fay, M.F.; Taylor, I.; Sayers, B. 804. Ophrys insectifera. Curtis’s Bot. Mag. 2015, 32, 51–62. [Google Scholar] [CrossRef]
  51. Borg-Karlson, A.K.; Groth, I.; Ågren, L.; Kullenberg, B. Form-specific fragances from Ophrys insectifera L. (Orchidaceae) attract species of different pollinator genera. Evidence of sympatric speciation? Chemoecology 1993, 4, 39–45. [Google Scholar] [CrossRef]
  52. Rankou, H. Ophrys insectifera. IUCN Red List Threat. Species 2011, e.T175957A7153465. [Google Scholar] [CrossRef]
  53. IUCN European Policy Plants 2011. Ophrys insectifera . Available online: https://www.iucnredlist.org/ (accessed on 8 March 2021).
  54. GBIF.org. Gbif Occurrence Download. Available online: https://doi.org/10.15468/dl.6ykdxh (accessed on 18 January 2021).
  55. Djordjević, V.; Lakušić, D.; Jovanović, S.; Stevanović, V. Distribution and conservation status of some rare and threatened orchid taxa in the central Balkans and the southern part of the Pannonian plain. Wulfenia 2017, 24, 143–162. [Google Scholar]
  56. Popatanasov, A. Ophrys insectifera L.—Update of the Status of a Critically Endangered Orchid in Bulgaria. J. Life Sci. 2018, 12, 83–91. [Google Scholar] [CrossRef] [Green Version]
  57. IUCN Standards and Petitions Committee. Guidelines for Using the IUCN Red List Categories and Criteria, Version 14; 2019, Prepared by the Standards and Petitions Committee. Available online: https://www.iucnredlist.org/documents/RedListGuidelines.pdf (accessed on 22 May 2021).
  58. Joppa, L.N.; Butchart, S.H.M.; Hoffmann, M.; Bachman, S.P.; Akçakaya, H.R.; Moat, J.F.; Böhm, M.; Holland, R.A.; Newton, A.; Polidoro, B.; et al. Impact of Alternative Metrics on Estimates of Extent of Occurrence for Extinction Risk Assessment. Conserv. Biol. 2016, 30, 362–370. [Google Scholar] [CrossRef] [PubMed]
  59. Dauby, G.; Stévart, T.; Droissart, V.; Cosiaux, A.; Deblauwe, V.; Simo-Droissart, M.; Sosef, M.S.M.; Lowry, P.P.; Schatz, G.E.; Gereau, R.E.; et al. ConR: An R Package to Assist Large-Scale Multispecies Preliminary Conservation Assessments Using Distribution Data. Ecol. Evol. 2017, 7, 11292–11303. [Google Scholar] [CrossRef] [Green Version]
  60. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 2017, 37, 4302–4315. [Google Scholar] [CrossRef]
  61. Zizka, A.; Silvestro, D.; Andermann, T.; Azevedo, J.; Duarte Ritter, C.; Edler, D.; Farooq, H.; Herdean, A.; Ariza, M.; Scharn, R.; et al. “CoordinateCleaner”: Standardized cleaning of occurrence records from biological collection databases. Methods Ecol. Evol. 2019, 10, 744–751. [Google Scholar] [CrossRef]
  62. Smith, A.B. Enmsdm: Tools for Modeling Species Niches and Distributions; R Package Version 0.5.1.5; CRAN, R Core Team: Cary, CA, USA, 2020. [Google Scholar]
  63. Aiello-Lammens, M.E.; Boria, R.A.; Radosavljevic, A.; Vilela, B.; Anderson, R.P. spThin: An R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography Cop. 2015, 38, 541–545. [Google Scholar] [CrossRef]
  64. Robertson, M.P.; Visser, V.; Hui, C. Biogeo: An R package for assessing and improving data quality of occurrence record datasets. Ecography Cop. 2016, 39, 394–401. [Google Scholar] [CrossRef] [Green Version]
  65. Varela, S.; Anderson, R.P.; García-Valdés, R.; Fernández-González, F. Environmental filters reduce the effects of sampling bias and improve predictions of ecological niche models. Ecography Cop. 2014, 37, 1084–1091. [Google Scholar] [CrossRef]
  66. Wu, T.; Lu, Y.; Fang, Y.; Xin, X.; Li, L.; Li, W.; Jie, W.; Zhang, J.; Liu, Y.; Zhang, L.; et al. The Beijing Climate Center Climate System Model (BCC-CSM): The main progress from CMIP5 to CMIP6. Geosci. Model Dev. 2019, 12, 1573–1600. [Google Scholar] [CrossRef] [Green Version]
  67. Hajima, T.; Watanabe, M.; Yamamoto, A.; Tatebe, H.; Noguchi, M.; Abe, M.; Ohgaito, R.; Ito, A.; Yamazaki, D.; Okajima, H.; et al. Description of the MIROC-ES2L Earth system model and evaluation of its climate–biogeochemical processes and feedbacks. Geosci. Model Dev. Discuss. 2019, 5, 1–73. [Google Scholar] [CrossRef]
  68. Yukimoto, S.; Kawai, H.; Koshiro, T.; Oshima, N.; Yoshida, K.; Urakawa, S.; Tsujino, H.; Deushi, M.; Tanaka, T.; Hosaka, M.; et al. The meteorological research institute Earth system model version 2.0, MRI-ESM2.0: Description and basic evaluation of the physical component. J. Meteorol. Soc. Jpn. 2019, 97, 931–965. [Google Scholar] [CrossRef] [Green Version]
  69. Riahi, K.; van Vuuren, D.P.; Kriegler, E.; Edmonds, J.; O’Neill, B.C.; Fujimori, S.; Bauer, N.; Calvin, K.; Dellink, R.; Fricko, O.; et al. The Shared Socioeconomic Pathways and their energy, land use, and greenhouse gas emissions implications: An overview. Glob. Environ. Chang. 2017, 42, 153–168. [Google Scholar] [CrossRef] [Green Version]
  70. Title, P.O.; Bemmels, J.B. ENVIREM: An expanded set of bioclimatic and topographic variables increases flexibility and improves performance of ecological niche modeling. Ecography Cop. 2018, 41, 291–307. [Google Scholar] [CrossRef] [Green Version]
  71. Jarvis, A.; Reuter, H.I.; Nelson, A.; Guevara, E. Hole-Filled SRTM for the Globe; Version 4: Data Grid. Web Publication/Site; CGIAR Consortium for Spatial Information: 2008. Available online: https://research.utwente.nl/en/publications/hole-filled-srtm-for-the-globe-version-4-data-grid (accessed on 1 March 2022).
  72. Hijmans, R.J.; van Etten, J. Raster: Geographic Analysis and Modeling with Raster Data; R Package Version 3.3.13. Available online: http://CRAN.R-project.org/package=raster (accessed on 23 February 2022).
  73. Evans, J.S.; Murphy, M.A. SpatialEco; R Package Version 1.2-0; 2021. Available online: https://github.com/jeffreyevans/spatialEco (accessed on 23 February 2022).
  74. Poggio, L.; De Sousa, L.M.; Batjes, N.H.; Heuvelink, G.B.M.; Kempen, B.; Ribeiro, E.; Rossiter, D. SoilGrids 2.0: Producing soil information for the globe with quantified spatial uncertainty. Soil 2021, 7, 217–240. [Google Scholar] [CrossRef]
  75. QGIS Development Team. QGIS Geographic Information System; Open Source Geospatial Foundation Project; Open Source Geospatial Foundation: Chicago, IL, USA, 2021. [Google Scholar]
  76. Naimi, B.; Hamm, N.A.S.; Groen, T.A.; Skidmore, A.K.; Toxopeus, A.G. Where is positional uncertainty a problem for species distribution modelling? Ecography Cop. 2014, 37, 191–203. [Google Scholar] [CrossRef]
  77. Dormann, C.F.; Elith, J.; Bacher, S.; Buchmann, C.; Carl, G.; Carré, G.; Marquéz, J.R.G.; Gruber, B.; Lafourcade, B.; Leitão, P.J.; et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography Cop. 2013, 36, 27–46. [Google Scholar] [CrossRef]
  78. Araújo, M.B.; New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 2007, 22, 42–47. [Google Scholar] [CrossRef]
  79. Araújo, M.B.; Anderson, R.P.; Barbosa, A.M.; Beale, C.M.; Dormann, C.F.; Early, R.; Garcia, R.A.; Guisan, A.; Maiorano, L.; Naimi, B.; et al. Standards for distribution models in biodiversity assessments. Sci. Adv. 2019, 5, eaat4858. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  80. Vignali, S.; Barras, A.G.; Arlettaz, R.; Braunisch, V. SDMtune: An R package to tune and evaluate species distribution models. Ecol. Evol. 2020, 10, 11488–11506. [Google Scholar] [CrossRef]
  81. Carlson, C.J. embarcadero: Species distribution modelling with Bayesian additive regression trees in r. Methods Ecol. Evol. 2020, 11, 850–858. [Google Scholar] [CrossRef]
  82. Muscarella, R.; Galante, P.J.; Soley-Guardia, M.; Boria, R.A.; Kass, J.M.; Uriarte, M.; Anderson, R.P. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol. Evol. 2014, 5, 1198–1205. [Google Scholar] [CrossRef]
  83. Roberts, D.R.; Bahn, V.; Ciuti, S.; Boyce, M.S.; Elith, J.; Guillera-Arroita, G.; Hauenstein, S.; Lahoz-Monfort, J.J.; Schröder, B.; Thuiller, W.; et al. Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography Cop. 2017, 40, 913–929. [Google Scholar] [CrossRef]
  84. Carlson, C.J.; Bevins, S.N.; Schmid, B. V Plague risk in the western United States over seven decades of environmental change. Glob. Chang. Biol. 2022, 28, 753–769. [Google Scholar] [CrossRef]
  85. Phillips, S.J.; Anderson, R.P.; Dudík, M.; Schapire, R.E.; Blair, M.E. Opening the black box: An open-source release of Maxent. Ecography Cop. 2017, 40, 887–893. [Google Scholar] [CrossRef]
  86. Velasco, J.A.; González-Salazar, C. Akaike information criterion should not be a “test” of geographical prediction accuracy in ecological niche modelling. Ecol. Inform. 2019, 51, 25–32. [Google Scholar] [CrossRef]
  87. Barbet-Massin, M.; Jiguet, F.; Albert, C.H.; Thuiller, W. Selecting pseudo-absences for species distribution models: How, where and how many? Methods Ecol. Evol. 2012, 3, 327–338. [Google Scholar] [CrossRef]
  88. Liu, C.; Newell, G.; White, M. The effect of sample size on the accuracy of species distribution models: Considering both presences and pseudo-absences or background sites. Ecography Cop. 2019, 42, 535–548. [Google Scholar] [CrossRef] [Green Version]
  89. Guevara, L.; Gerstner, B.E.; Kass, J.M.; Anderson, R.P. Toward ecologically realistic predictions of species distributions: A cross-time example from tropical montane cloud forests. Glob. Chang. Biol. 2018, 24, 1511–1522. [Google Scholar] [CrossRef]
  90. Konowalik, K.; Nosol, A. Evaluation metrics and validation of presence-only species distribution models based on distributional maps with varying coverage. Sci. Rep. 2021, 11, 1–15. [Google Scholar] [CrossRef]
  91. Somodi, I.; Lepesi, N.; Botta-Dukát, Z. Prevalence dependence in model goodness measures with special emphasis on true skill statistics. Ecol. Evol. 2017, 7, 863–872. [Google Scholar] [CrossRef] [Green Version]
  92. Leroy, B.; Delsol, R.; Hugueny, B.; Meynard, C.N.; Barhoumi, C.; Barbet-Massin, M.; Bellard, C. Without quality presence-absence data, discrimination metrics such as TSS can be misleading measures of model performance. J. Biogeogr. 2018, 45, 1994–2002. [Google Scholar] [CrossRef]
  93. Allouche, O.; Tsoar, A.; Kadmon, R. Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 2006, 43, 1223–1232. [Google Scholar] [CrossRef]
  94. Fielding, A.H.; Bell, J.F. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 1997, 24, 38–49. [Google Scholar] [CrossRef]
  95. Hirzel, A.H.; Le Lay, G.; Helfer, V.; Randin, C.; Guisan, A. Evaluating the ability of habitat suitability models to predict species presences. Ecol. Modell. 2006, 199, 142–152. [Google Scholar] [CrossRef]
  96. Liu, C.; Newell, G.; White, M. On the selection of thresholds for predicting species occurrence with presence-only data. Ecol. Evol. 2016, 6, 337–348. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  97. Liu, C.; White, M.; Newell, G. Measuring and comparing the accuracy of species distribution models with presence-absence data. Ecography Cop. 2011, 34, 232–243. [Google Scholar] [CrossRef]
  98. Sofaer, H.R.; Jarnevich, C.S.; Pearse, I.S.; Smyth, R.L.; Auer, S.; Cook, G.L.; Edwards, T.C.; Guala, G.F.; Howard, T.G.; Morisette, J.T.; et al. Development and Delivery of Species Distribution Models to Inform Decision-Making. Bioscience 2019, 69, 544–557. [Google Scholar] [CrossRef]
  99. Schwarz, J.; Heider, D. GUESS: Projecting machine learning scores to well-calibrated probability estimates for clinical decision-making. Bioinformatics 2019, 35, 2458–2465. [Google Scholar] [CrossRef] [PubMed]
  100. Márcia Barbosa, A.; Real, R.; Muñoz, A.R.; Brown, J.A. New measures for assessing model equilibrium and prediction mismatch in species distribution models. Divers. Distrib. 2013, 19, 1333–1338. [Google Scholar] [CrossRef]
  101. Signorell, A.; Aho, K.; Anderegg, N.; Aragon, T.; Arppe, A.; Baddeley, A.; Bolker, B.; Caeiro, F.; Champely, S.; Chessel, D. DescTools: Tools for Descriptive Statistics; R Package Version 0.99-40; CRAN, R Core Team: Cary, CA, USA, 2021. [Google Scholar]
  102. Broennimann, O.; Di Cola, V.; Guisan, A. Ecospat: Spatial Ecology Miscellaneous Methods; R Package Version 3.2; CRAN, R Core Team: Cary, CA, USA, 2021. [Google Scholar]
  103. Hammer, B.; Frasco, M. Metrics: Evaluation Metrics for Machine Learning; R Package Version 0.1.4; CRAN, R Core Team: Cary, CA, USA, 2018. [Google Scholar]
  104. Yan, Y. MLmetrics: Machine Learning Evaluation Metrics, R package version 1.1.1; CRAN, R Core Team: Cary, CA, USA, 2016. [Google Scholar]
  105. Zhu, G.; Fan, J.; Peterson, A.T. Cautions in weighting individual ecological niche models in ensemble forecasting. Ecol. Modell. 2021, 448, 109502. [Google Scholar] [CrossRef]
  106. Liu, C.; Berry, P.M.; Dawson, T.P.; Pearson, R.G. Selecting thresholds of occurrence in the prediction of species distributions. Ecography Cop. 2005, 28, 385–393. [Google Scholar] [CrossRef]
  107. Liu, C.; White, M.; Newell, G. Selecting thresholds for the prediction of species occurrence with presence-only data. J. Biogeogr. 2013, 40, 778–789. [Google Scholar] [CrossRef]
  108. Elith, J.; Kearney, M.; Phillips, S. The art of modelling range-shifting species. Methods Ecol. Evol. 2010, 1, 330–342. [Google Scholar] [CrossRef]
  109. Thuiller, W.; Lafourcade, B.; Engler, R.; Araújo, M.B. BIOMOD—A platform for ensemble forecasting of species distributions. Ecography Cop. 2009, 32, 369–373. [Google Scholar] [CrossRef]
  110. Arditti, J.; Ghani, A.K.A. Numerical and physical properties of orchid seeds and their biological implications. New Phytol. 2000, 145, 367–421. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  111. Pebesma, E. Simple features for R: Standardized support for spatial vector data. R J. 2018, 10, 439–446. [Google Scholar] [CrossRef] [Green Version]
  112. Wickham, H. ggplot2: Elegant Graphics for Data Analysis; Springer: New York, NY, USA, 2016; ISBN 978-3-319-24277-4. [Google Scholar]
  113. Carey, P.D. Changes in the distribution and abundance of Himantoglossum hircinum (L.) Sprengel (Orchidaceae) over the last 100 years. Flora 1999, 364, 353–364. [Google Scholar]
  114. Kull, T.; Hutchings, M.J. A comparative analysis of decline in the distribution ranges of orchid species in Estonia and the United Kingdom. Biol. Conserv. 2006, 129, 31–39. [Google Scholar] [CrossRef]
  115. Vogt-Schilb, H.; Munoz, F.; Richard, F.; Schatz, B. Recent declines and range changes of orchids in Western Europe (France, Belgium and Luxembourg). Biol. Conserv. 2015, 190, 133–141. [Google Scholar] [CrossRef]
  116. Reina-Rodríguez, G.A.; Rubiano Mejía, J.E.; Castro Llanos, F.A.; Soriano, I. Orchids distribution and bioclimatic niches as a strategy to climate change in areas of tropical dry forest in Colombia. Lankesteriana 2017, 17, 17–47. [Google Scholar] [CrossRef] [Green Version]
  117. Djordjević, V.; Tsiftsis, S. The Role of Ecological Factors in Distribution and Abundance of Terrestrial Orchids. In Orchids Phytochemistry, Biology and Horticulture; Mérillon, J.-M., Kodja, H., Eds.; Springer: Cham, Switzerland, 2020; pp. 1–71. ISBN 9783030112578. [Google Scholar]
  118. Zhang, S.; Chen, W.; Huang, J.; Bi, Y.; Yang, X. Orchid Species Richness along Elevational and Environmental Gradients in Yunnan, China. PLoS ONE 2015, 10, e0142621. [Google Scholar] [CrossRef]
  119. Dorland, E.; Willems, J.H. Light climate and plant performance of Ophrys insectifera; a four-year field experiment in The Netherlands (1998–2001). In Trends and Fluctuations and Underlying Mechanisms in Terrestrial Orchid Populations; Kindlmann, P., Willems, J.H., Whigham, D.F., Eds.; Backhuys Publishers: Leiden, The Netherlands, 2002; pp. 225–238. ISBN 9057821230. [Google Scholar]
  120. Preston, C.D.; Pearman, D.A.; Dines, T.D. New Atlas of the British and Irish Flora; Oxford University Press: Oxford, UK, 2002. [Google Scholar]
  121. Jacquemyn, H.; Brys, R.; Hermy, M.; Willems, J.H. Does nectar reward affect rarity and extinction probabilities of orchid species? An assessment using historical records from Belgium and the Netherlands. Biol. Conserv. 2005, 121, 257–263. [Google Scholar] [CrossRef]
  122. Damgaard, C.; Moeslund, J.E.; Wind, P. Changes in the abundance of Danish orchids over the past 30 years. Diversity 2020, 12, 244. [Google Scholar] [CrossRef]
  123. Wraith, J.; Pickering, C. A continental scale analysis of threats to orchids. Biol. Conserv. 2019, 234, 7–17. [Google Scholar] [CrossRef]
  124. Triponez, Y.; Arrigo, N.; Pellissier, L.; Schatz, B.; Alvarez, N. Morphological, ecological and genetic aspects associated with endemism in the Fly Orchid group. Mol. Ecol. 2013, 22, 1431–1446. [Google Scholar] [CrossRef]
  125. Waterman, R.J.; Bidartondo, M.I. Deception above, deception below: Linking pollination and mycorrhizal biology of orchids. J. Exp. Bot. 2008, 59, 1085–1096. [Google Scholar] [CrossRef]
  126. Shefferson, R.P.; Jacquemyn, H.; Kull, T.; Hutchings, M.J. The demography of terrestrial orchids: Life history, population dynamics and conservation. Bot. J. Linn. Soc. 2020, 192, 315–332. [Google Scholar] [CrossRef]
  127. Schweiger, J.M.-I.; Bidartondo, M.I.; Gebauer, G. Stable isotope signatures of underground seedlings reveal the organic matter gained by adult orchids from mycorrhizal fungi. Funct. Ecol. 2018, 32, 870–881. [Google Scholar] [CrossRef]
  128. Breitkopf, H.; Onstein, R.E.; Cafasso, D.; Schlüter, P.M.; Cozzolino, S. Multiple shifts to different pollinators fuelled rapid diversification in sexually deceptive Ophrys orchids. New Phytol. 2015, 207, 377–389. [Google Scholar] [CrossRef] [PubMed]
  129. Charitonidou, M.; Halley, J.M. What goes up must come down—Why high fecundity orchids challenge conservation beliefs. Biol. Conserv. 2020, 252, 108835. [Google Scholar] [CrossRef]
Figure 1. (a) Ophrys insectifera distribution map, compiled with data downloaded from IUCN [53], (b) map including the total occurrence data for O. insectifera retrieved from GBIF and other sources, and the study area defined as a buffered α-hull extent of occurrence (EOO). Both maps are designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System. Photo of the species taken in June 2019 at Mt. Tzoumérka (Epirus, NW Greece).
Figure 1. (a) Ophrys insectifera distribution map, compiled with data downloaded from IUCN [53], (b) map including the total occurrence data for O. insectifera retrieved from GBIF and other sources, and the study area defined as a buffered α-hull extent of occurrence (EOO). Both maps are designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System. Photo of the species taken in June 2019 at Mt. Tzoumérka (Epirus, NW Greece).
Biology 11 00497 g001
Figure 2. Current habitat suitability map for Ophrys insectifera, for the ensemble model using the environmental thinning procedure. Map is designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Figure 2. Current habitat suitability map for Ophrys insectifera, for the ensemble model using the environmental thinning procedure. Map is designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Biology 11 00497 g002
Figure 3. Current and future potential distribution maps for Ophrys insectifera EnvThin ensemble model. Left-hand panel: red-brown coloring indicates the cells the species currently potentially occupies. Each map is showing the transition from the present time-period to each respective GCM and SSP combination. Right-hand panel: future potential distribution maps of 2070 for the combinations of three GCMs (BCC, MIROC, and MRI) and two SSPs (SSP1-2.6 and SSP5-8.5) as ‘best’ and ‘worst’ case scenario. Grid cells with red coloring indicate the areas where the species is currently present but will not be in the future. Grey coloring represents cells where the species currently occupies and will continue to occupy in the future. White stands for the cells where the species is not currently present will not be in the future, while blue grid cells indicate the areas where the species is not currently present but will occupy in the future. All maps are designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Figure 3. Current and future potential distribution maps for Ophrys insectifera EnvThin ensemble model. Left-hand panel: red-brown coloring indicates the cells the species currently potentially occupies. Each map is showing the transition from the present time-period to each respective GCM and SSP combination. Right-hand panel: future potential distribution maps of 2070 for the combinations of three GCMs (BCC, MIROC, and MRI) and two SSPs (SSP1-2.6 and SSP5-8.5) as ‘best’ and ‘worst’ case scenario. Grid cells with red coloring indicate the areas where the species is currently present but will not be in the future. Grey coloring represents cells where the species currently occupies and will continue to occupy in the future. White stands for the cells where the species is not currently present will not be in the future, while blue grid cells indicate the areas where the species is not currently present but will occupy in the future. All maps are designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Biology 11 00497 g003
Figure 4. Median centroids for current and future projected distributions of Ophrys insectifera, for the EnvThin ensemble in 2070 time slice. Black point represents the current distribution’s median centroid, the blue points stand for the future best-case scenario (SSP1-2.6), and the magenta points for the future worst-case scenario (SSP5-8.5). GCMs differentiate by shape: BCC—circle, MIROC—diamond, MRI—square. Map is designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Figure 4. Median centroids for current and future projected distributions of Ophrys insectifera, for the EnvThin ensemble in 2070 time slice. Black point represents the current distribution’s median centroid, the blue points stand for the future best-case scenario (SSP1-2.6), and the magenta points for the future worst-case scenario (SSP5-8.5). GCMs differentiate by shape: BCC—circle, MIROC—diamond, MRI—square. Map is designed in QGIS v.3.18 ‘Zürich’, using ETRS89—Lambert Conformal Conic Coordinate Reference System.
Biology 11 00497 g004
Figure 5. Temperature seasonality variation for two selected locations of Ophrys insectifera distribution over the time period 1986–2100. UK (in blue) represents an area of range gain, while Greece (in red) an area of range loss. Two shared socioeconomic pathways (SSPs) are presented, as best-case and worst-case scenarios (SSP1-2.6 and SSP5-8.5, respectively). In all panels, each line represents a global circulation model (GCM): BCC-CSM2-MR (cyan line); MIROC-ES2L (yellow line); MRI-ESM2-0 (coral red line).
Figure 5. Temperature seasonality variation for two selected locations of Ophrys insectifera distribution over the time period 1986–2100. UK (in blue) represents an area of range gain, while Greece (in red) an area of range loss. Two shared socioeconomic pathways (SSPs) are presented, as best-case and worst-case scenarios (SSP1-2.6 and SSP5-8.5, respectively). In all panels, each line represents a global circulation model (GCM): BCC-CSM2-MR (cyan line); MIROC-ES2L (yellow line); MRI-ESM2-0 (coral red line).
Biology 11 00497 g005
Table 1. Percentages of area loss, gain, overall area change and current occurrences lost for future projections of Ophrys insectifera. For the future, 2 time slices, 2070 (2061–2080) and 2090 (2081–2100), are presented for all selected global circulation models (BC: BCC-CSM2-MR, MI: MIROC-ES2L, MR: MRI-ESM2-0) and two shared socioeconomic pathways (SSP1: best-case scenario, SSP5: worst-case scenario). The presented values are for the ensemble model with the environmental thinning procedure.
Table 1. Percentages of area loss, gain, overall area change and current occurrences lost for future projections of Ophrys insectifera. For the future, 2 time slices, 2070 (2061–2080) and 2090 (2081–2100), are presented for all selected global circulation models (BC: BCC-CSM2-MR, MI: MIROC-ES2L, MR: MRI-ESM2-0) and two shared socioeconomic pathways (SSP1: best-case scenario, SSP5: worst-case scenario). The presented values are for the ensemble model with the environmental thinning procedure.
Time SliceTransitionGCMArea Loss (%)Area Gain (%)Overall Change (%)Current Occurrences Lost (%)
2070Present to SSP1-2.6BC20.5817.45−3.148.27
MI16.4021.124.7210.54
MR38.2914.36−23.9322.03
Mean25.0917.64−7.4513.60
Present to SSP5-8.5BC51.4223.89−27.5438.46
MI20.0329.369.3313.37
MR44.0015.44−28.5531.13
Mean38.4822.90−15.5927.65
2090Present to SSP1-2.6BC20.0715.03−5.048.99
MI16.1322.916.7810.27
MR31.6413.25−18.3916.98
Mean22.6117.06−5.5512.08
Present to SSP5-8.5BC77.0325.65−51.3873.20
MI40.5123.41−17.0930.91
MR29.0615.00−14.0616.26
Mean48.8721.35−27.5140.12
Table 2. Mean altitude for present and future projections of Ophrys insectifera distribution. For the future, two time slices, 2070 (2061–2080) and 2090 (2081–2100), are presented for all selected global circulation models (GCMs) and two shared socioeconomic pathways (SSPs), representing ‘best’- and ‘worst’-case scenarios respectively. BC: BCC-CSM2-MR, MI: MIROC-ES2L, MR: MRI-ESM2-0. SSP1-2.6: ‘best’-case scenario, SSP5-8.5: ‘worst’-case scenario. The presented values are for the EnvThin ensemble model.
Table 2. Mean altitude for present and future projections of Ophrys insectifera distribution. For the future, two time slices, 2070 (2061–2080) and 2090 (2081–2100), are presented for all selected global circulation models (GCMs) and two shared socioeconomic pathways (SSPs), representing ‘best’- and ‘worst’-case scenarios respectively. BC: BCC-CSM2-MR, MI: MIROC-ES2L, MR: MRI-ESM2-0. SSP1-2.6: ‘best’-case scenario, SSP5-8.5: ‘worst’-case scenario. The presented values are for the EnvThin ensemble model.
Time SliceSSPPeriod/GCMMean Altitude (m)
Present Current542.6
2070SSP1-2.6BC484.7
MI520.3
MR467.5
SSP5-8.5BC428.6
MI447.8
MR447.0
2090SSP1-2.6BC504.6
MI505.0
MR498.3
SSP5-8.5BC278.4
MI442.9
MR492.7
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Charitonidou, M.; Kougioumoutzis, K.; Karypidou, M.C.; Halley, J.M. ‘Fly to a Safer North’: Distributional Shifts of the Orchid Ophrys insectifera L. Due to Climate Change. Biology 2022, 11, 497. https://doi.org/10.3390/biology11040497

AMA Style

Charitonidou M, Kougioumoutzis K, Karypidou MC, Halley JM. ‘Fly to a Safer North’: Distributional Shifts of the Orchid Ophrys insectifera L. Due to Climate Change. Biology. 2022; 11(4):497. https://doi.org/10.3390/biology11040497

Chicago/Turabian Style

Charitonidou, Martha, Konstantinos Kougioumoutzis, Maria Chara Karypidou, and John M. Halley. 2022. "‘Fly to a Safer North’: Distributional Shifts of the Orchid Ophrys insectifera L. Due to Climate Change" Biology 11, no. 4: 497. https://doi.org/10.3390/biology11040497

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop