We retrieved biogeographic information, i.e., the geographic range of ancestral species and the potential sequence of biogeographic events that led to the current distribution in Camptoloma and Canarina, from our own published studies: Culshaw et al. (2021) and Mairal et al. (2015), respectively. These studies used molecular phylogenies based on several plastid markers and multi-copy nuclear ITS to reconstruct phylogenetic relationships among the disjunct species and their closely related outgroups. Lineage divergence times were inferred using Bayesian relaxed molecular clocks. Biogeographic inference was performed using two parametric approaches that are based on continuous-time Markov-chain (CTMC) models of range evolution with states representing user-defined, discrete areas of distribution (Sanmartín, 2020). Mairal et al. (2015) employed the Bayesian Island Biogeography (BIB) model (Sanmartín et al., 2008), which assumes single-area ancestors and focuses on the sequence of migration events between areas. Culshaw et al. (2021) used the Dispersal Extinction Cladogenesis (DEC) model (Ree and Smith, 2008), which permits widespread ancestors and focuses on scenarios of range division at nodes. CTMC models were implemented in a Bayesian framework using software BEAST (Lemey et al., 2009; Mairal et al., 2015) and RevBayes (Landis et al., 2018; Culshaw et al., 2021), allowing inference of marginal posterior probabilities (mean and 95% high-posterior credibility intervals) for the CTMC parameters and ancestral ranges. In both Camptoloma and Canarina, the geographic range of the stem ancestor, i.e., the divergence between each genus and its sister-group, was inferred with high uncertainty. This is linked to long internodal branches subtending the stem- and crown-nodes in the two genera (Mairal et al., 2015; Culshaw et al., 2021). More details on these analyses can be found in the original studies.

ENMs for genera Camptoloma and Canarina were generated de novo for this study by following the below steps:

Figure 3A shows the continuous post-truncation ENM projections for Camptoloma and Canarina built upon the present-day bioclimatic, full and prec-temp datasets (the latter are presented as the “yearly” monthly averaged values). Figure 3B shows the corresponding post-truncation ENM projections under the truncation threshold value selected under the MinROCdist criterion, which oscillates between a value of 0.93 and 0.99. Supplementary Table 3 reports the performance statistics (error, accuracy, and precision, etc.) for each model under the selected threshold value. All models performed well for all statistics except for precision, which was consistently low across models. In our study, precision assumes that each genus that we have recorded is present in all possible locations with a suitable habitat, which is clearly not the case. Our presence datasets for Camptoloma and Canarina are incomplete because, as mentioned in the introduction, data collection is difficult, and the genera that we are modeling have a limited geographic distribution. Supplementary Table 4 presents the statistics for the strength of association between the full-ENM and the prec-temp and bioclimatic-ENM present-day projections. In the case of the pre-truncation threshold projections, the full and prec-temp datasets showed a “strong association” (according to the descriptor in Dancey and Reidy, 2007) in Camptoloma (Spearman’s rho = 0.45), and a very strong association in Canarina (Spearman’s rho = 0.93). The pre-truncation threshold projections based on the present-day full and bioclimatic datasets showed a weak association in Camptoloma (Spearman’s rho = 0.21), but strong in Canarina (Spearman’s rho = 0.78). The Simpson’s Similarity Index showed a strong to very strong association between the post-truncation threshold present-day projections (Supplementary Table 4) based on the full dataset and the prec-temp and bioclimatic datasets, suggesting that when an appropriate truncation threshold is used, even projections with weak associations become strong. In all, these results support our contention that the monthly mean temperature and precipitation variables included in the paleoclimate datasets hold enough information to make reliable predictions about habitat suitability, and thus can be used for ENM hindcasting over the paleo and preindustrial time slices (see below).

Figure 4 shows the hindcasted ENM projections, both pre-truncation (continuous) and post-truncation (presence/absence), based on the prec-temp dataset in the case of Camptoloma, and on the bioclimatic dataset in the case of Canarina and for each of the four time slices. Supplementary Figure 1 shows the post-truncation ENM projections for different threshold values in Camptoloma (Supplementary Figure 1A) and Canarina (Supplementary Figure 1B) using the same datasets as above. For the paleoclimate time slices I-III, a truncation threshold value of 0.12 was selected for Camptoloma and a value of 0.21 for Canarina. For the Preindustrial world time slice IV, a truncation threshold value of 0.95 was selected for Camptoloma and of 0.99 for Canarina. In both genera, the hindcasted ENM projections show a much larger extension of suitable habitat in the past compared to the present-day distribution. Furthermore, the paleoclimate projections predicted the existence of suitable habitat corridors connecting the currently disjunct species ranges during the Miocene and Pliocene: these involve central and northern Africa during the LMC for Canarina, and the LMC and MPWE for Camptoloma (Figure 4).

Supplementary Figure 2 shows the post-truncation hindcasted ENM projections across different threshold values based on the alternative climatic dataset for each genus, that is, ENM projections for the bioclimatic dataset in Camptoloma (Supplementary Figure 2A) and ENM predictions based on the yearly prec-temp dataset in Canarina (Supplementary Figure 2B). Comparisons indicated that for ENM projections based on these particular datasets and genera, it was not possible to find an optimal threshold value across the three paleoclimate time slices. Therefore, the bioclimatic and prec-temp ENMs for Camptoloma and Canarina, respectively, were excluded from further analysis. Finally, Figure 5 demonstrates the need to use a different truncation threshold value for different temporal scales: we compared post-truncation ENM projections for the three paleoclimate time slices I-III that were generated using either the threshold value for the Preindustrial time slice IV, calculated with the minimized distance ROC curve, versus the post-truncation ENM projections for the three paleoclimate time slices I-III that were generated using an independent threshold value calculated with the biogeographically informed truncation method described in Figure 2. Results, shown here for Camptoloma prec-temp and Canarina bioclimatic ENMs, respectively, indicate that using the MinROCdist method is too restrictive and fails to recover any spatial distribution pattern in the ENM hindcasts of the paleoclimate time slices I-III that fits the biogeographic inference results.

Supplementary Figure 3 and Supplementary Table 4 show the results of the comparison between the pre- (Supplementary Figure 3A) and post-truncation (Supplementary Figure 3B) ENMs generated with the non-inverted versus the inverted prec-temp dataset for Camptoloma (we did not investigate hemisphere effects in Canarina, since, as explained above, the prec-temp ENMs for this genus were excluded from further analyses). The truncation threshold value for the non-inverted prec-temp ENM projections was 0.97 for the preindustrial time slice and 0.1 for the paleoclimate time slices. A very strong association between the non-inverted and inverted prec-temp ENM projections was found for each of the four time slices [pre-truncation threshold: Spearman’s rho range = [0.72, 0.83]; post-truncation threshold: Simpson’s Similarity Index range = (0.78, 0.88)]. Moreover, the non-inverted and inverted post-truncation threshold ENM maps showed a similar history of range shift events, with climatic corridors in the LMC and MPWE time slices connecting the currently disjunct species distributions. This suggests that the hemisphere effect did not bias our yearly monthly averaged ENM projections in Camptoloma, and therefore had no effect in our predictions.

Figure 6iii reconstructs the history of geographic range shifts in Camptoloma and Canarina using information provided by the biogeographic inference analysis (Figure 6i) and the hindcasted ENM projections (Figure 6iii). The inferred ancestral ranges in Camptoloma (Figure 6Ai) matched the climatically suitable areas depicted in the post-truncation ENM projections based on the prec-temp dataset (Figure 6Aii). In Canarina, the hindcasted ENMs based on the bioclimatic dataset (Figure 6Bii) also showed a good match to the biogeographic inferences (Figure 6Bi).

In the original studies exploring Canarina and Camptoloma biogeographic history (Mairal et al., 2015; Culshaw et al., 2021), the currently disjunct distribution of the extant species was explained by ecological vicariance and climate-driven extinction linked to Neogene aridification cycles, in agreement with other Rand Flora lineages (Pokorny et al., 2015; Rincón-Barrado et al., 2021). However, the small phylogenies, with only three species each, and the deep evolutionary divergences, introduced large uncertainty in ancestral range inference; for example, the Western Asian origin of Canarina (Mairal et al., 2015) or the southern African stem-ancestor of Camptoloma (Culshaw et al., 2021), were associated with marginal probabilities lower than 0.1.

In our study, the introduction of additional information provided by hindcasted ENM projections was able to clarify some of these uncertain biogeographic inferences, in particular along the long internodal branches. For example, ENM projections support a restricted distribution of stem-Camptoloma, in southwestern Africa in the Early-Mid Miocene (21 Ma), while the crown-ancestor is reconstructed as occupying eastern Africa-southern Arabia and Macaronesia already in the Early Pliocene (4.6 Ma, Figure 6Aii). ENM projections (Figure 6Aii) support also the hypothesis of a northeastward dispersal of the ancestor of Camptoloma during the 15-million year branch length separating the stem and crown nodes (Figure 6Aiii, Culshaw et al., 2021): the ENM predictions for time slices I and II showed the presence of suitable habitat corridors over central-eastern Africa during the MMCO (17-14 Ma) and the LMC (12-5 Ma) time slices (Figure 6Aiii). Similarly, the biogeographic inference of a widespread distribution across northern Africa for the most recent common ancestor of Camptoloma canariense and C. lyperiiflorum (3.4 Ma, Figure 6Ai), Culshaw et al., 2021) is supported by the ENM projection for time slice III, which shows climatic corridors with suitable habitat during the MPWE (Figure 6Aii); these disappear in the preindustrial projection (Figure 6Aiii).

Mairal et al. (2015) inferred a stem-ancestor of Canarina distributed in central Asia around the Mid Miocene (13.9 Ma, Figure 6Bi), which migrated to eastern Africa during the Mid-Late Miocene (11-8 Ma). This scenario is supported by our hindcasted ENM projections (Figure 6Bii), which show suitable habitat in western Asia and central-eastern Africa during the MMCO interval. The ca. 8 million-year branch separating the disjunct sister species C. canariensis and C. eminii was explained in Mairal et al. (2015) by short-range dispersal over climatic corridors across the Sahara during the Late Miocene, followed by climatic extinction in the Pliocene (Figure 6Bi). The ENM projections for time slices II and III support these inferences, depicting corridors of suitable habitat across the Sahara during the wetter and colder climates of the LMC (Figure 6Bii), which disappear in the projections of the drier MPWE (Figure 6Bi).

A known difficulty in building ENMs for rare endemics is the selection of the optimal truncation threshold value (Diniz-Filho et al., 2010). The selection of a threshold value often relies upon model fit tests, such as the Kappa statistic (Monserud and Leemans, 1992; Fielding and Bell, 1997), the true skill statistic (TSS, Allouche et al., 2006), or the Receiver Operator Characteristic (ROC). All of these approaches rely on using a training dataset, subsetted from the original data, and evaluating model fit through examination of true and false predictions of species presences and absences (Liu et al., 2005; Jiménez-Valverde and Lobo, 2007). When the georeferenced data is complete and unbiased, and the ENMs are hindcasted upon shallow timescales, these truncation threshold criteria are reliable. However, when occurrence records are incomplete or biased toward certain regions, or when ENMs are hindcasted upon deep timescales, as in our study, these methods can lead to biased habitat suitability projections, as demonstrated in Figure 5.

Here, we propose an alternative approach that uses ancestral nodal ranges inferred through biogeographic inference as an independent source of information for selecting the truncation threshold value of ENM projections (Figure 2). Our approach is particularly suited for ENMs hindcasted into the distant past, where traditional criteria such as the MinROCdist result in very restrictive and conservative truncation threshold values and fail to recover any spatial pattern (Figure 5). Our study, thus, supports the importance of using an appropriate truncation threshold criterion when using absence/presence ENM models (Nenzén and Araújo, 2011), especially when comparing paleoclimate vs. predindustrial hindcasted projections.

There has been recent discussion on the importance of including fine-scale temporal resolution in the environmental data used for building ENMs (Kearney et al., 2012; Gardner et al., 2020; Perez-Navarro et al., 2020). The average-year bioclimatic variables obtained from WorldClim might be too coarse in temporal resolution to capture the climatic niche of short-lived species or species with large geographic ranges and contrasting seasonal regimes, for example, those distributed north and south of the equator (Kearney et al., 2012; Montalto et al., 2014). The importance of considering temporal resolution in climatic data has been emphasized also for paleoclimate data. Even if WorldClim bioclimatic variables correlate with physiological predictors, these relationships could break down when results are extrapolated into other time frames (Elith and Leathwick, 2009). Bova et al. (2021) recently showed that the mismatch between paleoclimate simulations and temperature reconstructions from geological records disappear when seasonal variation over time is employed instead of long-term averages.

In the case of Camptoloma, ENM projections built on the yearly averaged, monthly variables of the prec-temp dataset showed a strong association to those built on the full climatic dataset (Supplementary Table 4), suggesting that monthly mean temperature and precipitation contained enough information to infer ancestral tolerances when projected into the past; in contrast, ENMs built on the annual variables of the bioclimate dataset showed a weak association and we could not find an optimal threshold value for hindcasted projections (Supplementary Table 4 and Supplementary Figure 2). For Canarina, the opposite pattern was found: the ENM projections built on the annual bioclimatic variables showed a strong association with the full dataset, while the monthly averaged prec-temp projections proved to be less informative (Supplementary Table 4 and Supplementary Figure 2). These differences between the two genera may be explained by their different “rarity” character. Camptoloma fits Rabinowitz’s (1981) definition of a “narrow rare endemic”: a rare species with small populations restricted to climatic microrefugia. Camptoloma rotundifolium occurs in moist escarpments in the Namibian Brandberg Mountains, whereas C. canariense grows in shaded vertical crevices on the coastal slopes and cliffs of Gran Canaria and C. lyperiiflorum inhabits a special wetter microclimate in Socotra Island (Culshaw et al., 2021). In contrast, Canarina fits the definition of a “widespread rare endemic” (Rabinowitz, 1981): Canarina canariensis occurs in five of the seven Canary Islands, while the geographic range of C. eminii extends from Ethiopia in the north to Malawi in the south (Mairal et al., 2018). Moreover, Canarina species can spend unfavorable seasons throughout the year in a tuber/root form, which could buffer the influence of climate extreme events on species’ survival, growth and reproductive capacity (Hedberg, 1961). It might be that the use of fine-temporal resolution environmental data, i.e., the yearly averaged monthly values from the prec-temp dataset, is important for modeling the climatic preferences of a species with restricted habitats and little capacity for behavioral buffering, such as Camptoloma, whereas the annual bioclimatic variables from WorldClim work well for the more widespread Canarina. Our study thus supports the hypothesis that temporal resolution of environmental data when building ENMs should be at a scale adequate to understand physiological responses to climate change in the organism under study (Gardner et al., 2020), or in other words, at a scale over which biological responses directly influence their predicted presences in ENM models (Montalto et al., 2014).