Introduction

So far little is known about biogeography of Pleurothallidinae which is one of the largest taxon of Orchidaceae and currently includes over 5100 species (^{Karremans 2016}) distributed from USA (Florida) to southern Brazil. The critical gap in our knowledge concerns historical biogeography. So far macroremains of just two orchid genera, Dendrobium Sw. and Earina Lindl. (^{Conran et al. 2009}) were identified in fossil material. The secondary footprint of those plants, pollenbearing insect, preserved in amber was described in 2007 (^{Ramírez et al. 2007}). Insufficient palaeobotanical data limited numerous phytogeographical studies, especially research on post-glacial migration of orchids. Another obstacle in the research on orchids was related to the substantial deficiency in the information about their ecology, especially habitat requirements of tropical taxa. For a long time difficulties in gathering data on their preferred niches limited studies on the impact of future climatic changes on their potential distribution which is crucial for establishing long-term nature conservation plans. First comprehensive mathematical models of future climatic changes were developed in the mid-XX century, but their usage was possible only for taxa with recognized ecological tolerance, or at least climatic niche requirements.

The comprehensive biogeographical studies on Orchidaceae became possible with the development of species distribution models (SDMs). Initially the ecological studies on relationships between species and its habitat, relied largely on linear multiple regression and discriminant function analyses (^{Capen 1981}). SDMs improvement was related to new regression methods that provided coherent treatments for the error distributions of presence-absence and abundance data (^{Elith & Leathwick 2009}). The generalized linear models (GLMs) enabled pioneering regression-based SDMs which continue to be useful and are part of many current methods including maximum entropy models (MaxEnt; ^{Phillips, Anderson & Schapire 2006}). In orchid studies SDMs were first used to evaluate status of invasive taxa (^{Kolanowska 2013}a, Kolanowska & Konowalik 2014), but quickly they allowed studies on glacial distribution of various species (e.g. Kolanowska 2013b, 2015, ^{Szlachetko et al. 2014}). Recently, first estimates of phylogenetic niche conservatism and evolution of climatic tolerances within Apostasioideae were published (Kolanowska et al. 2016).

The goal of our study was to present comprehensive biogeographical data on the genus Diodonopsis ^{Pridgeon & M.W. Chase which was described in 2001} based on molecular studies results (Pridgeon & Chase 2001) and currently is broadly accepted (e.g. ^{Pupulin 2002}, Pridgeon 2005, Kolanowska & Szlachetko

2014). This taxon initially included five representatives of Masdevallia sect. Pygmaeae Luer, however, as suggested by Luer (2000) and proved in the subsequent studies, the genus in the treatment provided by ^{Pridgeon & Chase (2001}), is not monophyletic (^{Abele 2007}, ^{Matuszkiewicz & Tukallo 2006}). In 2006 Luer (2006) segregated from Diodonopsis two species, D. hoeijeri (Luer & Hirtz) Pridgeon & M.W. Chase and D. pterygiophora (Luer & R. Escobar) Pridgeon & M.W. Chase which were included into new genus named Pteroon Luer. The two taxa are basically distinguished based on ovary structure - soft, burr-like in Diodonopsis and ornamented with plate-like keels overlapping the bases of the sepals in Pteroon. Currently Diodonopsis is composed of just three species. All are caespitose plants with short ramicauls and slender peduncles. Their sepals are tailed and they form a short tube or cup. D. anachaeta (Rchb.f.) Pridgeon & M.W.Chase and D. erinacea (Rchb.f.) Pridgeon & M.W.Chase produce solitary flowers and D. pygmaea (Kraenzl.) Pridgeon & M.W.Chase is successively flowered. The aim of this research was to describe distribution patterns of Diodonopsis representatives, reveal the factors limiting their distribution as well as to evaluate differences in their preferred climatic niches. Moreover, using ecological niche modeling we reconstructed distribution of the glacial refugia of the three species and we estimated the possible changes in the potential range of particular species in the future.

Material and Methods

Distribution patterns and niche similarity.- Alpha diversity of Diodonopsis representatives was explored using DIVA-GIS version 7.5. The information about climate and vegetation type preferences of studied species were evaluated using Global Ecological Land Units developed by Esri and the U. S. Geological Survey (USGS) in 2015 (^{Sayre et al. 2015}). The global map of terrestrial ecoregions (^{Olson et al. 2001}) was used to estimate species compositional similarity. The matrix for this analysis was created based on binary coding systems. UPGMA algorithm was used to create dendrograms and the distance matrix was computed using Jaccard similarity for absence-presence data.

The differences between the climatic niches occupied by three studied species were evaluated using the niche identity and overlap test indexes: ^{Schoener’s D and I statistic as available in ENMTools v1.3 (Schoener 1968}, ^{Warren, Glor & Turelli 2008}, 2010) and illustrated using linear discriminant analysis (LDA). As input data for these analyses the climatic variables (Tab. 1) were used. Statistical computations (UPGMA, LDA) were performed using PAST v. 3.14 (^{Hammer, Harper & Ryan 2001}).

Ecological niche modeling.- The ecological niche modeling (ENM) was conducted using maximum entropy method implemented in MaxEnt version 3.3.3k (^{Phillips, Dudík & Schapire 2004}, Phillips, Anderson & Schapire 2006, ^{Elith et al. 2011}) based on the species presence-only observations. The list of localities was compiled based on examination of material deposited in COL, K, MO, as well as data provided by ^{Dodson & Luer (2009}) who additionally revised specimens in AMES, G, GUAY, K, SEL, and W. GoogleEarth software was used to gather coordinated for the occurrence data provided on herbarium labels. A total of 48 locations of Diodonopsis were gathered (Annex 1) which is more than the minimum number of records required to obtain reliable predictions in MaxEnt application (^{Pearson et al. 2006}, Wisz et al. 2008).

From 19 climatic variables (“bioclims”) in 2.5 arc minutes (± 21.62 km2 at the equator) developed by ^{Hijmans et al. (2005}) and provided by WorldClim (version 1.4 release 3, www.worldclim.org) we removed seven variables due to their significant correlation (above 0.9) as evaluated by the Pearsons’ correlation coefficient calculation computed using ENMTools v1.3. The following variables were excluded from the dataset: bio6, bio7, bio9, bio10, bio11, bio16 and bio17. All analyses was made for the area restricted to 37.83°S-26.125°N and 107.95-27.6°W. To reconstruct the distribution of suitable niches of studied species during the last glacial maximum (LGM, ca. 21 kya) the CCSM4 simulation was used. The original data was made available by available by Coupled Model Intercomparison Project Phase 5 and later downscaled and calibrated using WorldClim 1.4 as baseline “current” climate (http://worldclim.org/paleo-climate1; see ^{Braconnot et al. 2012}).

To estimate the impact of a hypothetical climate change climatic variables considered in the present models as well as data from simulations with the second and the third Generation Global Climate Model (CGCM2 and CGCM3) provided by the Canadian Centre for Climate Modelling and Analysis (CCCMA; spatial data available at http://ccafs-climate.org and http://climate-modelling.canada.ca/data/data. shtml) were used. Three various emission scenarios for 2080 (the longest time horizon available) were analysed: A1b (CCCMA-CGCM3 simulation, Delta Climgen dataset), A2a (CCCMA-CGCM2, Delta

Method IPCC AR3 dataset) and B2a (CCCMACGCM2, Delta Method IPCC AR3 dataset). The same approach was previously used to estimate the impact of climatic change on distribution of Dactylorhiza (Nacz & ^{Kolanowska 2015}) and Epipactis helleborine (Kolanowska 2013) in Europe. The characteristic of forcing agents in various future scenarios was included in the Intergovernmental Panel on Climate Change (IPCC) Special Report on Emission Scenarios (SRES). The A1b storyline describes a future world of very rapid economic growth, global population that peaks in mid-century and declines thereafter. In this scenario the balance across all energy sources is expected. Both A2a and B2a scenarios describe a “regionalisation” leading to a heterogeneous world development. The main driving forces in A2a are a high rate of population growth, increased energy use, land-use changes and slow technological change. In B2a a general evolution towards environmental protection and social equity is expected.

Figure 1 Diodonopsis species richness map. 

In all analysis the maximum iterations was set to 10000 and convergence threshold to 0.00001 forcing the program not to finish before the threshold was reached. The “random seed” option which provided random test partition and background subset for each run was applied. The run was performed as a bootstrap with 1000 replicates which is a valid method to assess model reliability with small sample sizes (Pearson et al. 2007, ^{Botero-Delgadillo, Páez & Bayly 2012}). The output was set to logistic. All operations on GIS data were carried out on ArcGis 9.3 (ESRI). The area under the Receiver Operating Characteristic curve (AUC) for each model was calculated to estimate the reliability of the analysis.

Results

Distribution patterns.- The largest geographical range has D. anachaeta which occurs from Bolivia to Ecuador. D. erinacea was reported from Ecuador, Colombia, Panama, and Costa Rica and D. pygmaea from Ecuador, Colombia, and Costa Rica. Populations of Diodonopsis are rather isolated and only in southern Ecuador, area near Ecuadorian-Colombian border, and northern Costa Rica two species occur close to each other (Fig. 1). The broadest elevation range has D. anachaeta (Fig. 2) which occurs in highlands (590 m) as well as in high mountains (3400 m) while the most restrict altitudinal range is observed in D. pygmaea (335-1800 m).

Figure 2 Elevation ranges of studied species. 

Species of Diodonopsis occur in six ecoregions: Bolivian Yungas, Eastern Cordillera Real montane forests, Isthmian-Atlantic moist forests, Northwestern Andean montane forests, Peruvian Yungas, and Talamancan montane forests. Diodonopsis anachaeta is found in four of them (Bolivian Yungas, Eastern Cordillera Real montane forests, Northwestern Andean montane forests, Peruvian Yungas), Diodonopsis erinacea in three (Eastern Cordillera Real montane forests, Isthmian-Atlantic moist forests, Talamancan montane forests), and Diodonopsis pygmaea in just two (Northwestern Andean montane forests, Talamancan montane forests). Two groups can be demarcated in the dendrogram presenting species composition similarity between ecoregions (Fig. 3). The first one includes Talamancan montane forests which is sister to the pair Isthmian-Atlantic moist forests and Eastern Cordillera Real montane forests.

Figure 3 Similarity of specific composition between ecoregions - dendrogram (UPGMA) based on the Jaccard similarity index. 

Figure 4 Current distribution of the suitable climatic niches of Diodonopsis species. A - D. anachaeta, B - D. erinacea, C - D. pygmaea. 

The second group is composed of Northwestern Andean montane forests and Yungas (Bolivian and Peruvian).

The highest variance of habitat is recorded in D. anachaeta which occur in grassland, scrub, or shrub, mostly needleleaf/evergreen forest, swampy or often flooded vegetation in cool, hot or warm wet mountains. D. erinacea grows exclusively in mostly needleleaf/ evergreen forest and D. pygmaea is observed in croplands as well as in mostly needleleaf/evergreen forest. Both these species occur in hot or warm wet mountains.

The potential range of Diodonopsis representatives calculated using MaxEnt is slightly larger than actually known geographical range of the particular species (Fig. 4). Suitable habitats of all three species are present in Guyana Shield, Guatemala, and southern Mexico. Additionally, D. anachaeta could also occur in Colombia and Costa Rica and D. pygmaea in Peruvian Andes and Nicaragua. The later species was not reported from Panama so far, however the coverage of suitable niches is high in this country.

Limiting factors and niche similarity.- The most important limiting factors for D. anachaeta are temperature seasonality (bio4, 24.1-29.6%), precipitation of the warmest quarter (bio18, 17.920.8%), and the annual mean temperature (bio1, 14.8-20.4%). Range of D. erinacea is limited by mean temperature of the wettest quarter (bio8, 17.920.8%), annual precipitation (bio12, 14-16.7%), and precipitation of the wettest month (bio13, 12.615.4%). The highest contribution to the models created for D. pygmaea gave temperature seasonality (bio4, 23.9-31.4%), maximum temperature of warmest month (bio5, 22.7-25.8%), and isothermality (bio3, 11-18.2%) (Tab. 2). The results of LDA analysis are presented in Fig. 5.

The most dissimilar niches are occupied by D. anachaeta and D. pygmaea (Tab. 3), however geographically the most separated are potential ranges of D. anachaeta and D. erinacea (Tab. 4). The range overlap was calculated with the suitability threshold for presence set as 0.5. The overlap between D. anachaeta and D. erinacea is only 0.0702, between D. anachaeta and D. pygmaea 0.560, and between D. erinacea and D. pygmaea 0.784.

Figure 5 Results of LDA analysis. The axes mark linear functions th discriminate among the populations of Diodonopsis species on the basis of climatic conditions. Axis1 - 96.03%, Axis 2 - 3,967%.at 

Glacial refugia and future potential range changes.- The glacial ranges of all studied species were wider than in the present time (Fig. 6). Refugia of D. anachaeta occurred east of south Peruvian and Bolivian Andes, in Peru Central Andes, on both slopes of Andes of Ecuador and southern Colombia. Moreover, suitable habitats were distributed in Guyana Shield and places which are currently under water, south west Brazilian coast and French Guiana coast. Glacial refugia of D. erinacea were more scattered, most of them were located in Ecuadorian and Colombian Andes and Venezuelan Gran Sabana. The potential habitats were also available in Sierra Madre de Chiapas, Costa Rican Talamanca range, Cordillera Guanacaste, Panamanian Darién region. Outside current land this species could survive in French Guiana coast. Suitable niches of D. pygmaea were very broadly distributed - in the Andes from southern Peru to Colombia, Santa Marta, in Venezuelan Gran Sabana, south-eastern Guyana Shield, Costa Rican Talamanca. As in previously mentioned species suitable niches were also located along French Guiana coastal line.

Future changes.- Generally the future climate changes will be harmful for Diosonopsis (Fig. 7). D. anachaeta will lost suitable niches in Bolivia and in Guyana Shield, in A1b scenario also southern Peruvian region will not be included in the species potential range. The geographical distribution of D. erinacea habitats will be generally the same, except of Guyana Shield. In other places the niches coverage will decrease. Surprisingly, additional niches preferred by D. pygmaea will become available in Sierra Madre de Chiapas, but in all other regions the potential range of this species will shrink.

Table 3 Similarity of climatic niches occupied by Diodonopsis representatives - results of niche identity tests (I and D). 

Figure 6 Models of glacial distribution (ca. 21 kya) of the suitable climatic niches of Diodonopsis species. A - D. anachaeta, B - D. erinacea, C - D. pygmaea. 

Figure 7 Models of future distribution of the suitable climatic niches of Diodonopsis species in 2080 based on A1b (balance between fossil and non-fossil energy sources), A2a (high energy requirements), and B2a (lower energy requirements) climate change scenarios. A, D, G - D. anachaeta, B, E, H - D. erinacea, C, F, I - D. pygmaea. 

Figure 8 Changes in the suitable climatic niches coverage of Diodonopsis species since LGM to 2080 according to various climate change scenarios (A1b, A2a, B2a). 

The coverage of suitable niches of all studied species are declining since LGM (Fig. 8). In the most harmful scenario of future climate change the potential range will be only 4.5% of the glacial range (Tab. 5). The highest habitat loss after glacial period is observed in D. anachaeta and in 2080 it will continue to decrease, even to 23% of the current potential range. D. pygmaea and D. anachaeta have greater chances to survive as in A2a and B2a scenarios the coverage of their suitable niches will not decrease for more than 60% of the current extent.

Discussion

Climatic niches of all Diodonopsis species are well-separated both in the aspect of their quality as well as in geographic distribution. ^{Givnish et al. (2015}) linked extraordinary diversity of orchids with the evolution of pollinia and, the epiphytic habit, CAM photosynthesis, tropical distribution, and pollination by Lepidoptera or euglossine bees, often by deceit pollination. However, we believe that the climatic niche tolerance should be studied more intensively to estimate its impact on orchids speciation.

As indicated by ENM analysis results currently several areas characterized by climatic conditions suitable for Diodonopsis species are located outside the known geographical ranges of the genus representatives. One of the reasons for which studied orchids do not grow there is that these locations, e.g. Guyana Shield was located too far from the centre of speciation and it was actually never inhabited. The other possibility that additional ecological limiting factors, e.g. lack of pollinator or mycorrhizal fungi, make Diodonopsis unable to survive in these climatically suitable regions.

The specific composition and similarity of ecoregions is a direct result of location of glacial refugia of studied species. The southernmost ecoregions, Bolivian Yungas and Peruvian Yungas served as refugial areas exclusively for D. anachaeta. The specific composition in the Northern Andes is more complicated. All three species could potentially occur along Andean slopes, but currently D. anachaeta and D. pygmaea occur in north-western Andes while D. anachaeta and D. erinacea grow in eastern Andean region. The only species with glacial suitable habitats in Panamanian Isthmus was D. erinacea. Talamancan region is characterized by a unique combination of species - both D. erinacea and D. pygmaea are present in this area. The former orchid shows huge distribution gap in the northern Andes despite it had potentially available niches in this region during LGM.

Table 5 Coverage of suitable niches of Diodonopsis species. 

The reconstruction of the history of organisms with narrow, scattered distribution and few observation records is laborious and generally two reasons of small geographical ranges are most often considered - range regression or very specific habitat requirements combined with effective dispersal capacity (^{Cassel-Lundhagen 2013}). Undoubtedly, not only climatic elements determine occurrence of Diodonopsis species. Numerous research aimed to recognize variables affecting the distribution and abundance of various orchids (e.g. ^{Bergman et al. 2006}, ^{O’Malley 2009}, ^{McCormick et al. 2009}, ^{Bunch et al. 2013}, ^{García-González et al. 2016}) indicating the importance of geography, soil chemistry, phorophyte specificity, presence of mycorrhizal fungi and pollinators, landscape heterogeneity, as well as landuse history. However, climate is the most important component of fundamental niche which is formed as a result of evolution. It is commonly accepted that the fundamental niche is broader than the realized niche (^{Severtsov 2012}). Also in our models of the present distribution of climatic niches of Diodonopsis species areas characterized by suitable conditions, where the genus representatives were not found, were recognized. Considering only potential occurrence, the results of our studies suggest that in case of Diodonopsis the known populations are only remnant of previously broader geographical ranges. The postglacial regression was recognized in other South American orchid - Vargasiella C.Schweinf. (^{Szlachetko et al. 2014}). ^{Bellard et al. 2012}) recognized three types of possible response to climate change: time (e.g. phenology), space (e.g. range) and organism itself (e.g. physiology), but so far only shifts in spatial distribution of orchids were indicated (e.g. ^{Reina-Rodríguez et al. 2017}). We assume that while suitable niches of all three species occurred in the Guyana Shield (and some of these exist also currently), the diversification centre of the genus was localized in the Andes and the short distance dispersal as well as significant geographical gap prevent it for expansion to this region. Unfortunately, the future climatic changes will not be beneficial for Diodonopsis and the whole genus will face huge risk of extinction. The possible places where the species can survive Andes on Ecuadorian/Colombian border and possibly in Costa Rican mountains (D. erinacea and D. pygmaea). However, the deforestation process in Ecuador is intensive (e.g. ^{Tapia-Armijos 2015}) and habitat fragmentation will further reduce the survival chances of Diodonopsis representatives.

From a phylogenetic point of view, it would be interesting to conduct niche modeling also for two Pteroon species and to combine it with results of molecular analyses to visualize the evolution of climatic tolerance within Pteroon-Diodonopsis group and to estimate its importance in the evolution.