Evolution of the climatic tolerance and postglacial range changes of the most primitive orchids (Apostasioideae) within Sundaland, Wallacea and Sahul

Data sources

A database of the sites of Apostasioideae representatives was prepared based on the examination of herbarium specimens stored in BM, K, MO, and P, as well as the electronic databases of AAU, BRI, CANBR, CNS, NOCC, SCBI, SING, US and the Swiss Orchid Foundation. The herbaria acronyms follow Index Herbariorum (Thiers, 2015). Only the sites that could be precisely placed on the map were used in our analyses. The process of georeferencing followed Hijmans et al. (1999). The geographic coordinates provided on the herbarium sheet labels were verified. If there was no georeferencing information on the herbarium sheet label, we followed the description of the collection place and assigned coordinates as precisely as possible.

We did not find enough material (occurrence data) to conduct analysis for Apostasia latifolia Rolfe, Apostasia parvula Schltr., A. ramifera S.C. Chen & K.Y. Lang, A. shenzhenica Z.J. Liu & L.J. Chen, Neuwiedia elongata de Vogel, N. griffithii Rchb.f., N. inae de Vogel, N. malipoensis Z.J. Liu, L.J. Chen & K. Wei Liu, or N. siamensis de Vogel. We gathered a total of 179 occurrence records of three Apostasia and three Neuwiedia species (Table S1): 27 of Apostasia nuda R. Br., 22 of A. odorata Blume, 80 of A. wallichii R. Br., eight of Neuwiedia borneensis de Vogel, 24 of N. veratrifolia Blume, and 18 of N. zollingeri Rchb. f. To reduce sample bias we applied spatial filtering (Boria et al., 2014). From a total of 179 records, we randomly removed sites of each species that were within 25 km of one another. The final database included 13 sites of Apostasia nuda, 13 of A. odorata, 41 of A. wallichii, four of Neuwiedia borneensis, 18 of N. veratrifolia, and 10 of N. zollingeri.

Ecological niche modeling

The ecological niche modelling was conducted using the maximum entropy method implemented in MaxEnt version 3.3.2 (Phillips, Dudik & Schapire, 2004; Phillips, Anderson & Schapire, 2006; Elith et al., 2011) based on the species presence–only observations. From 19 climatic variables (“bioclims”, Table 1) in 2.5 arc minutes (±20 km2 at the equator) developed by Hijmans et al. (2005) we removed seven “bioclims” due to their significant correlation (above 0.9) as evaluated by the Pearson’s correlation coefficient calculation computed using ENMTools v1.3. The following variables were excluded from the dataset: bio6, bio7, bio9, bio10, bio11, bio16 and bio17. As it was suggested that using a restricted area in ENM analysis is more reliable than calculating habitat suitability on the global scale (Barve et al., 2011), the area of our analysis was restricted to longitude 88°–190°E and latitude 53°S–47.8°N. The maximum iterations was set to 10,000 and the convergence threshold to 0.00001. The “random seed” option which provided the random test partition and background subset for each run was applied. The run was performed as a bootstrap with 1,000 replicates, and the output was set to logistic. The analogical settings were used in the modelling for LGM. The bioclimatic data for this time period were developed and mapped (CCSM4) by Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor, Stouffer & Meehl, 2012). All operations on GIS data were carried out on ArcGis 9.3 (ESRI).

The evaluation of the models was done using the most common metric—area under the curve (AUC; Mason & Graham, 2002) which was calculated by the Maxent application automatically. The niche overlap between the studied species was calculated using ENMTools v1.3 (Warren, Glor & Turelli, 2010).

Evolution of climatic tolerance and predicted niche occupancy (PNO)

To visualize niche evolution for Apostasioideae the ancestral tolerances were pictured on a phylogenetic tree. To construct the tree for taxa included in the ecological niche modeling sequences from ITS, matK and trnL available from GenBank (Table S2) were aligned in Mafft 6.833b (Katoh & Toh, 2008) and Gapcoder (Young & Healy, 2003) was used to code indels. Aligments were merged and a Bayesian phylogenetic analysis was performed in MrBayes 3.2.1 (Ronquist et al., 2012). In the nucleotide part for each marker separately, a model from the best selection according to AIC implemented in jModelTest 2.1.1 (Darriba et al., 2012) was used. For the binary coded gaps, a Jukes–Cantor model (Jukes & Cantor, 1969) was used. 15,000,000 generations were performed in two runs discarding the first 25% as the burnin fraction and sampling every 1000th tree. As an outgroup Cypripedium subtropicum was used. To estimate node ages, function chronos in package “ape” was used (Paradis, Claude & Strimmer, 2004) applying lambda set to 20. In this function, branch lengths elucidate mean numbers of substitutions and node ages are estimated using semi-parametric method based on penalized likelihood (Sanderson, 2002; Paradis, Claude & Strimmer, 2004). As calibration points splits between Apostasioideae and Cypripedium (80 Mya), and Apostasia and Neuwiedia (43 Mya) were used following published divergence times (Gustafsson, Verola & Antonelli, 2010). To reconstruct ancestral ecological tolerances and predicted niche occupany profiles (PNOs) the Phyloclim package was used (Heibl & Calenge, 2013) which implements the methodology originally developed by Evans et al. (2009). PNO profile takes into account species probability distribution derived from ecological niche modeling and compiles a response to a particular environmental variable for each species. Ancestral ecological tolerances are computed from the phylogenetic tree and PNO using nonparametric approach and ancestral character estimation (Paradis, Claude & Strimmer, 2004; Evans et al., 2009).

Models evaluation

All the created models received high AUC scores of 0.898–0.989. These results are consistent with the outcomes of previous studies which indicated the reliable performance of this method for developing ecological niche models based exclusively on presence-only data (Elith et al., 2006).

Glacial refugia

The ENM analysis indicated eastern Mindanao, Borneo, Sumatra, eastern Sunda, and northern Sahul as areas suitable for the occurrence of Apostasia nuda. The potential refugia of A. odorata are widely distributed along the coast of Sahul and Sunda, as well as between those shelves. The highly suitable niches of Apostasia wallichii were limited to the north-eastern part of Sunda. The created models indicated an even smaller potential coverage area of the glacial refugia of Neuwiedia borneensis. The most suitable niches were limited to numerous small islands north-east of Sahul. Some less appropriate habitats could be found in Borneo and south of this region, as well as along northern Sauhl. The refugia of N. veratrifolia were widely distributed in Sunda, Sumatra, Java, Mindanao, and northern Sahul. Suitable niches of Neuwiedia zollingeri during LGM were located in islands west of Sunda, eastern Mindanao, Sulawesi and smaller islands north-east of Sahul. Figures 1 and 2.

Current potential range

Currently, the most suitable niches of Apostasia nuda are located in Peninsular Malaysia, Borneo, Sumatra, western Mindanao, Sulawesi, Halmahera, north-western New Guinea and New Britain. The potential range of A. odorata is very wide within the studied area. No suitable climatic conditions were indicated by the created models in the areas of higher elevation in continental south-east Asia and the central part of New Guinea. In Australia the appropriate niches are located in Tasmania, in Cape York Peninsula, and along the northern, south-eastern and south-western coasts. Some areas characterized by the climatic conditions suitable for this species are also located in New Zealand. The potential range of A. wallichii includes the coasts of Thailand, Cambodia and Vietnam. It may also occur in Burma, southern Sumatra, Peninsular Malaysia, the Greater Sunda Islands, Borneo, Sulawesi, Muluccas, New Guinea, Cape York Peninsula, the Solomon Islands, Vanuatu, New Caledonia, and Fiji. Suitable niches of N. borneensis are located in Borneo, northern New Guinea, and Halmahera. Areas characterized by less appropriate climatic conditions are also found in Sumatra and the Solomon Islands. The potential range of N. veratrifolia includes Peninsular Malaysia, northern Sumatra, Palawan, Mindanao, Sulawesi, Moluccas, central New Guinea, and the Solomon Islands. Some less suitable niches are located in Borneo, Java, and Vanuatu. Habitats of N. zollingeri are distributed in Laos, Vietnam, Peninsular Malaysia, the Nicobar Islands, Sumatra, Borneo, the Philippines, Sulawesi, Moluccas, New Guinea, the Solomon Islands, Vanuatu, New Caledonia, and Fiji, as well as on the northern coast of Australia, and the Solomon Islands. Figures 3 and 4.

Limiting factors and niches overlap

The distribution of all the studied species is related mostly with the mean diurnal range (Table 2). The second important climatic factor limiting the occurrence of apostasioid orchids is the precipitation in the warmest quarter, which significantly influenced the models of the three Apostasia species and Neuwiedia borneensis. The geographical ranges of Neuwiedia veratrifolia and N. zollingeri are related with isothermality and the annual mean temperature. The latter variable has a somewhat less significant contribution to the models of four other orchids.

The predicted niche occupancy graphs (Fig. 5) showed that while most of the analyzed climatic variables had a similar influence on the distribution of the studied species, some differences are observed in preferences related to precipitation: the precipitation of the driest month (bio14), precipitation seasonality (bio15), and the precipitation of the coldest quarter (bio19). Apostasia species grows in regions with low precipitation values recorded in the driest month, while Neuwiedia representatives prefer rainfall of about 100–130 mm at this time of year. In the aspect of precipitation seasonality, the value of about 25 is preferred by most of the studied species. However, three studied Apostasia representatives may also grow in areas characterized by higher scores of this factor (about 80–85). While in the coldest quarter Apostasia species prefer low amounts of rainfall, Neuwiedia representatives are adapted to higher precipitation. The niches occupied by the studied Apostasia species did not overlap significantly during LGM (Table 3). A. wallichii and A. odorata received the highest scores of D = 0.603 and I = 0.849. In the glacial period, Neuwiedia representatives shared slightly more habitats. The highest overlap (D = 0.600, I = 0.832) was calculated for N. zollingeri and N. veratrifolia (Table 3). While within Apostasia the same pattern of niche commonality is observed in the present time, the habitat overlap within Neuwiedia has changed (Table 4). Currently, the highest niche overlap is observed for N. borneensis and N. veratrifolia (D = 0.577, I = 0.854).

Interestingly, the tests indicated that the potential ranges of several of the studied species overlap more significantly with representatives of other genera than their closest relatives. During LGM N. borneensis and N. veratifolia shared most of their potential range area with A. nuda and N. zollingeri with A. odorata. The same pattern is also observed at the present time. The overlap of the suitable niches of some species has decreased since LGM, e.g., A. odorata and A. nuda, A. nuda and N. zollingei, A. wallichii and A. odorata, N. borneensis and A. odorata, N. veratifolia and A. odorata, N. veratifolia and A. wallichii, N. zollingeri and A. wallichii. In other cases the potential ranges are currently more similar, e.g. A. wallichii and A. nuda, A. nuda and N. borneensis, A. nuda and N. veratifolia, N. zollingeri and A. odorata, N. borneensis–and A. wallichii.

Evolution of climatic tolerance and predicted niche occupancy (PNO)

In this study the intra-specific variation was not addressed as it would require denser sampling throughout the whole region. However, when accepting taxonomical backbone and treating species as being monophyletic, the assumption is made that there is no potential gene flow among taxa or possible diverse adaptations that are present within the intraspecific lineages. This may lead to some bias but currently available phylogenetic studies (Judd, Stern & Cheadle, 1993; Kocyan et al., 2004) suggest that it should be of minimal importance since there is no gene flow between species nor any visual signs of local adaptations. The phylogenetic tree with marked estimated divergence times is presented in Fig. S1.

The evolution of climatic tolerance (Fig. 6) of all variables within the studied taxa indicated that Apostasia odorata developed somewhat broader ecological tolerance for four analyzed variables (bio1, bio5, bio8, bio13), in the case of another four (bio2, bio3, bio4, bio12) two species, A. odorata and Neuwiedia zollingeri, have distinctly broader tolerance. In relation to several climatic factors (bio14, bio18, bio19) all the studied species present similar tolerance. It is noteworthy that A. nuda shows similar preferences in the aspect of many climatic variables (bio2, bio3, bio4, bio12, bio14, bio15 and bio19) to N. borneensis and N. veratifolia. In most cases, the tolerance of the studied species’ ancestors diverged a long time ago into two lineages: the first included Apostasia species and the other Neuwiedia, and among these two groups the diversification of their ecological tolerances began relatively recently. Furthermore, the presented predicted niche occupancy profiles show that most species have very similar climatic requirements in respect to the studied variables, and deviations may be noted only in a few cases (Fig. 5).

Postglacial range changes

The ENM analysis did not indicate any significant changes in the distribution of the suitable niches of most of the studied taxa. The highest congruence is observed in the potential glacial and current ranges of Neuwiedia borneensis and Neuwiedia zollingeri. ENM analysis indicated a minor loss of areas characterized by suitable niches in Apostasia odorata and Neuwiedia veratrifolia. A slight extension of the potential range is observed in Apostasia nuda and A. wallichii. The second species probably occupied niches of lower suitability during the LGM.

Information provided in the atlas of palaeovegetation (Adams & Faure, 1997) indicated that all the studied species occupied areas covered by tropical rainforest and monsoon or dry forest during the LGM. Apostasia odorata and A. wallichii could have also occurred in relatively low, usually deciduous tropical woodlands. Potential glacial refugia of Apostasia nuda and Neuwiedia veratrifolia were also located in savannas. Apostasia odorata could also have survived the glacial period in montane tundra, grassland and semi-arid temperate woodlands or scrublands. The available data on the current distribution of the studied species are consistent with the habitat preferences during the LGM, and no significant shift is observed in any of these orchids. This would partially explain the lack of extensive migration after the glacial period.

Potential vs documented geographical ranges

The potential range of the studied species was compared with the current knowledge on their geographic distribution. While the models overlap mostly in the regions from which those orchids have been reported, we found several incongruences between ENM outcomes and floristic data. Smaller islands of Melanesia and Polynesia were indicated in the analysis as regions suitable for Apostasia odorata, A. wallichii and Neuwiedia zollingeri; however, we did not find any information on their occurrence in those areas. Additionally, we did not confirm the existence of Apostasia nuda or Neuwiedia borneensis in New Guinea, Neuwiedia zollingeri in the Philippines, or Apostasia odorata in Australia. Populations of Apostasia nuda were reported, inter alia, from Bangladesh, India, Cambodia and Vietnam, but ENM analysis did not indicate continental south-east Asia as suitable for this species. This may be a result of the lack of occurrence data from this area included in the input dataset, and it would suggest that the continental populations grow in different climatic conditions from those distributed along the islands and in the Malay Peninsula.

The most probable reason for the inconsistencies between the potential and recognized geographical ranges of those species is the geographical and/or climatic barriers that prevent migration to suitable habitats across the sea (Hosner et al., 2014). On the other hand, we cannot exclude the extinction of some populations in their refugia localized on small islands after deglaciation and sea level rising (Dávalos & Russell, 2012). In this situation, there would be no progenitors in some parts of the potential ranges of the studied species. It should also be emphasized that the floristic data from numerous island groups are relatively obscure, and the populations of Apostasioideae representatives may be found in those regions in the future.

Niche conservatism and evolution of climatic tolerance

The niche conservatism of Apostasioideae is relatively high, and the tolerances of the studied climatic variables of the representatives of this group are rather narrow. The only species which may survive in a broader range of climatic factors is Apostasia odorata. Our studies indicated that the distribution of both Apostasia and Neuwiedia is limited by the same variable—the mean diurnal range. Moreover, within the examined orchids the overlap of their potential ranges is high, and in some cases representatives of Apostasia shares more habitats with Neuwiedia than with other congeners.

The reconstructed evolution of climatic tolerances of the studied species indicated that after the initial divergence of separated lineages of Apostasia and Neuwiedia the differentiation of the climatic niches within both genera began. Separation between the two genera occurred approximately 43-41 Ma ago (Gustafsson, Verola & Antonelli, 2010) while the diversification of currently recognized species started around 21 to 10 Ma ago (Guo et al., 2012). This pattern was also observed in other flowering plants, e.g., palaeotropical genus Pseuduvaria (Annonaceae; Su & Saunders, 2009). In the time of divergence of Apostasioideae from other orchids, about 45 Ma ago (Hall, 2002) Australia began to move rapidly northward. Subduction of the Australian plate resumed beneath Indonesia, causing widespread volcanism at the active margin and producing chains of the Sunda Arc stretched from Sumatra, through Java and the north arm of Sulawesi, and continued into the western Pacific to join the Philippines-Halmahera Arc (Lohman et al., 2011). This island chain could be colonized during several events and as proven by Tänzler et al. (2016) the in situ diversification on different segments of the Sunda Arc occurred. The speciation within studied genera was most probably a consequence of the collision of Sahul and Sunda in the mid-Miocene that resulted also in the emergence of much of the existing land east of Wallace’s Line (Morley, 1998). Moreover, a second collision event initiated in western New Guinea, creating the Central Range orogeny (Cloos et al., 2005; Baldwin, Fitzgerald & Webb, 2012).

The niche conservatism within Orchidaceae remains poorly recognized and previous studies were focused mostly on invasive species (e.g., Kolanowska, 2013; Kolanowska & Konowalik, 2014). The degree to which plants and animals retain their ancestral ecological traits plays essential role in speciation; however, in many orchids this pattern may be misinterpreted and the phylogenetic niche conservatism signal may be disturbed by other factors promoting evolution of new species, especially the role of pollinators. As assumed by Ramírez et al. (2007), pollination in Apostasioideae must be accidental and the speciation within this group was apparently not related to pollinator specificity. For this reason, apostasioid orchids are ideal model for testing niche conservatism using ENM and molecular data. While the revealed pattern cannot be extrapolated into other Orchidaceae, it brings important information on the ecological differentiation in the earliest stage of orchid evolution.