Global warming not so harmful for all plants - response of holomycotrophic orchid species for the future climate change

Current and expected changes in global climate are major threat for biological diversity affecting individuals, communities and ecosystems. However, there is no general trend in the plants response to the climate change. The aim of present study was to evaluate impact of the future climate changes on the distribution of holomycotrophic orchid species using ecological niche modeling approach. Three different scenarios of future climate changes were tested to obtain the most comprehensive insight in the possible habitat loss of 16 holomycotrophic orchids. The extinction of Cephalanthera austiniae was predicted in all analyses. The coverage of suitable niches of Pogoniopsis schenckii will decrease to 1–30% of its current extent. The reduction of at least 50% of climatic niche of Erythrorchis cassythoides and Limodorum abortivum will be observed. In turn, the coverage of suitable niches of Hexalectris spicata, Uleiorchis ulaei and Wullschlaegelia calcarata may be even 16–74 times larger than in the present time. The conducted niche modeling and analysis of the similarity of their climatic tolerance showed instead that the future modification of the coverage of their suitable niches will not be unified and the future climate changes may be not so harmful for holomycotrophic orchids as expected.

The cluster analysis UPGMA of the studied holomycotrophic species in terms bioclimatic differences. Diagram was prepared in PAST 2.14 (https://folk.uio.no/ohammer/past/). suitability of climatic niche in areas where the species is known to occur. Models of P. schenckii, U. ulaei, W. calcarata, A. montana, G. lindleyana, and A. macranthus showed additional areas of high suitability in places where no records of these species were reported so far, however it is commonly accepted that realized niche is smaller than the climatic niche.
Factors limiting occurrence. The most limiting distribution factor of studied orchids was precipitation of the driest month (bio14). This variable was crucial for occurrence of six studied species: H. spicata, W. calcarata, D. hamiltonianum, E. cassythoides, N. nidus-avis, and E. aphyllum. The precipitation of the coldest quarter (bio19) significantly contributed in the models for C. austiniae, U. ulaei, and G. sesamoides. Occurrence of two species, C. wisteriana and P. schenckii, depends mostly on the maximum temperature of the warmest month (bio5). The annual precipitation (bio12) was the most important limiting factor for A. montana, the precipitation of the warmest quarter (bio18) for G. lindleyana, and the mean diurnal range (bio2) for A. macranthus. The complete list of the variables contributing to the models of studied species is given in Supplementary Table S2.
Variation of climatic preferences. The bioclimatic variation of the studied orchids visible on the PCA analysis diagrams was large. The specimens have created almost linear pattern of grouping along environmental gradient (Fig. 1A). The first axis distinguished G. lindleyana, E. aphyllum, H. spicata (above the axis) from P.      schenckii and L. abortivum (below the axis). Furthermore, the second axis separated E. aphyllum and C. wisteriana from P. schenckii, U. ulaei, W. calcarata, A. montana, D. pallens, D. hamiltonianum, E. cassythoides, G. sesamoides, and A. macranthus on the left side of the diagram. Variables with the greatest contributions enhanced in the overall analysis were: the temperature seasonality (bio4); the annual precipitation (bio12); the precipitation of the wettest month (bio13); the precipitation of the warmest quarter (bio18); the precipitation of the coldest quarter (bio19). In turn, the canonical variate analysis (CVA) showed a highly statistically significant differentiating value for the holomycotrophic species with respect to the bioclimatic factors (Wilks λ = 0.001; F (180;12807) = 75.381; P < 0.0001). The canonical variate axes (52.96 and 22.35% of explained variance, respectively) have distinguished individual species from each other along environmental gradient obtained on the basis of the bioclimatic factors (Fig. 1B).
The major share in the discrimination of the studied species had the temperature seasonality (bio4) and the annual precipitation (bio12) which discriminated C. austiniae and C. wisteriana. Isothermality (bio3 = −1.197), the precipitation of the coldest quarter (bio19 = 0.704), as well as, the annual mean temperature (bio1 = 3.137), the mean diurnal range of the temperature (bio2 = 1.971) and the maximum temperature of the warmest month (bio5 = −2.983) had also the significant meaning in the performed discrimination of the studied orchids. The bioclimatic differences between the pairs of species were evaluated by permutation tests of the climatic factors. The interaction term in the MANOVA (Pillai's trace = 3.843; F (180;16380) = 42.88; P < 0.0001) was significant. The      Two clades were recognized in the UPGMA dendrogram ( Fig. 2) illustrating the similarity between analysed orchids in the preferred niches. Both included species from different geographical regions.    Climatic niche loss will be observed in case of G. lindleyana in A2a and B2a scenarios while E. aphyllum will lose niches exclusively in B2a scenario.
The negative effect of climate change will be remarkable for C. austiniae and P. schenckii. The extinction of C. austiniae was predicted in all analyses. The coverage of suitable niche of P. schenckii will decrease to 1-30% of its current extent. The reduction of at least 50% of climatic niche coverage of E. cassythoides and L. abortivum will be observed. On the other hand, the coverage of suitable niches of H. spicata, U. ulaei and W. calcarata may be even 16-74 times larger than in the present time.

Discussion
Predicting the response of biodiversity to climate change became an important aspect of nature conservation 2,25-28 . The computed models are crucial for warning decision makers to potential future risks and support the establishment of active strategies to reduce negative impacts of climate change on organisms 26,29,30 . Climate change scenarios are the base for predictive analysis and these depend on a various socio-economic storylines for greenhouse gas emissions 26 and on a broad range of General Circulation Models used to calculate climate change for given trajectories of greenhouse gas emissions 30 . For this reason the projections of habitat loss can gave a contrasting results depending on the choice of combinations of emissions scenarios and climate models 31 . Bellgard & Williams 32 partitioned the global climate change into four putative time episodes. According to their proposal the long-term impacts (over 21-50+ years) are related with increased temperatures and CO 2 that will destabilize global rainfall patterns, soil properties and plant ecosystem resilience. Because of dependence on their host for C-supply, orchid mycorrhizas and all heterotrophic mycorrhizal groups will be immediately impacted through loss of habitat and plant-hosts.
In our study three different scenarios of future climate changes were tested to obtain the most comprehensive insight in the possible climatic niche loss of 16 species of holomycotrophic orchids. No general trend in the response of holomycotrophic orchid for climate changes was identified. We did not recognize any correlation between climatic factors currently limiting the distribution of the studied orchids with their response to climate changes. The first of two clades recognized in the UPGMA dendrogram included 12 species from different geographical regions (Europe, Asia, Australia, South and North America) which suitable niche coverage will be variously modified in the future. The second clade included species which potential niche coverage will increase as a result of climate change -two South American species (W. calcarata and U. ulaei), African A. macranthus and Asian D. pallens. This analysis indicate that species similar in general climatic preferences will respond to the global warming in different way.  Obviously, the actual future range of studied orchids may be slightly different than presented in our models. The availability of mycorrhizal fungi and presence of pollinators may restrict the distribution of these species while the niche shift may allow them to occupy new habitats. The current state of knowledge on the specify of holomycotrophic orchid symbionts and their dependence on particular pollen carriers is unfortunately too incomplete to include these elements in the analyses. Moreover, such complex ecological process are not possible to evaluate using currently available statistic approaches. However, we do believe that climate is the key factor limiting occurrence of most organisms by determining characteristics and distributions of natural and managed systems. Hereby modeling of distribution of climatic niches is the first step to evaluate the response of particular species to the possible climate changes.
It is not clear how and whether the function of mycorrhizal fungi in orchid germination and growth will be maintained with rising temperature, erratic rainfall and reduced moisture. The research conducted on photosynthetic, terrestrial orchids from Europe -Anacamptis morio (L.) R.M. Bateman, Pridgeon & M.W. Chase and Dactylorhiza fuchsii (Druce) Soó indicated that orchid rarity and persistence are not necessarily related to fungal diversity and that other factors may be more important in determining orchid survival chances 33 . Furthermore, the modifications of the soil structure which is especially important for terrestrial plants are extremely difficult to evaluate due to multiple factors which can influence the ground properties. Studies on European plants growing in understory layer indicated that numerous species are able to shift their realized niches, also for soil nutrients and pH 34 . The climate change (i.e. CO 2 enrichment, arise temperature, altered precipitation, increased N-deposition) impact the soil-rhizosphere, plant and fungal physiology and/or ecosystem(s) directly and indirectly. Direct effects include changes in resource availability and change in distribution of mycorrhizas. Indirect effects include changes in below ground allocation of C to roots and changes in plant species distribution 32 .
As mentioned before the pollinators availability can limit the long term viability of orchid populations, but this factor could not be included in our analyses due to the lack of sufficient data on pollinator specificity of studied plants. We know at least five species of insects that can pollinate E. aphyllum in Poland and the Czech Republic 35 , but the complete set of species which are able to transfer pollinia of this orchid remains unknown, in turn, the recent study proved that the turnover in the composition of insect communities may occur in the future 36 . On the other hand, almost no data on pollinators of tropical Orchidaceae has been published. The estimation of the impact of future climate changes on the pollinators distribution in not the only problem because, as indicated by Robbirt et al. 37 climate changes can also led to asynchrony in flowering of orchids and insect flight periods. The future modification in insects phenology are related with early adult emergence [38][39][40] , earlier larval emergence 41 , and earlier migration 40,42 . Additionally, the latitudinal 43,44 and altitudinal 45,46 shifts are predicted as a result of climate change. Also, the global expansion of tropical species into temperate areas is possible 47,48 .
Orchids as a group show evolutionary flexibility whereby diversification in the family is often related to habitat complexity and fragmentation. Throughout their history, Orchidaceae have been able to cope with face of climatic change caused by shifting continents, mountain uplifting, fluctuating sea levels and temperatures 49 . All these phenomena occur today but the rate of change seems to be occurring faster than the detectable past with available methods 50 . Some species experiences loss of habitat, population reduction and cannot keep pace with the climate change 51,52 . Liu et al. 53 estimated that populations of at least 15% of the orchid species in a diverse region of southwestern China will be threatened with extinction over the next two centuries given projected climate changes. Also, Nadkarni & Solano 54 and Olaya-Arenas et al. 55 indicated drying trend in some cloud forests, including Costa Rica, which might affect those species such as the hundreds of Neotropical Lepanthes representatives that depend on cool and wet conditions.
Our studies indicated that the future climate changes may be not so harmful for all holomycotrophic orchids as expected 8 . Out of the 16 studied species, the potential suitable niches was not negatively affected in all the models. In our opinion the most significant threat for the plants is the direct human activity that results in habitat destruction and fragmentation or modification of the ecosystem net. The other conclusion of our study is that generalization of the effect of global warming is misleading. Not all organisms, even related, occurring in the same geographical region and preferring similar climatic conditions will face the same modifications of their suitable niches. The object of our study were holomycotrophic orchids that in theory should respond in a similar way for a global warming. The conducted niche modeling and analysis of the similarity of their climatic tolerance showed instead that the future modification of the coverage of their suitable niches will not be unified. Over 1500 localities of studied orchids were compiled based on available herbarium material, literature data and electronic databases (Supplementary Dataset S1). Only those localities which could be precisely localized on the map were used. From the database the duplicate records were removed. To reduce the sampling bias we removed also closely lying localities. Ecological niche modelling. While a broad range of algorithms 58,59 and platforms (i.e. BIOMOD 60 , ModEco 61 , openModeller 62 ) can be used to produce species distribution models, in this study the maximum entropy method implemented in MaxEnt v. 3.3.2 63-65 based on the species presence-only observations was used. This application has been proved to provide the most robust response across the number of environmental variables tested 66 and it has been shown to work better with small number of samples than other approaches 67 68 and provided by WorldClim (v. 1.4 release 3; www.worldclim.org) we removed seven variables due to their significant correlation. For the same reason we did not use altitude as input data. The variables included in the analysis are listed in Supplementary Table S4. To improve model performance target-group background was selected for each continent 69 . Additionally for P. schenckii we reduced background to avoid significant over-estimation of the model (Supplementary Table S5).

Species selection and list of localities.
In all analyses the maximum iterations was set to 10000 and convergence threshold to 0.00001. The "random seed" option was applied, this provides random test partition and background subset for each run was applied. A bootstrap procedure with 1000 replicates was applied and the output was set to logistic. All operations on GIS data were carried out on ArcGis v. 9.3 (ESRI).
To estimate the impact of hypothetical climate changes on the coverage of suitable niches of studied plants we used the same climatic variables as considered in the present models. The layers based on the coupled global climate models were used (http://ccafs-climate.org). Three various emission scenarios for 2080 were analysed: A1b (CCCMA-CGCM3 SRES simulation), A2a (CCCMA-CGCM2 SRES) and B2a (CCCMA-CGCM2 SRES). This approach was used in numerous recent studies that focused on the climate change impact on the distribution of various organisms [70][71][72] . While we are aware that recently several new models predicting future climate change have been published 73 , we used aforementioned datasets to compare our results with previously published studies focused on changes in coverage of climatic niche of Orchidaceae 18 .

Multivariate analysis.
Principal components analysis (PCA) was performed to explain the general variation pattern among the studied holomycotrophic species, based on 12 bioclimatic factors used in ENM analysis (Table S4) and altitudinal data. To determine the bioclimatic factors which differentiate the studied orchids the most, the canonical variate analysis (CVA) was applied, in order to reduce the data set by selecting only the factors that showed strongest discrimination. The significance of bioclimatic differences between the studied species was evaluated using Wilk's λ and Goodall's F-ratio with 1000 permutations. In turn, the Mahalanobis distance matrix was prepared to classification of examined species and was used in the cluster analysis UPGMA (unweighted pair-group average method) 74 . Statistical computations were performed with the program PAST v. 2.14 75 .
SCieNTiFiC REPORts | 7: 12704 | DOI:10.1038/s41598-017-13088-7 Compliance with Ethical Standards. The authors declare that there are no conflicts of interest. This article does not contain any studies with human participants or animals performed by any of the authors.