The SSP was dated by the following secondary calibrations obtained from Eiserhardt et al. (2011b): core north Andean clade at 11 My and crown of genus at 28 My, both using normal distribution with a standard deviation of 1, a log-normal relaxed clock, and a birth death branching prior.

Dating of the SCPs was performed by penalized likelihood using the function chronos in “ape” R package v. 5.4–1 (Paradis et al., 2004; Paradis, 2013) (lambda smoothing parameter = 2, model = correlated) and node calibration for the crown Aiphanes (age.min = 27 and age.max = 29), and a second calibration age resulting from population parameter estimation in the coalescent. This divergence age was estimated for lindeniana + linearis clades in SNAPP (Bryant et al., 2012) implemented in Beast 2.6 (Bouckaert et al., 2019) on CIPRES Science Gateway version 3.3 (Stamatakis et al., 2008; Miller et al., 2010). For this age estimation, we included all sampled individuals of A. hirsuta, A. lindeniana, A. linearis, and A. concinna in the SCP. We ran 5 chains with a randomly resampled matrix of 800 nuclear SNPs for 10 million generations and then checked them in Tracer v (Rambaut et al., 2018) for chain convergence, sampling efficiency of the priors, likelihood and posterior as well as mutation rates, and theta parameters (these two were left to estimate). Time was obtained from node height of the maximum clade credibility tree after 10% burn-in by conversion using the mutation rate for corn (m = 2.61 × 10^–9) (Gaut et al., 1996).

Species ranges were estimated by extracting species occurrence data from GBIF using the “rgbif” package (version 2.2.0) (Chamberlain et al., 2022) and other data sources (Herbaria: HUA, JAUM, and UTMC). To clean up common errors in the occurrence data, we used QGIS (version 10.2). The coordinates were compared with maps made by experts (Henderson et al., 1995; Borchsenius and Bernal, 1996; Dransfield et al., 2008; Bernal and Borchsenius, 2010; Galeano and Bernal, 2010), as we explain below. Expected distributions were taken from drawn polygons from the existing literature (cited above) and were overlaid to each species coordinates to check for outliers. When we found a dubious point, it was checked by RB (coauthor). If the point could not be verified (or the taxonomy was dubious, meaning the point could have been erroneously assigned to another species), it was discarded.

Georeferencing precision was assessed following Escobar et al. (2016) and verifying with official Base Cartography that the sites described coincided with the coordinates of the biological records. Multicollinearity was evaluated on bioclimatic layers from where the accessible area (m area for species) was cut using an “Extract by mask” algorithm (ArcMap algorithm). Accessible areas (m) were defined by selecting Olson ecoregions that intersected with biological records (Barve et al., 2011). We also performed attribute filtering to eliminate incomplete information, duplicate records, and data from years before 1979. Finally, we stored this information in a geodatabase using the ArcCatalog software (version 10.5).

We used the 19 bioclimatic variables developed by Karger et al. (2017) and available online. We chose CHELSA layers because they incorporate corrections for the effect of wind, the boundary layer, and exposures in mountain valleys (Karger et al., 2017). The layers Bio1-Bio19 incorporate data from the time period of January 1979 to December 2013 and are available at 1 km (30 arcs) resolution. We extracted all climatic layer values for each curated occurrence point using the accessible areas (M) as a mask. This M area was based on the delineation of ecoregions (Olson et al., 2001) and the species range proposed by different authors (Henderson et al., 1995; Galeano and Bernal, 2010; Galeano et al., 2015). We eliminated the correlated predictors using the variance inflation factor (VIF) with the R package usdm v 1.18 (Naimi et al., 2014). Seven bioclimatic variables did not exhibit collinearity (VIF < 10) (Supplementary Table 3). All used coordinates fell within 20° N and 18° S and 60° W and 81° W.

Species distribution models were generated with Maxent maximum entropy algorithm v 3.4.1 (Phillips et al., 2017). This algorithm has been widely used, generating adequate results in exploration of niches and species distribution (Altamiranda-Saavedra et al., 2017; Calixto-Pérez et al., 2018). Maxent has the advantage of providing an evaluation of omission/commission, response curves, and analysis of variable contributions, which is highly useful for understanding the outputs of the model. The background points were generated on the M area to avoid inflation of the AUC. We generated 10 repetitions per species applying the bootstrapping technique and randomly partitioned the species data in each replicate (85% training and 15% validation). All the models were run using default settings (10,000 background points, 500 maximum iterations with a 10^–5 convergence threshold, regularization multiplier of 1, and duplicate occurrence removal). The predictive capacity of each model was assessed using the area under the curve (AUC) value that is generated using the ROC (receiver operating characteristic) technique performed by Maxent. The best model was kept ensuring a test AUC ≥ 0.9 for each species. Finally, the best model of each species was projected to geography and reclassified using a threshold that represents the lowest rate of omission in the data training and testing (<15%). The result was a binary polygon (1/0: presence/absence) (Elith et al., 2011; Phillips et al., 2017). For species with few records and known as endemic to only a few pixels, manual models were produced by reclassifying the digital elevation model Alos Palsar (pixel size: 12.5 m × 12.5 m) to the elevation reported in the coordinates of the few herbarium samples available, corresponding to the narrow distribution of the species. Then, the reclassified raster was cut with watersheds that included all the coordinates.

We tested for niche overlap following the environmental-PCA method proposed by Broennimann et al. (2012) and implemented in the R package ecospat v. 3.1 (Di Cola et al., 2017). A PCA on the selected predictors by the VIF was computed, and the resulting environmental space for the study area was gridded at 1 km² cell resolution. Then, a smooth kernel density function was applied to the occurrence records plotted on the gridded environmental space. The observed niche overlap score for each species pair was estimated with Schoener’s D metric, which ranges from 0 (no overlap) to 1 (complete overlap) (Warren et al., 2008). Then, statistical tests for niche similarity and niche equivalence hypotheses (Warren et al., 2008) were performed. The first test evaluated whether the ecological niche occupied by two lineages were identical and the second assessed whether the ecological niches of two entities were more or less similar than expected by chance. We repeated each test 100 times, returning a null distribution of overlap values to which the observed niche overlap (D) was compared. If the observed D value fell outside of the 95th percentile of expected D values, the null hypothesis of random equivalency/similarity was rejected. Finally, we evaluated if the mean D values per clade was negatively correlated to the age of each clade [Pearson method implemented in cor.test function of the stats package in R Core Team (2020)].

The climatic data for statistical analysis were extracted from the previously curated occurrence points overlaid on the CHELSA layers. The temperature-related (Bio1-Bio11) and precipitation-related variables (Bio12-Bio19) were explored by principal component analyses on temperature and precipitation biplots. Four variables were chosen based on their orthogonality in the PCA, on Sanín et al. (2016), and on our knowledge of factors that could be related to topography and that could determine plant growth; these were: mean annual temperature (Bio1), temperature during the coldest month (Bio6), annual precipitation (Bio12), and precipitation seasonality (Bio15). Although we could have reconstructed the seven variables from VIF, this analysis was exploratory and independent from the SDMs; our aim was to see how they evolved on the phylogenies. These variables were reconstructed as continuous traits on the SCP and SSP phylogenies using the Rphylopars R package v 0.3.2 (Bruggeman et al., 2009) under the Brownian motion default function of continuous trait evolution by Ho and Ané (2014) and using the mean of each species of all values extracted from all occurrences available for each species. We also used these values to conduct a PGLS (phylogenetic generalized least squares) test for phylogenetic signal of two variable sets: the temperature-related and precipitation-related variables using the pgls and pgls.profile functions of the package caper v 1.0.1 (Orme et al., 2012).

The SSP (Figure 2) shows seven clades that agree with the SCPs, although these clades are not all well-supported [posterior probabilities (PP) of 1 for clades acaulis, weberbaueri, and horrida; PP of 0.97 for clades parvifolia and linearis; PP of 0.72 for clade lindeniana; PP of 0.54 for clade simplex]. Support for this seven-clade backbone is strengthened in the SCPs that we discuss below (for support values in all trees, refer to Table 2). Thus, the SSP provides a backbone showing several important relationships: (a) Aiphanes killipii as sister to the rest of the genus, (b) the horrida clade is supported, with three allopatrically distributed species, (c) A. grandis as sister to the mainly Andean species placed by Burret (1932) in subgenus Brachyanthera, and (d) the main pattern inside the Brachyanthera clade. This main pattern shows a well-supported weberbaueri clade, a well-supported acaulis clade, and the remaining predominantly Colombian species organized into four clades (lindeniana, linearis, parvifolia, and simplex clades).

The SCPs provide additional resolution to groups that the Sanger tree did not resolve. The targeted capture method provided a total of 2,557,512 high-quality sequenced bases, with an average coverage of 29.2× per sample. The Maximum Likelihood concatenated RAxML (SCP) and ASTRAL (SCP) topologies reveal similar relationships (Figure 3 and Supplementary Figures 1, 2), with a few exceptions of infraspecific accessions. Support values in these two topologies are significantly higher than in the SSP in the mainly Colombian species (Table 2). Both samplings are not equivalent, as the SSP includes more species and the SCP includes more and complementary accessions for the poorly supported and less studied clades parvifolia, lindeniana, simplex, and linearis; only A. bio and A. cogollo were left unsampled altogether.

There were several important topologic dissimilarities between the trees obtained from different sequencing strategies, with the most relevant being not the species included in clades but the relationship between clades (i.e., the backbone relationships between the clades, not clade definition). Thus, we chose to discuss the clade relationships obtained by the SCPs, as the resulting supports were higher and sequencing more extensive in the genome, and used the SSP for clade definition (i.e., species placement into clades) because of the inclusion of A. verrucosa and A. grandis, as well as A. graminifolia Galeano & R. Bernal, A. weberbaueri Burret, and A. spicata Borchs. & R. Bernal in the sampling. An important inconsistency was related to two of the species forming monospecific clades, A. pilaris and A. macroloba. A. pilaris was sampled in both approaches and was recovered as sister to the parvifolia complex in the SSP but not in the SCP, where it was recovered as sister to everything but A. killipii Burret, the horrida and weberbaueri clades. In the SCP, A. macroloba was placed as sister to both the lindeniana and the linearis clades, whereas in the SSP it was placed as sister only to the linearis clade. Despite these differences, the strategies were complementary by providing species assignment into the clades (SSP) plus relationships between the clades (SCP).

Using the two different calibration schemes for the SSP and SCP resulted in similar divergence age estimates, concentrating in the Miocene-Pliocene. The SNAPP analyses resulted in an estimated divergence time between the lindeniana and linearis clades of 2-1 million years before the present, and tree inconsistencies regarding the relationship among the four included species (Supplementary Table 4).

The first branching A. killipii has a narrow distribution in the Eastern Cordillera of Colombia. Further splits in the phylogeny produce the widely distributed horrida clade, spanning the Tropical Andes and the Caribbean, a weberbaueri clade from mostly Peru and Ecuador, the Colombia and Ecuador Pacific acaulis clade, and then the four most nested and predominantly Colombian Andean linearis, lindeniana, parvifolia, and simplex clades.

The species of Aiphanes present a very small overlap of distribution ranges, with most species in clades showing a less than 10% overlap. Average range overlap was highest in the deltoidea and linearis clades (38 and 31%, respectively), followed by the acaulis and lindeniana clades (18 and 16%, respectively), and with the other clades showing a less than 10% average overlap. Most of the overlapped geographical ranges include a species that is microendemic (known from one or a few close localities) embedded in a range of a widely distributed species (i.e., 99% of A. spicata’s range is within A. weberbaueri’s but only 1% vice versa). Overlap of 31–92% did occur between widely distributed Andean sister species pairs: A. concinna/A. lindeniana, and A. hirsuta/A. linearis. The geographical distributions are available in Supplementary Figure 3 (distribution point maps) and Supplementary Figure 4 (SDMs used to estimate geographical overlap); the table with all the paired overlap estimates is available in Supplementary Table 5.

The earliest diverging members of the genus (except A. minima) occur today in mid-elevation mesic environments (A. killipii, A. horrida, and A. eggersii). Diversification in the genus appears to have occurred recurrently in different temperatures (Figure 6A), and to have been related to the early conquest of areas with higher precipitation and lower precipitation seasonality, and the conservation of this trait (Figures 4B,C, 6B,C and Supplementary Figure 5). These precipitation-related adaptations were conserved both at higher elevations in the northern Andes (linearis and lindeniana clades) and at lower elevations in the cis- and trans-Andean lowlands (acaulis clade in Choco and A. ulei and A. deltoidea in the Amazon). This is shown by the high phylogenetic signal of the joint precipitation-related variables but not of the temperature-related variables reflected in the lower values of Blomberg’s K and Pagel’s Lambda for the former variable set (Table 3 and refer to Supplementary Figure 7).

Both the geographic distribution of species and clades and the climatic evolution of the genus highlight its close relationship with Andean mountain building from the origin of the genus, as shown by the character state reconstruction indicating mid-elevation temperatures from the ancestor of the genus and subsequent nodes (Figure 4). Also, only A. hirsuta reaches Panama and Costa Rica, and A. minima Burret the West Indies; the rest of the species are on the Andean Mountains or surrounding lowlands. The mean annual temperature of the nodes in Aiphanes is reconstructed as mild, which can be understood as mid-elevations at 1,000–2,000 m. The MAT during the late Miocene to Pliocene, when Aiphanes diversified (Eiserhardt et al., 2011b; this study), was probably warmer than today (Zachos et al., 2001; Hansen et al., 2013); thus we expect mild temperatures to be related to paleoelevation. Furthermore, as expected by increasing mountain building, high elevation species belong to the more nested clades of lindeniana and linearis, as well as the species A. spicata and A. pilaris.

The progression from mid to higher elevations probably follows the process of mountain building and points at mid-elevation areas in the Colombian Andes (possibly the Eastern Cordillera where A. killipii is found today and is near two other microendemic species from different clades, A. graminifolia and A. suaita) already available during the Miocene (Anderson et al., 2015). Despite the absence of paleoelevation constraints for the Central and Western Cordilleras, thermochronological and geomorphological considerations suggest that these Cordilleras may have also experienced exhumation and uplift since the Miocene (León et al., 2018; Noriega-Londoño et al., 2020). Nevertheless, the concentrated origin of several early-diverging species in the Eastern Cordillera suggests that this area was either uplifted or connected earlier in time. The high-elevation A. linearis from the Western and Eastern Cordilleras is nested in the phylogeny and is among the last to form. Diversification in the Colombian cordilleras and the Choco (the acaulis, parvifolia, lindeniana, linearis, and simplex clades), as well as tendency to gradually occupy areas with higher precipitation and less seasonality and a wide variety of temperatures, hints at how the evolution of Aiphanes and the North Andean chains are intertwined.

In simple terms, the tracking of these palms of continuously moist environments allowed them to diversify, occupying different temperatures available at the places where they became geographically isolated by mountain building (Figure 7; the elevation ranges for all species are available in Supplementary Figure 8). High level of precipitation and low precipitation seasonality have been the backdrop for species persistence in the dynamic Neogene evolution of Northwestern South America (Hoorn et al., 2010), and are conserved traits in the phylogeny. On the other hand, the repeated temperature differentiation in the different clades is the result of the evolving topography, isolating species in valleys and mountains at different elevations.

Geographical divergence is important in different Andean birds (Hazzi et al., 2018) and plants (Vargas, 2018; Vargas et al., 2020) and is reflected in the low mean geographical overlap between species and their relatively small ranges. Also, it can, in some cases, be concomitant with niche adaptation (Benham and Witt, 2016). Indeed, the temperature differentiation here follows a spreading pattern shown by the parallel laddered occupation of different temperatures of same-clade species (Figure 4, refer to Bio1 occupation for clades horrida in light blue, acaulis in purple, parvifolia in blue, and lindeniana in green). Climatic space differentiation was also shown for an Andean endemic palm genus, Ceroxylon (Sanín et al., 2016), and Andean vertebrates (Patterson et al., 1998; Cadena et al., 2011). In the case of Aiphanes, topography plays a dual role by isolating populations and then truncating their potential occupation along the temperature gradient (Benham and Witt, 2016).

It is widely acknowledged that realized distributions underestimate species’ physiological tolerance (Soberon and Peterson, 2005; Feeley and Silman, 2010), and that physical barriers to dispersal might amplify this effect. The resulting physically constrained extension of species’ natural ranges enforces niche truncation (Bush et al., 2018). In Aiphanes, the rarely overlapping occupied temperature niches between same-clade species might be the effect of dispersal limitations affecting climatic niche realization because of, in particular, topographic barriers. Thus, our estimation of potential niche occupation, which relies on occurrences, might be ultimately biased by how well the species has succeeded in “sampling” the topographically fragmented landscape. Conversely, it is clear that the occupation of precipitation-related climatic space is much more phylogenetically constrained, as it is conserved despite geographical circumscription.

The overall pattern appears as diversification by geographical isolation in moist areas with different temperatures, which has made Aiphanes species-rich on the regional scale but locally rare (Borchsenius and Bernal, 1996) and locally species-poor. Ultimately, this recoils to the genus’ aforementioned paucity. Either the relative scarcity of pollinators or fluctuations in their specificity or other biotic factors, such as poor ability to recolonize areas due to poor seedling recruitment (Svenning, 1998), may have limited the population expansion of most species of Aiphanes. The reproductive biology of this genus has only been studied in A. erinacea and in A. chiribogensis Borchs. & Balslev, where hoverflies and gnats/midges, respectively, were reported as the main pollinators in Ecuador (Borchsenius, 1993). Our field observation is that infructescences and ripe fruits are seldom found in wild populations, but this has not been systematically assessed.

This study could be expanded in different ways. First, we lack natural history knowledge of the physiological tolerance of compound variables like energy-water balance in many tropical plant groups. Second, we lack knowledge of actual dispersal abilities mediated by biotic interactions in tropical taxa. Third, we should further understand geographical isolation: is it slope, topographic complexity, or the temporal and spatial extent of a physical barrier that is related to biodiversity patterns determined by isolation? A limitation of this study stems from these knowledge gaps.

As seen in the previous sections, Aiphanes provides an example of a diversification pattern by the occupation of the elevation gradient following high precipitation with low seasonality. Isolation into different corners of Andean relief leads to temperature differentiation the same way repeatedly across the different clades, and to high abundance of narrowly distributed species. The phylogenetic conservation of precipitation-related adaptations is coupled with possibly fundamentally larger tolerance to different temperatures. However, wider tolerance to different temperatures is not reflected in species occupation of climatic space (i.e., their realized climatic niches), because it is truncated by physical occupation of mountains and valleys from where the species are unable to disperse. This pattern might be common in Andean plant groups but understudied because of lack of taxonomic understanding, collections, sequencing, and validated distribution data for plant groups that include a high proportion of rare and narrowly distributed species.

In Aiphanes, future research should aim at discovering how this often-unnoticeable genus maintains populations and effective reproduction in time despite its elusive presence in mountain forests. This, combined with the expected discovery of additional microendemic species, will help to understand the evolution of this fascinating Andean gem.