Modelling sexually deceptive orchid species distributions under future climates: the importance of plant–pollinator interactions

Biotic interactions play an important role in species distribution models, whose ignorance may cause an overestimation of species' potential distributions. Species of the family Orchidaceae are almost totally dependent on mycorrhizal symbionts and pollinators, with sexually deceptive orchids being often highly specialized, and thus the interactions with their pollinators are expected to strongly affect distribution predictions. We used Maxent algorithm to explore the extent of current and future habitat suitability for two Greek endemic sexually deceptive orchids (Ophrys argolica and Ophrys delphinensis) in relation to the potential distribution of their unique pollinator (Anthophora plagiata). Twelve climate change scenarios were used to predict future distributions. Results indicated that the most important factors determining potential distribution were precipitation seasonality for O. argolica and geological substrate for O. delphinensis. The current potential distribution of the two orchids was almost of the same extent but spatially different, without accounting for their interaction with A. plagiata. When the interaction was included in the models, their potentially suitable area decreased for both species. Under future climatic conditions, the effects of the orchid-pollinator interaction were more intense. Specifically, O. argolica was restricted in specific areas of southern Greece, whereas O. delphinensis was expected to become extinct. Our findings highlighted the significant role of plant–pollinator interactions in species distribution models. Failing to study such interactions might expose plant species to serious conservation issues.

range. However, it is well known that species distribution and abundance is not only determined by abiotic factors but also by other ones [11][12][13] .
The potentially suitable area for species whose existence is exclusively dependent on other organisms (e.g., plant-pollinator interaction) is expected to be overestimated when SDMs are used without taking into consideration these interactions. As a result, biotic interactions are increasingly used directly in the SDMs or in combination with their outcomes [14][15][16][17] .
The family Orchidaceae is one of the richest and most fascinating among the families of the plant kingdom, consisting of 736 genera 18 and approximately 28,500 species 19 . Their distribution is not only driven by climate or other site characteristics (e.g., geology) but also by their interactions with specific biotic factors (e.g., mycorrhizal fungi, pollinators) [20][21][22] . Many species of this family are disappearing worldwide 23 , mostly due to habitat loss and factors such as climate change and the resulting shifts in species distributions 24,25 . Furthermore, orchids could also be negatively affected by biotic factors such as the concurrent reduction of their pollinators. These negative effects are expected to be exacerbated by climate change during the twenty-first century.
A considerable number of orchids are pollinated by specific pollinators 26 . According to Tremblay 26 , 62% of the studied orchids had a single pollinator and this percentage varied in the different subfamilies of the Orchidaceae. The most characteristic orchid pollination system where species have a high plant pollinator specificity is sexual deception 27,28 . This pollination system characterizes the European genus Ophrys and the species Orchis galilaea (Bornm. and M. Schulze) Schltr., two African Disa species (Jersáková et al. 29 , and references therein), more than 150 orchid species-classified in at least nine genera-in Australia 28 , as well as six genera of Central and South America 29,30 . Sexual deception by orchids imposes a high level of specialization, as insect pheromones are usually species-specific. Moreover, a large proportion of Ophrys species pollinators are solitary bees belonging to the oligolectic group that specializes in a particular pollen species 31 . Specialization varies from species attracted by several pollinator taxa to species pollinated by only one pollinator species. Orchids' specialized pollinators often require specific conditions (e.g., specific nesting sites), further denoting the fragile plant-insect interactions 31,32 . Due to the high pollinator specificity, it is expected that these orchids cannot form reproductive populations in areas where their pollinators are not present. As a result, although they can spread through seed dispersal and colonize areas far away from their stabilized populations (some terrestrial orchids are known to have colonized areas c. 250 km away from their nearest known populations, whereas seed dispersal over distances of 5-10 km seems to be very common 33 ), they will only be able to survive for a number of years and finally they will disappear.
SDMs have been widely used to predict the potential distribution of orchid species at the global, European or regional scale (e.g. 5,16,[34][35][36][37][38][39][40]. Although most of these studies used orchids per se and were focused on species ecology and/or conservation, only one of them explored plant-pollinator co-occurrence and how the predictions of SDMs were affected by these interactions 16 . Their study focused on the current distribution of several orchids, but the results of SDMs were not projected to future climates. Orchids and their pollinators might react in different ways to global warming. This could have serious conservation consequences for the dependent species (orchids) due to the possibility of future divergence between the potential distributions of the plants and their pollinators. We did not find any study on the effects of plant-pollinator interactions on future orchid distribution under climate change scenarios in the literature.
Based on this background, the aim of the present study was to explore how orchids' potential distribution was affected by plant-pollinator interactions under current and future climatic conditions. Thus, the major objectives were (a) to compare orchids' projections based on a set of environmental predictors and (b) to compare orchids' projections with the projections of an insect pollinator. Specifically, we hypothesized that the potential distribution (under current and future climatic conditions) of two Ophrys taxa (O. argolica and O. delphinensis), would be affected by their common and unique pollinator-Anthophora plagiata 41 -and consequently the predicted distribution of the pollinator would restrict the potential distribution of these two Ophrys taxa.

Results
The results of the Jackknife test about variable importance are presented in Table 1 Ophrys argolica, per se, had its optimum in almost the whole area of the Peloponnese except for the northcentral parts (Fig. 1a). However, when its current potential distribution was combined with the modelled distribution of A. plagiata, its potentially suitable area was reduced, mostly restricted in areas of lower altitude and close to the sea (Fig. 1b), whereas the percentage of the area of Greece that was potentially suitable for this species fell from 11 to c. 7% (Fig. 2a).
Under most of the future climatic scenarios, the potential distribution of O. argolica was greatly increased compared to the current potential distribution, reaching c. 20% of the total area of Greece (Fig. 2a, Supplementary Information 3), and only areas towards the highest mountain tops of the Peloponnese seemed unsuitable.    (Fig. 5a). Moreover, although the model identified potentially suitable areas in the northern part of Greece, these areas could not be colonized by the specific orchid for biogeographical reasons. However, when taking into consideration the potential distribution of its pollinator (A. plagiata) the areas of high suitability values were greatly reduced, mostly restricted to the northern part of the Peloponnese and the southern part of Sterea Hellas (Fig. 5b).
Surprisingly, under almost all the future climatic scenarios, the most suitable sites for the distribution of O. delphinensis were restricted in the central-east part of mainland Greece, as well as in the north-east part of Peloponnese, whereas in these areas the habitat suitability values fell below the P10 threshold ( Supplementary  Information 4). However, even this restricted distribution was not congruent at all with the future potential distribution of A. plagiata ( Supplementary Information 2), demonstrating that the species will not be able to survive (Figs. 2b, 6; Supplementary Information 5).
Anthophora plagiata was overlapping to a greater degree with O. argolica (I statistic = 0.453) compared to O. delphinensis (I statistic = 0.186), under current climate conditions ( Table 2). It is worth mentioned that in two out of three climate change models (CCSM4 and HadGEM2-ES), niche overlap between these two species (O. argolica and A. plagiata) was increasing from 2050 to 2070 under the medium greenhouse gas emission scenario

Discussion
Conservation actions often require knowledge of species' realized niche, which is defined as a subset of species' fundamental niche under specific restrictions (e.g., biotic interactions) 4,42 . Contrary to the realized niche, fundamental niche is independent of any biotic interaction and is referred to the potential distribution of a species, which is determined by the SDMs on the basis of the environmental conditions of the presence locations. In our study, the distribution of the two orchids seemed to be determined by different factors. In the case of O. argolica, its distribution was mostly determined by the effect of precipitation, whereas geological substrate was not important. The opposite trend had been found in the case of O. delphinensis. This could be due to the fact that the studied orchids have different preferences regarding the geological substrate of the sites where they have been recorded so far. Although O. delphinensis is mostly found in limestones, O. argolica is more generalist and can be found in a great variety of substrates 43 . Moreover, contrary to other orchids whose existence strongly depends on vegetation (e.g. 5 ), both orchids can be found in several types of habitats (e.g., grasslands, scrubs, forests), even in anthropogenic habitats, such as olive groves which are common in the southern parts of Greece 43 .
Although proper selection of the environmental variables used in the SDMs is of crucial importance, species interactions should not be ignored. Our results showed that biotic interactions should be considered as a key component in modelling species' distributions, because potentially suitable areas could be overestimated when they are not included in the models. Orchids depend on pollinators and mycorrhizal symbionts for their survival, and some of them present a high degree of specialization 20,22 . However, orchid mycorrhizal fungi are often present at sites where orchids are lacking, suggesting that fungi are not particularly limiting orchid distribution 22,44 . On the contrary, most terrestrial orchids can survive pollinator deficiency for only 3-10 years 32 . However, it is expected that climate change will severely affect pollinators in the future. Unfavourable future climatic conditions can either cause a great reduction in insects' populations or can change the time of their emergence, or finally, can cause a range shift towards areas of favourable climate 45,46 . In turn, these changes will cause a temporal or/ and spatial mismatch between plants and their pollinators and will affect plant populations and their viability 47 . Therefore, plant-pollinator interactions could heavily affect the outcome of SDMs, especially for orchid species with high pollinator specificity.
The analysis of the current potential distribution of the two orchid species, both with and without plant-pollinator interactions, demonstrated the significance of such interactions in the outcomes of SDMs. Specifically, although the potential distribution of both orchids was almost of the same extent (10.74% for O. argolica and 8.79% of the area of Greece for O. delphinensis), after calculating habitat suitability values on the basis of their interaction with A. plagiata, the potentially suitable area decreased by c. 52% and c. 73%, respectively. These results are consistent with the findings of Duffy and Johnson 16 who found that pollinators constrain the distribution of the orchids which are characterized by a highly specialized pollination system. www.nature.com/scientificreports/ Our study also explored the consequences of global warming to the future potential distribution of the two studied orchids. There are several ways in which organisms may respond to climate change. Species either can adapt to the new ecological conditions, by shifting their niche, or they can move to other areas (e.g., northwards). However, the responses of species that have strong and specialized interactions with others are always very complex. Robbirt et al. 48 found that Ophrys sphegodes and its pollinator, male bees of the species Andrena nigroaenea, showed phenological responses to climate warming, by reducing the frequency of pseudocopulation and, as a result, pollination success and seed production. Another eventuality is that species that interact with others can migrate and colonize different areas under the effects of climate change. In doing so, their area of co-occurrence could decrease, and in some cases the hitherto interacting species could be even spatially separated 17,49 .
In our study, the two orchids were characterized by entirely different trends under future climate change scenarios. O. argolica was expected to retain its distribution range under the CCSM4 climate change model when its interaction with A. plagiata was considered, whereas on the contrary, it was expected to narrow its range under all other scenarios. The fact that habitat suitability values of the co-occurrence model O. argolica × A. plagiata were high for a large number of grid cells under most scenarios, demonstrated that both species can co-occur successfully for a relatively extensive area, as both are characterized by high suitability values.
The worst case scenario is that species might be unable to adapt to the changing environment 50 . In our case, O. delphinensis was expected to become extinct in the future, as it was predicted that O. delphinensis and A. plagiata will not occur in the same areas. Although this is a pessimistic scenario, it was the simple result of an algorithm, Maxent model, which is not able to consider limitations such as species evolution and niche shift 51,52 . Ophrys www.nature.com/scientificreports/ delphinensis, as all Ophrys species that are highly specialized, attracts its unique pollinator (A. plagiata) using sexual deception. The orchid emits volatile chemical compounds similar to the pheromones of female A. plagiata, so attracting its pollinators who attempt copulation (pseudocopulation) with the orchid lip 31,53 . Sexual deception is an evolutionary trait of Ophrys taxa and Schlüter et al. 54 found that small changes in the genes involved in scent production could cause the attraction of a different pollinator. This allow for the hypothesis that O. delphinensis could be able to attract another insect-pollinator as the result of mutations or genetic drift. The fact that O. delphinensis is considered as a stabilized hybrid between O. argolica and a taxon of the Ophrys oestrifera group, two taxa that are pollinated by different bees (reviewed in Ref. 43 ), offers further support to this hypothesis. Our findings suggested that the effect of the interaction of the two orchids with their unique pollinator is very important under both current and future climate conditions. Although the study dealt only with plant-pollinator interactions, future research should also focus on the potential effect of their mycorrhizal symbionts. By doing this, we will be able to explore whether orchids can shift their distribution under climate change using a more complete spectrum of environmental and biotic constraints. Given that many orchids, either in Greece or in Europe 55,56 , have declined dramatically during the twentieth century and continue to decline, better insights into the processes affecting the distribution and abundance of orchids will be of high value for conservation management and the restoration of orchid populations. This study provides information on the distribution of Ophrys argolica and O. deplhinensis, whereas can serve as a basis for other studies dealing with the potential distribution of orchid species and their pollinators. Moreover, although the orchid datasets used in the analyses are good from a quantitative and qualitative perspective, a systematic survey for the presence of A. plagiata in other parts of Greece could provide valuable ecological information and improve the results of the SDMs and consequently the effective conservation of the studied orchids.  43,57 . Based on their known distribution, the study area was set between 36° 05′-41° 37′ N and 19° 30′-24° 40′ E. Species of the European genus Ophrys usually present a high degree of pollinator specificity as many of them are being pollinated by a single pollinator 27 . Both, O. argolica and O. delphinensis are sexually deceptive orchids and are exclusively pollinated by the same pollinator, Anthophora plagiata (Illiger, 1806) 41,43,58,59 , an insect which is distributed in several European countries but is rather rare in the southern part of Europe 60 .

Material and methods
Species distribution models. The extent of current and future habitat suitability for O. argolica and O. delphinensis in Greece in relation to the potential distribution of their unique pollinator (A. plagiata) was explored using Maxent software (version 3.4.1) 61,62 . Maxent is considered as an appropriate technique for modelling species distributions, even in the case of very small sample size 7,63 , and thus is the most widely used method. The model estimates species distribution based on presence-only data, using the maximum entropy principle, and computes a probability distribution based on environmental variables spread over the entire study area. However, the probability value of each grid cell corresponds to the habitat suitability value of the specific grid cell.
In total, three datasets of species occurrences were created, one for each plant and pollinator species. The distribution data for the two orchid species used in the respective models, were based on the database that was built for the purposes of the Orchid Flora of Greece project 64 . The third dataset corresponded to the occurrences of the pollinator of the two studied orchids. However, there aren't available detailed data about the distribution of A. plagiata in Greece (either in published sources or in GBIF-www.gbif.org). To overcome this problem, we used a combined database of orchid occurrences, consisted of the total number of records of the two studied orchids, as well as all the known records of the taxon Ophrys olympiotissa Paulus, also pollinated by A. plagiata 57 , as a surrogate. It is worth mentioned that in the vast majority of these sites, pollinated flowers of the three orchids have been detected, clearly indicating that the pollinator (A. plagiata) is present. So, in this way by using this approach we assumed that the three orchids and their pollinator were present in the same sites. www.nature.com/scientificreports/ Model settings and selection. Together with the species distribution data, a second input in Maxent software corresponded to the environmental variables that were used to predict the potential distribution of the studied species. Initially, 20 environmental variables were selected as predictors in species distribution modelling. Nineteen of them were bioclimatic variables and one corresponded to the geological substrate of Greece. The bioclimatic variables were obtained from the WorldClim database 65,66 in a 30-s resolution (approximately 1 km 2 ), whereas the environmental layer corresponding to the geological substrate was at a scale of 1:500,000 67 . Although the map of the geological substrate was in vector format, the layer was converted in raster format at the same resolution and extent with the layers of the bioclimatic variables. To account for multicollinearity between the 19 bioclimatic variables and avoid overfitting, Pearson correlation coefficients were calculated for all pairwise interactions, for each species. To eliminate highly correlated variables, only one (i.e. the one with the higher percent contribution and training gain) was selected among any pair of those with a correlation coefficient r >|0.75| 68,69 . Specifically, in modelling the potential distribution of the orchid species, the non highly inter-correlated bioclimatic variables and the geological substrate were used. On the contrary, in modelling the pollinator of the Ophrys spp. only the selected bioclimatic variables were used. Moreover, having in mind that A. plagiata has a rather wide distribution in Europe 60 and we used distribution data only from the southern part of mainland Greece, a bias file was used, with which the background samples exploited by the algorithm were restricted to areas that fall within a distance (20 km) from the presence locations 70 . The idea behind that data handling was that biases in background samples will account for similar biases in presence locations 70 . Five (5) models for each species were run with Maxent using the auto-features mode and the default settings, as suggested by Phillips and Dudík 61 . Specifically, we used 10,000 background points and we increased to 5,000 the number of maximum iterations. This was done to provide adequate time to the model for convergence. Moreover, we removed duplicate presence records to reduce the effect of spatial auto-correlation. Maxent models were run using cross-validation as a form of replication. Model performance was assessed using the Akaike information criterion (AICc) because it greatly outperforms BIC and AUC based methods, especially when sample size is small 71 . climate projection models. To explore the future potential distribution of the two studied orchid species were used. These three climate change models have shown good performance in simulating future climate, and as a result they have been widely used in the last decade in studies of range shifts 17,72,73 . The four available pathways differ in the greenhouse gases that are to be emitted in the future. Specifically, the RCP2.6 is a low emission scenario, RCP4.5 and RCP6.0 are medium greenhouse gas emission scenarios, whereas the RCP8.5 is a maximum emission scenario. Here we used the maximum greenhouse gas emission scenario (RCP8.5) and one out of the two medium greenhouse gas emission scenarios (RCP 4.5).
Co-occurrence models. Except of the direct output of the Maxent models, the main aim of the present study was to identify the areas where the two studied orchids can potentially occur in relation to their pollinator. In probability theory, the probability of co-occurrence of two species (A and B) is given by the product of P(A) × P(B) 74 . However, this is accurate and correct only if species occur independently and are not interacting. In general, species distributions derived by species distribution models is assumed to be independent, and these distributions are not influenced from any ecological interaction among species 75 .
To calculate the probability of suitable environmental conditions for both species (species co-occurrence) we used only the habitat suitability values that were above the 10th percentile training presence threshold (P10) for each species separately, whereas values lower than the P10 were set to zero. Among the several thresholds that can be calculated by Maxent software, we used the 10th percentile training presence threshold, as it provides a better ecologically significant result in comparison with others 61 .
Based on the above-mentioned assumption, the probability (habitat suitability value) of orchid-pollinator co-occurrence (species co-occurrence model: plant × pollinator) can be calculated using the formula: where, P(A i ') is the habitat suitability value for each pair of co-occurring species (Ophrys argolica-Anthophora plagiata and Ophrys delphinensis-Anthophora plagiata, respectively), P(A i ) is the habitat suitability value of species A for the grid cell i www.nature.com/scientificreports/ under the future climatic scenarios. This statistic uses suitability scores and has been widely used in the past (e.g. 77,78 ). It measures niche overlap between two species distribution models and its values range from 0 (no overlap between the two distributions) to 1 (identical distributions). Prior the calculation of the "I statistic", habitat suitability values below the P10 threshold for each species were set to zero. All analyses were performed in R version 3.5.2 (R Foundation for Statistical Computing) using the packages 'raster' , 'rgdal' , ' dismo' , 'ENMeval' and 'SDMTools' , whereas the distribution maps were generated in ArcMap 10.1 79 .

Data availability
Species records that were used in developing niche models are available from the corresponding author upon request.