The genus Epipactis contains a large number of terrestrial orchids which vary greatly in distribution area and habitat type (Sramkó et al., 2019; Evans and Jacquemyn, 2020). Some species (e.g., E. dunensis and E. albensis) have very localised distributions and are restricted to particular habitats such as coastal dunes and beech forests, whereas others (e.g., E. helleborine and E. atrorubens) are widespread and can tolerate a relatively wide range of habitat conditions. There are several ecotypes of E. helleborine that can be found in specific habitats such as coastal dunes and forests (Jacquemyn et al., 2018). Species are autogamous, allogamous, or facultative allogamous (Claessens and Kleynen, 2011; Brys and Jacquemyn, 2016). The numerous seeds produced by Epipactis species are very small, dispersed by wind, and rely on the presence of mycorrhizal fungi in the soil to germinate and establish (Bidartondo and Read, 2008; Smith and Read, 2010; McCormick and Jacquemyn, 2014; Jacquemyn et al., 2018; Xing et al., 2020). Differences in mycorrhizal communities between localities may contribute to reproductive isolation and spatial distribution of Epipactis species and populations (Ogura-Tsujita and Yukawa, 2008; Jacquemyn et al., 2016, 2018; but see Těšitelová et al., 2012).

Records of each species’ occurrence from 2000 to 2020 on the continent of Europe were obtained from the online database GBIF¹ (Supplementary Material). We discarded records with missing GIS coordinates, ambiguous species identification or with coordinates with a spatial resolution lower than 100 m. This resulted in between 31 (Epipactis lusitanica) and 45,354 (E. helleborine) occurrences per species. Records for each species were aggregated into 10 km² grid cells to reduce the effects of spatial clustering resulting from sampling bias, by extracting the centre coordinates of each grid cell in which the species was recorded (Supplementary Table 1). Processing of occurrence data was performed in QGIS v3.4.9 (QGIS Development Team, 2019).

Previous studies have shown that land cover, bedrock, precipitation, and temperature are important variables predicting the distributions of some Epipactis species (Tsiftsis et al., 2008; Djordjević et al., 2016; Evans and Jacquemyn, 2020). We therefore used nine raster-format predictor variables with <0.5 correlation with one other. Two of the 19 bioclimatic variables available at the WorldClim v2 online database (Fick and Hijmans, 2017²) were used in our model, mean annual temperature and annual precipitation, projected for the near-present climate (1970–2000). These two variables were chosen because they are the most representative of the mean climate of an area, and are therefore appropriate for a continent-wide distribution study such as this. We also obtained the mean annual temperature and annual precipitation rasters predicted for the years 2061–2080 predicted by two Shared Socio-economic Pathways (SSPs), SSP 2-4.5 and SSP 5-8.5 from WorldClim.³ SSP 2-4.5 models the climate in a scenario where greenhouse gas emissions are at their highest (∼44 GT CO₂) in 2040 and then decrease to 9.6 GT in 2100, while in SSP 5-8.5, emissions increase steeply until the year 2080 (∼130 GT) before starting to stabilise and decrease (Riahi et al., 2017). Maps of the distribution of temperature and precipitation values in Europe were created by calculating the mean temperature for each cell of a 50 km² cell grid of Europe and summarising the values in QGIS.

The other seven variables used were the same as those used to model Epipactis species in Evans and Jacquemyn (2020). These include the first two components of two PCAs run on two topsoil datasets (physical and biochemical measures) acquired through the European Soil Data Centre (ESDAC) (Hiederer, 2013; Ballabio et al., 2019), dominant bedrock from the ESDAC database (Van Liedekerke et al., 2006), Corine Land Cover (CLC) from the Copernicus programme of the European Environmental Program (Heymann, 1994) and elevation (Amatulli et al., 2018). All raster processing was performed in RStudio v4.0.2 (R Core Team, 2021).

Defining and quantitatively comparing plant niches can be achieved using ecological niche models (ENMs). Ecological niche modelling has been applied successfully to numerous species to investigate ecological niches and to assess the impacts of climate change and land use on species ranges (Guisan and Thuiller, 2005). We used the programme Maxent v3.4.1. (Phillips et al., 2017) to model the effects of predicted climate change on species’ habitats. Maxent is a popular ENM tool that uses species occurrence data and environmental rasters to calculate a Gibbs value for each pixel of the study area, or the probability that the pixel has suitable habitat conditions for the species (Phillips, 2005) and performs well in comparison to other modelling methods (Elith et al., 2006; Phillips and Dudík, 2008; Valavi et al., 2021). Maxent creates habitat suitability maps over the study area from these data, as well as a table of the contribution of each predictor variable to the distribution of suitable habitat for each species. The choice of Maxent settings was informed by Barbet-Massin et al. (2012) and Merow et al. (2013). Each model was run using a random seed and 100 bootstrap replicates with 75% of the data used to train the model and 25% to test it. The rest of the settings were left as the default (convergence threshold of 0.00001, regularisation threshold of 1 and a maximum of 10,000 background points) and allowed for linear, quadratic, product and hinge features to be chosen automatically, producing a cloglog output. The models were run for the current climatic features and projected onto the SSP climate data to produce separate environmental niche map outputs for the current and future climate conditions.

For each habitat suitability map, the mean Gibbs value with standard error was calculated for every latitudinal interval of 0.5 decimal degree of the study area (Europe) using the Zonal Statistics tool in QGIS. We ran a Wilcoxon signed rank test on these mean Gibbs values multiplied by their corresponding latitudes to test whether the suitable habitat of each species will shift in latitude in future climate scenarios. Before running the Wilcoxon tests, the set of Gibbs values for each climate scenario was centred by dividing by the mean to eliminate the influence of different mean Gibbs values between climate scenarios and test only for shifts in latitude.

The continuous probabilistic maps produced by Maxent were converted into binary presence maps using the Maximum Test Sensitivity plus Specificity (MTSS) value of each species as a threshold. The pixels with a Gibbs value of greater than the MTSS were extracted and plotted as a new map for each species, with each pixel representing the species being present at that location. The total numbers of pixels occupied by these habitat distribution maps were compared between current and future climate scenarios. The overlapping pixels between the current and future distributions (i.e., pixels for which occurrence equalled one for both maps) of each species were extracted and counted to provide a measurement of the area suitable for a species if it were unable to expand into any newly available areas created by future climate change.

Mean species occurrence per climate scenario and per species, for both the continuous Gibbs values and the number of pixels occupied of the binary maps, were compared between the climate scenarios using Kruskal–Wallis and Dunn tests with a Holm correction or ANOVA and Tukey tests if the data were normally distributed.

We calculated Levins’ B2 values of niche breadth in geographic (B2_geo) and environmental space (B2_env) for each species using the functions raster.breadth and env.breadth, respectively, in the ENMTools 1.0.5 R package (Warren et al., 2021). B2 ranges from 0 to 1, with values closer to 0 representing narrow (specialised) niche breadth and values closer to 1 representing wide (generalised) niche breadth. Finally, we investigated whether species niche breadth predicts changes in distribution in response to climate change by dividing the range change of a species from current to future scenarios, converting to the proportional change for each species, and comparing these values to each species’ B2 value using ordinary least-squares regression. Ordinary least-squares regression was also used to compare the Levins’ B2 values between geographic and environmental space. The difference in niche breadth between species with positive range changes and negative range changes in response to habitat change was investigated using Kruskal–Wallis tests. The effect of mating system on response to climate change was tested by comparing the mean changes in latitudinal habitat distribution and proportional range size between autogamous, allogamous and facultative autogamous species, using Kruskal–Wallis tests. All analyses were performed in RStudio v4.0.2 (R Core Team, 2021).

Although there was no change in the mean Gibbs value of the 14 species between current and future climate scenarios, the Gibbs values for the majority of the species individually was predicted to increase in 2061–2080. The area of suitable habitat increased into the north for some species and decreased in the south for most species in the future, resulting in a mean northern shift in habitat. When expansion into the north (climate tracking) was restricted, the area of habitat decreased by up to 95% for all except two small-range species.

Despite the expectation that species with narrow environmental tolerances are most threatened by climate change, in the case of Epipactis, the habitats of most of the small-range localised species that we investigated were predicted to increase with future climate change. Some northern hemisphere herbaceous species benefit from increased temperatures at the northern edge of their distribution through increased population growth, which in turn can lead to an increase in geographic range at this edge (Bremer and Jongejans, 2009). This includes orchids such as H. hircinum where climatic changes in the United Kingdom were shown to be partially responsible for the species’ expansion between 1991 and 2001, as well as for projected future scenarios (van der Meer et al., 2016). Similarly, Ongaro et al. (2018) predicted that the habitat range for nine orchid species will increase by 2070 on the island of Sardinia, although the probability of presence in the newly colonised habitats was not predicted to increase. However, Vogt-Schilb et al. (2015) found that the distributions of many orchids in Western Europe have declined in the last two decades due to land-use change, particularly in the northern parts of their distributions. If land-use continues to change in more northern latitudes, this could limit the areas into which Epipactis can move in response to climate change. The results of our study provide further support for the potential for orchid ranges to increase at their leading edges in response to climate change, but go further to demonstrate that this does not necessarily mean an increase in available habitat, particularly if they cannot move into the northern habitats in time.

A widely supported paradigm is that the maximum range limits of a species coincide with its ecological niche limits and that, given the opportunity to disperse, range limits will shift to match the geographical extent of the niche under climate change (Reed et al., 2021; Shay et al., 2021). The pattern of species with wider niche breadths demonstrating greater latitudinal shifts in response to climate change has been documented in a number of terrestrial plant taxa (Thuiller et al., 2005; Alarcón and Cavieres, 2018). This was also demonstrated in orchids by Geppert et al. (2020) where the distributions of generalist orchid species and those inhabiting forests and semi-natural grasslands tended to be less affected than the more specialised and rare species in subalpine, natural grassland and wetland habitats, whose rear and leading edges shifted upward. This corresponds with our finding that Epipactis species with wider niche breadths (generalists) experience greater change in habitat area in response to changing climate than specialists. If we were to assign species to the groups of specialist and generalist based on current spatial ranges and values of B2 in geographic space, E. fageticola, albensis, kleinii, and lusitanica would be considered the most specialist (relative to the other species in this study), followed by E. tremolsii, dunensis, and leptochila as moderately specialist (Supplementary Table 2). E. muelleri, microphylla, and purpurata are moderately generalist, while E. helleborine, palustris, and atrorubens could be considered generalists. However, E. purpurata had a low B2 value but a fairly large spatial distribution and E. phyllanthes a high B2 value and small range, which is in contrast to this pattern. E. purpurata was predicted to experience a significant decrease in suitable habitat under climate change which may indicate that species with relatively large current ranges may still have fairly narrow niches which are nonetheless currently common in the environment, but are under threat from changing climate.

The distributions of all investigated species, even generalists, tended to lag behind climate warming, without being able to fully track the upward shift in suitable climate resulting in a range contraction, in both our study and Geppert et al. (2020). Plant species inhabiting forests may be somewhat buffered from the effects of climate warming (De Frenne et al., 2013; Zellweger et al., 2020) and those in grasslands tend to have high thermal ranges because of the lack of this buffering (Geppert et al., 2020). Similarly, Vogt-Schilb et al. (2015) found higher rates of disappearances in wetland orchid species in Western Europe than those in grassland, and more appearances in forest. There did not seem to be any clear pattern in response to climate change and habitat-use in our species (other than with niche breadth), with woodland species such as E. muelleri decreasing in suitable habitat and E. microphylla increasing. However, our study used a broad-scale specification of land cover, while more may be revealed at a finer resolution that captures microclimate gradients.

An important caveat to consider when carrying out niche breadth studies is that the metric used to describe niche breadth can greatly affect the results. Levin’s B2 is the reciprocal of Simpson’s diversity index (Levins, 1968) and has been a popular metric of niche breadth for more than 50 years. However, it has been noted that the traditional calculation of this metric is in geographic space (see Peterson and Soberón, 2012) for more on geographic and environmental space) and more accurately represents the “flatness” of the geographic distribution of suitable habitat (Warren et al., 2019), which may be a useful measure of spatial habitat-use but is not niche breadth in terms of specificity of resource-use. This is demonstrated clearly in our results, where B2 in geographic space was consistently correlated with measures of the size of the habitat distribution and the mean Gibbs value. B2 in environmental space as proposed by Warren et al. (2019) and developed in Warren et al. (2021) filters the geographic habitat suitability distribution through the set of environmental variables that was used to create the Maxent model, resulting in a B2 value that is closer to the concept of niche breadth as being the specificity in environmental conditions of a species’ habitat. It is important to note, however, that although closer to what we understand to be niche breadth, B2 in environmental space is still dependent on the availability of habitats in geographic space (Petraitis, 1979; Warren et al., 2019). Although the values of B2 in geographic and environmental space were not correlated, B2 in environmental space showed some evidence for having the same relationship with habitat changes as B2 in geographic space. This indicates that B2 in environmental space has the potential to be a useful representation of niche breadth for Epipactis in Europe, but further study is required to conclude this.

Although the abiotic characteristics discussed here are important for predicting orchid ranges, biotic interactions and species-specific characteristics are also essential contributors to the realised niche, and including these interactions can improve the accuracy and performance of niche models (Flores-Tolentino et al., 2020; Phillips et al., 2020). Orchids rely on insect pollinators and mycorrhizal fungi to reproduce and germinate (Rasmussen, 2008; McCormick and Jacquemyn, 2014). As with other pollinator-reliant plants, allogamous Epipactis will only persist and track shifting climate if their pollinators are also able to disperse (Benning and Moeller, 2019; Shay et al., 2021), such as has been predicted for a Neotropical orchid bee which is predicted to persist and increase its habitat range under future climate change (Silva et al., 2015). Mating system was not significantly associated with changes in habitat distribution in response to climate change, indicating that in the specific case of these species, autogamous and allogamous species did differ in response to predicted climate change. This is not surprising, considering that mating system was not significantly associated with niche breadth or range size for Epipactis species in previous studies (Evans and Jacquemyn, 2020) and niche breadth in geographic space is directly linked to range size. The presence of soil microbes has also been linked to the ability of plants to expand into newly available habitats (David et al., 2019; Bueno de Mesquita et al., 2020; Benning and Moeller, 2021; Shay et al., 2021). The diversity of mycorrhizal fungi is linked to latitudinal gradients for some orchid species (Duffy et al., 2019), but it is unclear whether the northern shifts in orchid distributions will be supported by the lower diversity of fungi in more northern latitudes. Our understanding and predictions of orchid distribution changes in response to climate change would be greatly improved with the addition of pollinator and fungal symbiont distributional data.

Studies that model the ecological niches of species are useful for conservation planning, particularly for identifying newly accessible areas available to plants (more so than predicting range contractions) and assessing the risk faced by populations as a consequence of their range size (Schwartz, 2012; Shay et al., 2021). This study provides a useful estimate of new areas into which Epipactis can expand, which in conjunction with more information on predicted land change in these areas, could be used in conservation schemes to allow the genus to flourish under climate change. However, it is important to assess the results in light of individual patterns in addition to drawing general conclusions. This is demonstrated in the contrast between mean change in Gibbs value (no change) and the change in Gibbs value for individual species, where a number of species were predicted to increase in the future and, the increase in habitat area for some species and the decrease for other. This disparity between general vs. individual species patterns has also been demonstrated in Geppert et al. (2020), who showed high interspecific variation among orchids grouped by habitat preference. We also show how some species with large areas of habitat such as E. purpurata should not be considered immune to the detrimental effects of future climate change as they may suffer considerable range reductions if they are not able to sufficiently disperse northwards.