Phylogenetic climatic niche conservatism and evolution of climatic suitability in Neotropical Angraecinae (Vandeae, Orchidaceae) and their closest African relatives

Selection of taxa

Only the species included in the phylogenetic studies of Carlsward et al. (2006) and Micheneau et al. (2008) are included in this study. The localities of the populations of the following species are listed: Angraecum chevalieri Summerh. A. cultriforme Summerh. (=Angraecoides cultriforme (Summerh.) Szlach., Mytnik & Grochocka), A. eichlerianum Kraenzl. (=Eichlerangraecum eichlerianum (Kraenzl.) Szlach., Mytnik & Grochocka), A. erectum Summerh. (=Angraecoides erecta (Summerh.) Szlach., Mytnik & Grochocka), Campylocentrum fasciola (Lindl.) Cogn., C. lansbergii (Rchb. f.) Schltr., C. micranthum (Lindl.) Rolfe, C. pachyrrhizum (Rchb. f.) Rolfe, C. poeppigii (Rchb. f.) Rolfe, C. tyrridion Garay & Dunst., Dendrophylax barrettiae Fawc. & Rendle, D. fawcetti Rolfe (=Polyrrhiza fawcetti (Rolfe) Cogn.), D. filiformis (Griseb.) Benth. ex Fawc. (=Harrisella filiformis (Sw.) Cogn.), D. funalis (Sw.) Fawc., D. lindenii (Lindl.) Benth. ex Rolfe (=Polyradicion lindenii (Lindl.) Garay), D. porrectus (Rchb. f.) Carlsward & Whitten (=Harrisella porrecta (Rchb. f.) Fawc. & Rendle), D. sallei (Rchb. f.) Benth. ex Rolfe (=Polyradicion sallei (Rchb. f.) Garay) and D. varius Urb. Based on morphological studies (Kolanowska, Pérez Escobar & Parra Sánchez, 2012; Kolanowska & Szlachetko, 2013; Szlachetko & Kolanowska, 2013a; Szlachetko & Kolanowska, 2013b) we consider C. sullivanii Fawc. & Rendle as conspecific with C. fasciola as these taxa are morphologically indistinguishable. This is also the conclusion of Bogarín & Pupulin (2010).

List of localities

The occurrence data were obtained from the information recorded with the specimens deposited in herbaria AMES (Orchid Herbarium of Oakes Ames), AMO (Asociación Mexicana de Orquideología), BM (Natural History Museum, London), COL (Universidad Nacional de Colombia), F (Field Museum of Natural History), FLAS (Florida Museum of Natural History), K (Royal Botanic Gardens), MO (Missouri Botanical Garden), NY (New York Botanical Garden), and P (Muséum National d’Histoire Naturelle), as well as from the literature and the original protologues. Herbarium acronyms follow Index Herbariorum (Thiers, 2015). The georeferencing process followed is that used by Hijmans et al. (1999). The geographic coordinates on the herbarium sheet labels were verified. If there were no geographic coordinates on the herbarium label, we used the description of the locality where the plant was collected and assigned coordinates as precisely as possible to this location. The Google Earth (ver. 6.1.0.5001, Google Inc.) application was used to validate all the information gathered. A total of 162 localities that could be precisely located were used in the ENM (5–46 localities per species; Table S1). Dendrophylax filiformis (Griseb.) Benth. ex Fawc. (= Harrisella filiformis (Sw.) Cogn.) and Dendrophylax funalis (Sw.) Fawc. were not included in this analysis because of the lack of a precise locality for them.

Ecological niche modelling and niche similarity

Input data for the ecological niche modelling were 35 bioclimatic variables with a 10 arc minute resolution obtained from the CliMond dataset (Kriticos et al., 2012). To choose an appropriate modelling extent (i.e., the species movement limits described in the M of the BAM diagram; see Barve et al., 2011) terrestrial ecoregions where species occur were selected. Biotic regions serve as a reliable estimate of the area that is accessible to a species and the easiest way to accurately designate the distribution limits of a species (Soberón & Peterson, 2005; Barve et al., 2011). For this purpose, the world map of terrestrial ecoregions was used (Olson et al., 2001). The bioclimatic maps were clipped to include only those regions where it was confirmed that the species occurred. In addition, some island archipelagos and island-like regions were removed from the study area (i.e., Cape Verde archipelago, Ascension Island, Saint Helena, Nile Delta and some smaller features in North Africa north of 19°N). To account for the multicollinearity nature of the initial variables, a principal component analysis (PCA) was performed in R (R Core Team, 2014).

To select an appropriate number of PCA derived maps (i.e., PC axes) the Kaiser-Guttman criterion, which is based on a mean of all eigenvalues was used by selecting only PC axes with eigenvalues larger than this mean. Following this criterion, maps based on the first six principal components were selected (Table S2). The ecological niche modelling was done using a maximum entropy method implemented in Maxent version 3.3.3k (Phillips, Dudík & Schapire, 2004; Phillips, Anderson & Schapire, 2006; Elith et al., 2011), which is based only on species presence. The maximum iterations were set to 10,000 and convergence threshold to 0.00001. Random seeds were used to keep 1,000 bootstrapping runs from using replicate test and training samples. For each run 20% of the data were set aside as test points (Urbina-Cardona & Loyola, 2008). In order to maintain a sufficient sample size for taxa with a small number of occurrences (below 10) duplicate presence records were not removed and presence points were duplicated within a given grid cell. This procedure was only used for taxa that are narrowly distributed insular endemics with a small number of occurrences for which it would be difficult carry out climatic filtering (Varela, Anderson & García-Valdés, 2014). To simplify the interpretation of the probability of a species being present the logistic Maxent output format was chosen.

The model was evaluated using the most common metrics, area under the curve (AUC; Mason & Graham, 2002) and true skill statistic (TSS; Peirce, 1884). AUC was calculated using the Maxent application automatically based on the training localities. Whilst some authors suggest that it may be misleading (Lobo, Jiménez-Valverde & Real, 2008), it seems to be a valid metric for determining the reliability of the fit of the ENM (Warren, Glor & Turelli, 2008). TSS was calculated using maximum training sensitivity plus specificity.

The similarity between the niches occupied by the species studied was measured using Schoener’s D (D; Schoener, 1968) and I statistic (I; Van der Vaart, 1998; Warren, Glor & Turelli, 2008; Warren, Glor & Turelli, 2010) implemented in ENMTools package for R (Warren, 2016), using the methods for calculating environmental distances proposed by Broennimann et al. (2012). Schoener’s D statistic uses direct measures of species density, which in this study are changed to measures of densities of occurrence modelled in environmental space. ‘I’ statistic is based on the modified Hellinger distance that compares two probability distributions. These two metrics range from 0 (no similarity) to 1 (high similarity). To test the importance of distance the niche identity test was used and calculated by the same function that produced the D and I statistics in the ENMTools package for R (Broennimann et al., 2012; Warren, 2016). To visualize the suitable climatic niche of each species, a PCA of the raw climatic conditions occurring at given locations was performed using R (R Core Team, 2014).

Operations on GIS data were carried out on ArcGis 9.3 (ESRI) and R (R Core Team, 2014).

Phylogenetic analysis

To construct the phylogenetic tree, sequences from ITS and trnL-F (Table S3) were aligned using Mafft 6.833b (Katoh & Toh, 2008). Gapcoder (Young & Healy, 2003) was used to code indels. The Alignments were then merged and a Bayesian phylogenetic analysis was performed using MrBayes 3.2.1 (Ronquist et al., 2012). For the nucleotide part, the best model according to the 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 burning in fraction and sampling every 1,000th tree. To estimate node ages function chronos in package “ape” was used as described previously (Kolanowska et al., 2016). Lambda was set to 20 and ages were estimated using a semi-parametric method based on penalized likelihood where branch lengths indicate mean numbers of substitutions (Sanderson, 2002; Paradis, Claude & Strimmer, 2004). As this calibration indicates splits between Campylocentrum and Dendrophylax (max. 12.72 Mya, min. 4.64 Mya), and Angraecum and Polystachya (max. 32.48 Mya, min. 22.2 Mya) we used published divergence times (Givnish et al., 2015; Andriananjamanantsoa et al., 2016).

To reconstruct ancestral climatic suitabilities the Phyloclim package was used (Heibl & Calenge, 2013) which implements the methods originally developed by Evans et al. (2009). Predicted niche occupancy (PNO) was reconstructed and together with the phylogenetic tree were used to infer ancestral climatic suitabilities. PNO integrates species probability distributions (derived using MaxEnt) with respect to climatic variables. Ancestral climatic suitabilities are the PNOs projected onto the phylogenetic tree. They were reconstructed based on the distribution of climatic suitabilities (PNOs) using maximum likelihood and Brownian motion assumption to plot them at each interior node of the tree.

In addition, a Mantel test was used to verify the correlation between the genetic and niche distances (Serra-Varela et al., 2015). This analysis was performed using the Mantel test available in the ADE4 package of R (Dray & Dufour, 2007) and its significance was assessed by performing 9,999 replications. Genetic distances were generated in PAUP* (Swofford, 2002) using the distances predicted by the GTR model and as niche distances Schoener’s D indices were used. GTR distances were chosen based on the best AIC score in the jModelTest 2.1.1 (Darriba et al., 2012).

Age range correlation (ARC, Fitzpatrick & Turelli, 2006) was performed in R using package phyloclim (Heibl & Calenge, 2013). ARC is equivalent to phylogenetic independent contrasts (PIC) and explains how niche similarity (in this case measured using D and I statistics) change over time, which is represented by the nodes of the phylogenetic tree (Fitzpatrick & Turelli, 2006). Monte Carlo resampling with 3,000 replicates was used to assess the statistical significance. This test was used to verify whether any of the observed differences could be explained by phylogeny or whether these are more probably a result of “ecological drift” than “ecological specialization”.

Phylogenetic relationships

The Bayesian tree, based on sequences of both ITS and trnL-F, indicate that Campylocentrum and Dendrophylax form well-supported clades, which originated from African representatives of Angraecum (Fig. 1). These analyses support the results of the previous molecular studies of Carlsward et al. (2006), Micheneau et al. (2008) and Szlachetko et al. (2013).

Evaluation of the modelling of the ecological niche

All projected niche models had high AUC scores ranging from 0.690 to 0.990 (Table 1) indicating the model is good in distinguishing presence data from background data (Phillips, Anderson & Schapire, 2006). The TSS scores showed the same trend but were lower than the AUC scores and ranged from 0.499 to 0.974. The low scores of AUC and TSS are mainly for taxa that are narrow endemics (e.g., those occurring on Caribbean islands). As these scores are close to that of a random prediction they should be treated with caution and may indicate that climate is not a significant factor influencing their distribution (most probably they originated by allopatric speciation, which is not necessarily reinforcing climatic differences).

Distribution of suitable niches and niche similarity

The distribution of the suitable niches of all species studied is presented in Figs. 2–4. This indicates there is a relatively low cover of suitable habitats for Dendrophylax fawcettii and Campylocentrum lansbergii.

The Mantel test did not reveal a correlation between genetic and environmental distance (r =  − 0.14, p = 0.87) supporting the hypothesis that niche evolution (i.e., niche difference) is not related to phylogeny and that phylogenetically distant taxa can occupy similar niches. In addition, the results of this analysis (Fig. 5) did not indicate a correlation between the similarity of the niches occupied by the species studied and their phylogenetic relationships (Fig. 1). Overall, the most similar niches are occupied by Angraecum cultriforme and Dendrophylax porrectus (I = 0.99, D = 0.90), whilst the greatest differences were recorded for Dendrophylax fawcettii and Angraecum erectum (I = 0.10, D = 0.03). The general similarity of the niches occupied by Angraecinae is also visible in the PCA graph (Fig. 6). In this diagram the suitable niches of many clades overlap significantly. It is noteworthy that the niches of all Dendrophylax species fall within those of Angraecum and Campylocentrum. The climatic suitabilities of these particular species also overlap each other and there are no clear boundaries between the genera (taxa) studied in terms of their preferred climatic conditions.

Ancestral climatic suitability

There are no significant differences in the evolutionary history of the taxa studied in their climatic suitabilities and tolerance of species overlap (Fig. 7). Species of both the Neotropical genera included in this study appeared to have slightly wider climatic suitabilities than their closest African relatives. The climatic suitability of Angraecum eichlerianum measured in terms of the climatic factors analyzed has changed very little over time. The niches of the other species of Angraecum studied (A. erectum, A. chevalieri and A. cultriforme) diverged and stabilized much later. The climatic suitability of the Neotropical species evolved in various directions and the climatic niche was apparently not highly conserved within particular genera. Climatic suitability of the African Angraecum cultriforme and A. erectum changed and became more similar to that of Neotropical Angraecinae.

Similar results were obtained from the age range correlation using both D and I for which the slope is negative and intercept is higher than 0.5, which could indicate sympatric speciation (Fitzpatrick & Turelli, 2006; Warren, Glor & Turelli, 2008). These results imply that recently diverged nodes are more similar than more ancestral nodes and recent species occupy similar niches (Fig. 8). ARC was not significant (Table 2), which indicates that climatic niche differentiation did not play a role in the diversification of Angraecinae. The evolution of Angraecinae corresponds to Brownian motion. The species in this group do not exhibit ecological specialization within their climatic suitabilities.

Phylogenetic niche conservatism

Traditionally, natural selection and adaptation result in lineages inhabiting different environments. However, the similarity between the niches occupied by Neotropical angrecoid orchids and their African relatives indicated by this research suggests that they have retained the general climatic suitability of their common progenitor. The recorded shift in climatic suitability of two ancestral African angrecoid orchids could have occurred relatively recently, after the stabilization of niche preferences within New World species of this subtribe. Another possible explanation is that their climatic suitabilities were derived from that of the ancestor of these taxa.

Wiens (2004) suggests that limited adaptation to environmental conditions is a crucial factor in promoting the divergence in the initial origin of lineages and less important in the subsequent divergence of these lineages. The results of this study reveal that climatic niche diversification did not significantly influence the speciation of the Neotropical angrecoid orchids. This is the third study showing that PNC in Orchidaceae is very high (Kolanowska & Szlachetko, 2014; Kolanowska et al., 2016) and that their diversity cannot be explained by niche differentiation. Our results confirm Wiens’s assumption about the low effect of niche modification in speciation.

Evolutionary implications

The divergence of Neotropical angrecoid orchids from their African relatives took place ca. 11.5–16 Mya (Micheneau et al., 2010), in the Miocene, but the diversification within the American representatives of Angraecinae began about 6–7 Mya. The late Miocene marked the start of a paroxysm resulting in the uplift of the northeastern Andes (Hoorn et al., 1995) and a significant increase in sediment rates between ca. 7.9 and 6 Mya (Uba, Strecker & Schmitt, 2007). The Andean sediments reached the Atlantic coast via the Amazon drainage system, and the Amazon River became fully established at about 7 Mya (Hoorn et al., 2010). The rise of the Andes influenced climatic conditions in South America and the range began to constitute a solid migration barrier that limited the dispersal of Campylocentrum. Uplift of the Lesser Antilles arc and the associated platform initiated during Early-Middle Miocene formed a barrier to sedimentation between the two basins by the late Miocene. After the Miocene, most deposition in the Greater Antilles was confined to present coastal and offshore areas (Khudoley & Meyerhoff, 1971). This heterogeneity in topography and environmental conditions probably led to rapid divergence between populations, as it reduced the dispersal and gene flow between adjacent populations. Based on the geological events that occurred at the time of diversification of Neotropical Angraecinae and the result of the ENM analysis we assume that adaptation to various niches was not the main factor determining speciation within the taxa studied.

The climatic-based potential ranges of the Neotropical species studied are considerably larger than the observed distributions of these orchids, indicating that abiotic conditions are not the main factor limiting the occurrence of Campylocentrum and Dendrophylax. Whilst the species of these two genera are characterized by similar climatic suitabilities and they share the same evolutionary and geological history, their known geographical ranges differ significantly. Populations of Campylocentrum occur on both sides of the Andes as well as in Mesoamerica and the West Indies. The distribution of Dendrophylax is restricted to Mesoamerica and the West Indies. Since there is no report of a firm relationship between any of the American angrecoid orchids and a specific phorophyte, we believe that only two factors could promote the diversification of these two genera: pollinator specificity and/or mycorrhizal specialization. Unfortunately, little is known about both of these two aspects of angrecoid orchid biology. Preliminary studies reveal that some species of Campylocentrum form mycorrhizal associations with Ceratobasidiaceae (Otero, Ackerman & Bayman, 2002) and it is possible that the same symbiotic association occurs in Dendophylax lindenii (Chomicki, Bidel & Jay-Allemand, 2014). Undoubtedly, leafless species rely on endomycorrhizal associations, but the diversity of their mycobionts remains unknown.

Pollinator specificity is suggested as a promoter of evolution and speciation in angraecoid orchids by Dressler (1981). Both Capylocentrum and Dendrophylax are entomophilous, although flowers of Dendrophylax are reported to be pollinated by sphignid moths and the pollen of Campylocentrum is also transferred by halictid and meliponini bees (Singer & Cocucci, 1999; Singer, 2003). Whilst species of Campylocentrum grow sympatrically (Pupulin, 1998; Pupulin, 2001; Bogarín & Pupulin, 2010), there are no records of hybrids between any the Neotropical angrecoid orchids. The lack of any significant differences in climatic suitability is also recorded in the African Angraecum, a close relative of Neotropical Angraecinae. The African angrecoids are highly specialized in terms of moth pollination. Most of Angraecum species produce long-spurred, white, nectariferous flowers. These floral characters are consistent with hawkmoth pollination (Darwin, 1862; Grant, 1985; Nilsson et al., 1985; Haber & Frankie, 1989; Micheneau et al., 2008). However, long-spurred species, especially those in Madagascar, are recognized as ancestral in the angraecoid orchid group (Nilsson et al., 1985). The studies on the pollination of species of Angraecum have mainly been on the Malgasy (Nilsson et al., 1985; Arditti et al., 2012) and Mascarene Islands (Micheneau, Fournel & Pailler, 2006; Micheneau et al., 2008). Moreover, the research of Micheneau et al. (2008) indicates that species of section Hadrangis, which have atypical short-spurred, scentless flowers are pollinated by birds and crickets. Unfortunately, the data on pollination are still very scanty and the relationship with pollinators of Angraecum species on Mainland Africa has not been studied. Most likely differences in pollinator composition provides a strong barrier to gene flow.

In conclusion, we found that niche conservatism in Angrecinae is very strong and that the African and Neotropical species in this group have similar climatic suitabilities. Moreover, climatic niche differentiation does not appear to be an important factor in the speciation of Dendrophylax and Campylocentrum. We suggest that pollinator specificity or restricted mycorrhizal associations played a crucial role in the development of the diversity of species in the taxa studied.