Broadening Predictive Understanding of Species’ Range Responses to Climate Change: The Case of Aloidendron dichotomum

Brodie, Lara P.; Grey, Kerry-Anne; Bishop, Jacqueline M.; Midgley, Guy F.

doi:10.3389/fevo.2021.715702

ORIGINAL RESEARCH article

Front. Ecol. Evol., 04 October 2021

Sec. Biogeography and Macroecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.715702

This article is part of the Research Topic Predicting and Managing Climate-Driven Range Shifts in Plants View all 12 articles

Broadening Predictive Understanding of Species’ Range Responses to Climate Change: The Case of Aloidendron dichotomum

$\r\nLara P. Brodie$ Lara P. Brodie¹

Kerry-Anne Grey^1*

Jacqueline M. Bishop²

Guy F. Midgley¹

¹Global Change Biology Group, Department of Botany and Zoology, Stellenbosch University, Stellenbosch, South Africa
²Institute for Communities and Wildlife in Africa, Department of Biological Sciences, University of Cape Town, Cape Town, South Africa

Concerns have been raised about attribution of species range shifts to anthropogenic climate change. Species paleo-range projections are emerging as a means to broaden understanding of range shifts and could be applied to assist in attribution. Apparent recent range contraction in the Quiver Tree (Aloidendron dichotomum (Masson) Klopper and Gideon F.Sm) has been attributed to anthropogenic climate change, but this has been challenged. We simulated the paleo- and future geographic range of A. dichotomum under changing climate using species distribution models (SDMs) to provide a broader perspective on its range dynamics. Ensemble modelling of the Last Glacial Maximum (LGM), mid-Holocene, current, and projected 2070 time periods simulates a paleo-historical poleward expansion of suitable bioclimatic space for this species under natural climate change post-LGM, and projects an eastward shift towards 2070. During the LGM, suitable bioclimatic space for A. dichotomum was simulated to be restricted to the equatorward part of its current range. During the Pleistocene/mid-Holocene climate transition period, the species’ range is predicted to have expanded significantly polewards at an average rate of 0.4 km per decade, assuming constant tracking of its optimal climatic niche. By 2070, suitable bioclimatic space is projected to expand further eastward into the summer rainfall region of South Africa, and contract in its equatorward reaches. Simulated post-LGM shifts roughly match expectations based on preliminary phylogenetic information, further supporting the attribution of current population declines to anthropogenic climate change drivers. Equatorward populations are required to migrate south-eastwards at a rate roughly 15 times faster than that calculated for the LGM/mid-Holocene climate transition period to avoid local extirpation. A preliminary analysis of range-wide genetic variation reveals a cline of variation, with generally higher levels in the central and more northerly part of the species distribution, as expected from the proposed paleo-range of the species. A more detailed analysis of the species’ phylogeographic history could be used to test the proposed paleo-range dynamics presented here, and if confirmed, would provide strong support for the use of this species as an indicator of anthropogenic climate change and a powerful case study for testing the implementation of conservation actions.

Introduction

African terrestrial surface air temperatures have warmed by roughly double the rate (Engelbrecht and Engelbrecht, 2015) of the global surface air temperature rise of between 0.8 and 1.2°C since the preindustrial era (IPCC, 2014). There is high certainty that this warming will continue into the future (IPCC, 2018) with increases in Africa projected to exceed 5°C century^–1 in the dry subtropical regions (Engelbrecht and Engelbrecht, 2015). A large body of research and observation work shows that biological systems and species are already responding to anthropogenic climate change (Parmesan, 2006; IPCC, 2014), but the study locations have a strong northern hemisphere bias. A common prediction for plants is that they will shift their geographic ranges in response to warming, either poleward or upward in elevation (Parmesan and Yohe, 2003; Parmesan, 2006; Midgley and Thuiller, 2007; Chen et al., 2011; Midgley and Thuiller, 2011; Manes et al., 2021), with many plant species already having undergone such changes (Chen et al., 2011; Feeley and Rehm, 2012; Telwala et al., 2013).

Range shifts are a well-known response of species to changes in regional and global climates in the distant geological past (Hewitt, 2004; Urrego et al., 2010; IPCC, 2014). Nonetheless, certain plant functional types have shown the early signs of adverse impacts as a result of recently accelerated anthropogenic climate change (Allen et al., 2015). Among these, trees are arguably the most vulnerable, with mortality potentially resulting in rapid and long term alterations of ecosystem types and associated ecological processes (Breshears et al., 2005), and with long generation times slowing the establishment rate of new populations (Pitelka, 1997). Globally, trees are critical for maintaining the climate regulating service of carbon sequestration (Canadell and Raupach, 2008). They typically also play keystone roles in the structuring and functioning of many habitats by providing habitat, shelter, food and mutualistic relationships with many animal, plant, and microbial species (Aitken et al., 2008). The responses of long-lived trees to past changing climates are therefore extremely important for ecological understanding of the systemic impacts of future projected climate change, and the accurate attribution of climate-driven impacts on these long-lived sentinels can support the development of appropriate conservation responses.

Over the past two decades, the Quiver Tree (Aloidendron dichotomum (Masson), formerly Aloe dichotoma), an iconic long-lived stem succulent of the arid western regions of southern Africa (Figure 1), has undergone substantial demographic change (Foden et al., 2007; Hoffman et al., 2010; van der Merwe and Geldenhuys, 2017). Population surveys have shown enhanced mortality of established individuals in the warmer parts of its range (central and equatorward populations), and enhanced population recruitment in the cooler parts (generally the more poleward populations) (Foden et al., 2007; Midgley et al., 2009; Hoffman et al., 2010; van der Merwe and Geldenhuys, 2017). Repeated observations of A. dichotomum population status over several decades, analysis of demographic patterns in relation to expected climate drivers, and consideration of alternative explanations (i.e., protocols meeting the requirements of attribution of range changes to anthropogenic climate change; see Taheri et al., 2021) strongly supports attribution of current population dynamics of A. dichotomum to anthropogenic climate drivers (see Foden, 2002; Foden et al., 2007; van der Merwe and Geldenhuys, 2017). Coupled with projected south and south-easterly range expansion in cooler parts of its range (Foden, 2002; Guo et al., 2016), these observations may suggest that a range shift is underway, although this interpretation has been contested (Jack et al., 2016).

FIGURE 1

Figure 1. Geographic range of Aloidendron dichotomum. Presence localities obtained from Foden et al. (2007) Diversity and Distributions, Foden (pers. obs.) and the Global Biodiversity and Information Facility (GBIF). Mean annual temperature (mean of all monthly temperatures) was obtained from WorldClim (Hijmans et al., 2005) at 2.5-min resolution.

Multiple lines of evidence indicate that A. dichotomum is sensitive to changes in climate (Foden et al., 2007; Hoffman et al., 2010; Guo et al., 2016; van der Merwe and Geldenhuys, 2017). Observations indicate that many populations have not had significant recruitment in the last 50 years, and at sites where limited recruitment has occurred, the timing of seedling establishment correlates with timing of rare favourable climatic windows (Foden, 2002; Hoffman et al., 2010). Further, Guo et al. (2016) the geographic range of A. dichotomum is constrained by temperature-based variables, indicating the potential vulnerability of this species to changing climate conditions, while the range of the species appears insensitive to rainfall seasonality, with positive population growth rates in both the winter-rainfall and summer-rainfall zones of south-western Africa (Foden et al., 2007; van der Merwe and Geldenhuys, 2017).

Further evidence for the species sensitivity to climate change can be gained through insight into the role of glacial-interglacial climate changes and the corresponding range shifts. Past climatic conditions are known to have affected speciation, extinction, and migration dynamics. Quaternary glacial and interglacial cycles have been shown to have caused drastic shifts in species distributions (Graham et al., 1996; Hewitt, 1996; Riddle, 1998; Avise, 2000; Bellard et al., 2012; Lyam et al., 2020; Łabiszak et al., 2021). Using both molecular and paleo-climatic data, Łabiszak et al. (2021) propose that the European peat bog pine (Pinus uliginosa) underwent a significant range contraction during the Last Glacial Maximum (LGM), reflected by a substantial reduction in genetic diversity during this period. Similar trends have been inferred for sub-Saharan savanna trees, with Senegalia Senegal proposed to have undergone a range reduction during the LGM and subsequently expanding its range out of climatic refugia through the mid-Holocene, to its current wider distribution (Lyam et al., 2020).

In many cases, geographic processes have led to distinct signals of past distributional changes through the patterning of genetic differentiation across populations, subspecies, and species (Comes and Kadereit, 1998; Hewitt, 2004; Habel et al., 2005; Łabiszak et al., 2021). In general, however, southern Africa has only sparse paleo-environmental archives to provide a basis for understanding how climates and vegetation have interacted. Development of Global Climate Models (GCMs), which enable reconstruction of past environmental conditions and prediction of future conditions, provide the potential for insight into the temporal distribution of vegetation types (Svenning et al., 2011). Additionally, the use of species distribution modelling (SDMs) (Guisan and Zimmermann, 2000; Guisan and Thuiller, 2005), provides an effective way of reconstructing the likely spatial extent of species’ ranges during glacial periods. And the combined use of SDMs and population genetic methods, broadens our understanding of past and current population dynamics, providing a stronger base for both modelling and interpreting future projected population responses to climate change (Łabiszak et al., 2021).

The primary aim of this study was to reconstruct the paleo-historical macroclimate of southern Africa to determine regions in which A. dichotomum may have persisted during natural cycles of climate change to broaden our predictive understanding of how it may respond to future projected anthropogenic climate change. The objective of this approach is to establish the potential value of this species as an indicator (sentinel) for anthropogenic climate change, and as a subject for considering, developing and testing pragmatic conservation solutions to conserve genetic diversity under anthropogenic climate change. Our hypothesis is that the reconstructed macroclimate would reflect a poleward shift of suitable bioclimatic space for A. dichotomum post-LGM, from a more equatorward distribution during the late Pleistocene. Further, we predict that the current geographic patterns of genetic variation will reflect higher genetic diversity in the regions that remained climatically suitable for the persistence of A. dichotomum during the LGM, with relatively lower variation outside of this space. We predict that in response to future warming, suitable habitat of A. dichotomum will shift further polewards as surface air temperatures increase in southern Africa. We finally predict that the migration rate required to keep pace with optimal climate space under anthropogenic climate change would exceed the rate needed to track post-LGM climate change (e.g., see IPCC, 2014).

Materials and Methods

Species and Climate Data

Species Occurrence Data

To construct a distribution data set for A. dichotomum, three separate data sources were used for presences, including data from Foden et al. (2007), occurrence points documented by W. Foden in March 2018, and presence localities from the Global Biodiversity and Information Facility Database (GBIF¹). All presence data were collated into a single presence dataset with any duplicate or outlier coordinates removed. The final dataset consisted of 367 presences (Figure 1).

Climate Data

Nineteen bioclimatic variables were obtained from WorldClim² (version 1.4; Hijmans et al., 2005) at 10 min resolution. The variables consisted of summary statistics for temperature and rainfall at different temporal resolutions, and which represent average climatic conditions from the period of 1950–2000, interpolated from weather station data. These data were used as “current” climate data to develop the SDM for the species (Hijmans et al., 2005).

The palaeoclimatic datasets (Last Glacial Maximum and mid-Holocene) and future climate scenario (2070) were also obtained from WorldClim, downscaled data from Global Climate Models (GCMs). Ten-minute spatial resolution datasets were used for the Last Glacial Maximum (22,000 years ago), the mid-Holocene (6,000 years ago) and 2070. For both the LGM and mid-Holocene, the Community Climate System Model (CCSM4) (Gent et al., 2011), the Max-Planck-Institute Earth System Model (MPI-ESM) (Giorgetta et al., 2013) and the Model for Interdisciplinary Research on Climate (MIROC-ESM) (Watanabe et al., 2011) are used. These GCM models are used as they are the three GCM’s available on Worldclim for all three time periods (LGM, mid-Holocene, and 2070). All three models were used in the SDM simulations to allow for comparison between models and to assess the agreement in the species ranges predicted. The climate dataset used to represent the future scenario in 2070 are the same GCMs as used for the past reconstructions, taken from the IPCC 5th Assessment Report climate projections (IPCC, 2014), using the representative concentration pathway (RCP 4.5) scenario. The RCP 4.5 is an intermediate “stabilisation pathway” in which radiative forcing is curbed at 4.5 W/m² after 2100 where CO₂ levels stabilise at 650 ppm (Clarke et al., 2007; Moss et al., 2008). The pathway was used as it represents a mid-range global mitigation response (Thomson et al., 2011).

WorldClim data show western reaches of southern Africa as arid during the LGM while robust proxies (Stuut and Lamy, 2004) strongly suggest an increase in rainfall during this time, caused by a northern displacement in the westerly winds. As a consequence, we focused primarily on temperature predictors due to uncertainty relating to rainfall reconstruction for the LGM, and incorporated precipitation only via composite bioclimatic variables. Temperature variables are known to be important from previous work that has shown that temperature plays an important role in limiting A. dichotomum distributions (Foden et al., 2007; Guo et al., 2016). The final variables were selected using principal components analysis (PCA), visualised using correlation circle plots and informed by Pearson’s correlations between pairs of variables. The variables selected were those which both minimised the correlation coefficient and have been identified as biologically important for the species (Guo et al., 2016).

Variables selected included temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, temperature annual range, mean temperature of the wettest quarter, mean temperature of the driest quarter, mean temperature of the warmest quarter, mean temperature of the coldest quarter and precipitation of the coldest quarter.

Modelling Framework for Current Climate

Simulations of the distribution of A. dichotomum were conducted using an Ensemble Platform for Species Distribution Modelling, “Biomod2” in R version 3.5 (Thuiller et al., 2016, 2009; R Core Development Team, 2018). Four algorithms were used including three regression models (generalised linear model, GLM; generalised additive models, GAM; and generalised boosted model, GBM) as well as a tree-based method (Random Forest, RF). The GLM was run with quadratic terms and first order interactions, the GAM was run using the mgcv package (Wood, 2017), and the GBM was limited to 1000 trees. The models were calibrated using 80% of the presence data and then evaluated on the remaining 20% (Araújo and New, 2007). In addition, 2000 pseudo-absences were used, and a threefold cross validation was done to yield an average model for each algorithm (4 model algorithms, 3 cross-validations, and 4 repetitions). Therefore, in total, 48 model simulations were run.

The models were then evaluated according to two metrics, namely the area under the Receiver Operating Characteristic (ROC) Curve and the True Skill Score Statistic (TSS) (Allouche et al., 2006). Models are considered to be credible when ROC scores > 0.8 (Swets, 1988) and TSS scores > 0.65. TSS scores varying between 0.4 and 0.8 are considered as fair to good performances (Li et al., 2016; Shrestha et al., 2018). Sensitivity and specificity were also calculated to assess predictive accuracy of the models. Sensitivity is the proportion of observed presences that are predicted as presences and specificity is the proportion of observed absences that are predicted as absences (Allouche et al., 2006). The models were calibrated onto “current” conditions using the same environmental predictors that were used to build the models. Binary (presence/absence) projections were made using the threshold that maximises the ROC score, providing a process which allows the projections to be mapped onto the African climate.

Models often yield different predictions and may vary depending on the differing scenarios. A solution for this variation is the use of multiple models within an ensemble approach, which may then allow a consensus scenario to be achieved (Araújo and New, 2007; Marini et al., 2009). Ensemble models were conducted using two methods, namely committee averaging (CA) and weighted means (WM). These models incorporate all algorithms, all pseudo-sampling and all cross-validations, subsequently producing a coefficient of variation in order to demonstrate whether or not the predictions are consistent. Ensemble models allow for better informed projections to be made when compared to a single model (Araújo and New, 2007). The ensemble models were then projected onto the current African climate which produces both continuous projections (habitat suitability) and binary projections. The ensemble models were also evaluated using TSS and ROC scores and the values were found to be higher in the ensemble modelling in comparison to the individual models (Tables 1, 2).

TABLE 1

Table 1. Statistical validation of the four models using two metrics: Area under the ROC curve (ROC) and True Skill Statistic (TSS).

TABLE 2

Table 2. Evaluation scores of the two ensemble models (CA: committee averaging and WM: Weighted Mean) based on ROC scores.

Palaeoclimate Reconstruction and Future Climate Projections

For the past and future distributional projections, the WorldClim paleo- and future climate data were sourced (version 1.4; Hijmans et al., 2005; see text footnote 2). In order to delineate the range changes of the A. dichotomum at each time period an ensemble approach was used. The ensemble models that were calibrated under the current conditions were then projected to the LGM, mid-Holocene and 2070 using the nine bioclimatic predictors from the three GCMs. These models incorporated the GLM, GAM, and RF, in addition to all pseudo-absence sampling and all cross-validations. The models were evaluated using the ROC score and were only kept if the score > 0.8. Two types of ensemble models were used, as above, namely CA and WM. From the models, binary projections were created. This allows for a visualisation of the presence and absence of A. dichotomum across the range (Figure 2). From this we can determine where range was lost and gained at each time projection.

FIGURE 2

Figure 2. Projected ranges of Aloidendron dichotomum for the (A,D) Last Glacial Maximum (22,000 years BP), (B,E) mid-Holocene (6,000 years BP), and (C,F) 2070. Ensemble models were computed using the (A–C) weighted means and (D–F) committee averaging methods using the CCSM4 global circulation model. Orange regions denote areas where the species is absent; dark blue where it has gained range; light blue where it has remained present when compared to current range projections; and red where range is unsuitable when compared to current range.

All SDM predictions were visualised in Quantum-GIS 3.2.1³ (QGIS). To calculate migration rate, coordinates were obtained from the centroid position of each range for each GCM at each time frame. Range shifts were calculated as the distance between current centroid and centroid of projected range (LGM, mid-Holocene, or 2070) and subsequently, migration rates were estimated. Migration rates were estimated using the distance between centroids as a function of the difference in time between current (present) and the other modelled time periods (LGM/mid-Holocene/2070).

Sample Collection and Amplified Fragment Length Polymorphisms Analysis

To generate a preliminary assessment of genetic variation in A. dichotomum and test concordance with models of the species’ past distribution, we sampled multiple localities across the species range allowing for inference of broad-scale patterns (detailed in Foden et al., 2007). Thinly sliced sections from the tips of young leaf samples were taken from 3 to 5 adult plants per sampling locality and collected into silica gel for drying prior to DNA extraction. Total genomic DNA was extracted using a standard CTAB extraction method for plant material (Doyle and Doyle, 1987) after grinding the dried sample to a fine powder with liquid nitrogen. A final total of 57 individuals from 17 localities were included in the analysis (Table 4).

TABLE 3

Table 3. Aloidendron dichotomum latitudinal range shift values for the Last Glacial Maximum (LGM), mid-Holocene, and 2070, in relation to the current range.

TABLE 4

Table 4. Sampling locations and estimates of genetic diversity across the range of Aloidendron dichotomum.

Amplified fragment length polymorphisms (AFLP) reactions were performed following the general method of Vos et al. (1995) using the Applied Biosystem AFLP Plant Mapping Kit Protocol. Approximately 100 ng of DNA from each sample was used in the AFLP reactions. After digestion with EcoR1 and Mse1, adaptors were ligated on to both ends of the DNA fragments and a two-step selective amplification was carried out. We tested four primer combinations for the final selective amplification. Reactions were repeated 3–5 times for five test samples to assess the number, quality, and reproducibility of AFLP fragment peaks. Final analysis of samples was based on amplification with EcoR1-ACA/Mse1-CAT (FAM). The amplified products were run with GeneScan 500 ROX size standard on an ABI 377 Genetic Analyzer (Applied Biosystems) and analysed using GeneScan v3.2 software (Applied Biosystems).

The genetic diversity estimators% polymorphic fragments (P) and Nei’s gene diversity (h) (Nei, 1973) were calculated in Genalex v6 (Peakall and Smouse, 2006) and AFLP-SURV v1.0 (Vekemans et al., 2002). The 5% criterion was applied before analysis; loci with fragment frequencies of >95% and <5% were removed. Nei’s h was then mapped to the sampling area, interpolating values onto the map space between geographical sampling locations using inverse weighting of the distance to the three nearest neighbour geographical centres in ArcView GIS (ESRI). Known populations (from Foden et al., 2007) are indicated on Figure 3.

FIGURE 3

Figure 3. The distribution of genetic variation in Aloidendron dichotomum based on Nei’s gene diversity h (NGD) calculated from amplified fragment length polymorphisms (AFLPs). Triangles represent genetic sampling locations; circles represent known populations surveyed in Foden et al. (2007).

Results

Climatic Variable Selection

The current distribution of A. dichotomum could be modelled using nine bioclimatic variables selected through PCA and correlation coefficients. The relative importance of climate variables was evenly spread for the GLM, GAM, and the RF algorithms (Figure 4). “Mean temperature of the coldest quarter” and “Temperature seasonality” were important in all three of these models. The GBM, however, was heavily driven by one variable, namely “mean temperature of the coldest quarter,” with only “minimum temp of the coldest month,” and “temperature seasonality” having some importance (Figure 4).

FIGURE 4

Figure 4. The relative importance of the nine bioclimatic variables for each model algorithm used to model the current distribution of Aloidendron dichotomum. The four algorithms used include the generalised linear model (GLM), generalised additive model (GAM), generalised booted model (GBM), and random forest (RF).

Current Climate Projections Using Individual Models and Model Performance

Results of the independent models had high congruence. All ROC and TSS values were found to be above 0.8 and 0.65, respectively (Table 1 and Supplementary Figure 1). TSS scores indicate that GAM and RF are the models that have the highest accuracy and ROC scores indicate that GBM and RF are the most accurate. The TSS score is often seen as a more practical and realistic method of evaluating the accuracy of SDM methods (Shabani et al., 2018). However, ROC scores were used as their accuracy result was found to be higher than the TSS score.

The realised current distribution of A. dichotomum was generally well predicted by the models, except for the GBM. Mapping the current model projections using ROC scores (to align with most previous work published for southern Africa) showed that the GBM model significantly overestimates the potential current distribution of A. dichotomum (Figure 5), likely the due to the severe overestimation of the range, as a result of being dominated by only two bioclimatic variables, while the RF model is more constrained. The GBM was thus removed from further analysis.

FIGURE 5

Figure 5. Current modelled geographic range of Aloidendron dichotomum. Ensemble models of the generalised linear, generalised additive, and random forest models were computed using the (A) committee averaging and (B) weighted means methods. Colour scale indicates density estimates.

Current Climate Projections Using Ensemble Models

Models with ROC scores of > 0.8 were used to build ensemble models using committee averaging (CA) and weighted mean (WM) from h. Model evaluation scores for both CA and WM were >0.98 and thus performed better than the individual models in projecting the current distribution of A. dichotomum (Table 2), noting some apparent range overestimation in the southern and eastern ranges of the species relative to recorded presences at these areas (Figure 2). Distant disjunct bioclimatic conditions identified by the ensemble model with no known occurrence data for the species were not considered further.

Palaeoclimatic Distribution

Using the ensemble model predictions, suitable bioclimatic space for A. dichotomum was modelled for LGM climate, and found to be strikingly different from the current modelled distribution. More than 90% of the species’ current range was excluded as bioclimatically suitable, suggesting a significant poleward latitudinal shift for this species since the LGM informed by modelled bioclimatic suitability (extending to 5.8 degrees south of LGM modelled range). During the LGM, models predict that 69% of suitable habitat was outside the current range, mostly inside Namibia (Figure 2). Using an average latitudinal distance from the centroid of the current projected range, suitable bioclimatic space for A. dichotomum was around 650 km to the north of the current modelled range and around 320 km to the west (Table 3). This implies that, provided the species has tracked changes in climate suitability, A. dichotomum has migrated 650 km in the poleward direction in the last 18,000 years. As such, its poleward migration rate since the LGM has been ca. 0.4 km/decade.

The mid-Holocene model predicts that suitable bioclimatic space for A. dichotomum was also contracted in comparison to the current range (≈62%). When comparing modelled shifts in climate suitability between the LGM and mid-Holocene, it is clear that suitable habitat for the species shifted significantly poleward (5.9°S) subsequent to the end of the LGM (Table 3). In comparison to the current range, suitable habitat for A. dichotomum during the mid-Holocene, despite contracting, is suggested to be found in similar locations and predicted to be on average only 40 km north and 52 km to the west of the current range.

Future Distribution

By 2070 suitable bioclimatic space for A. dichotomum is projected to shift 191 km (average of the three GCMs) eastwards (Table 3). The northern reaches (above 23.7 degrees) of the current range are expected to become unsuitable in the future and this range will likely be lost. Populations to the west of 14.7 degrees are also projected to be lost (Figure 2). The models differed slightly in the latitudinal shift with the CCSM4 and MIROC-ESM models projecting a poleward range shift, whereas MPI-ESM projection suggests a slight equatorward shift (Table 3). Nevertheless, on average, the distance that the species will be required to migrate is 42 km in a 70-year period. This implies that the migration rate required for A. dichotomum to reach its predicted optimal range by 2070 is 6 km per decade.

Spatial Patterns of Genetic Variation

The final data set comprised 251 AFLP loci analysed as a binary matrix. Mean% polymorphic loci (P) across the data set was 39.8% (std. dev. ± 13.3; range 17.1–65.7%) and mean Nei’s h was 0.36 (std. dev. ± 0.03; range 0.29–0.41) (Table 4). The spatial representation of Nei’s h indicates generally higher levels of genetic variation in the central and more northerly (equatorward) part of the species distribution, with lower levels characterising the more southern and south-easterly regions (Figure 3). At the most northerly extent of the species range, individuals at higher altitude on the Brandberg Massif in Namibia support greater levels of genetic variation as measured by Nei’s h, than their nearest neighbour counterparts on the lower lying plains.

Discussion

Unequivocal attribution of observed range shifts to anthropogenic climate change is rare, as pointed out by Taheri et al. (2021). An argument has been made that agreement between model predictions and observation is now sufficient to support high confidence in the conclusion that anthropogenic climate change is driving species range shifts (Parmesan et al., 2013). However, the incidence of false attribution for individual studies due to range shifts that mimic anthropogenic climate as a driver (Hockey and Midgley, 2009) is not known. In long-lived species, for example, the legacy of historical climatic, or other non-anthropogenic, events can be used as evidence for their influence and challenge the attribution of observed population trends to anthropogenic climate change. For A. dichotomum, matching of demographic rates for multiple populations with observed climate shifts (e.g., Foden et al., 2007) in addition to observed geographical patterns of mortality serves to reduce uncertainty in attribution, but nonetheless alternative explanations for these patterns have been proposed (Jack et al., 2014). Understanding and predicting species’ responses to ongoing and future climate change can be enhanced through an understanding of species’ biological responses to deeper historical climatic trends (Savolainen et al., 2007; Lawing and Polly, 2011). It is therefore extremely useful to use the example of A. dichotomum to develop a paleo-historical view on recent observed demographic changes, and ongoing and potential range shifts.

During periods of glaciation, exemplified most recently by the Last Glacial Maximum (LGM; approx. 18,000–22,000 years BP), it has now been well established that species ranges across the world shifted and/or contracted equatorwards, persisting in areas now often referred to as climatic refugia (e.g., Stewart et al., 2010; Lyam et al., 2020). Range expansions occurred in many species when temperatures warmed rapidly (e.g., during the end Pleistocene, 18,000–9,000 years BP) and more slowly towards the mid-Holocene thermal optimum (approx. 6,000 years BP) (Hewitt, 1999; Lyam et al., 2020). While studies investigating the role of the palaeoclimatic oscillations in shaping modern geographic and genetic patterns are generally biased towards the Northern hemisphere (Rödder et al., 2013; Roces-Díaz et al., 2018; Łabiszak et al., 2021), patterns that they identify indicate the general evolutionary fingerprint that may be left by paleo-range shifts. For example, Łabiszak et al. (2021) attribute current phylogenetic patterns of P. uliginosa to an ancient bottleneck event that took place around 26,400 years ago as a result of a range contraction during the LGM. Roces-Díaz et al. (2018) showed that Castanea sativa contracted its range to Western Europe during the LGM, followed by a subsequent expansion during the mid-Holocene, a shift that was corroborated by phylogeographic analysis.

In southwestern Africa, surface air temperatures likely rose by about 5°C over 10,000 years during glacial to interglacial warming at the end of the Pleistocene (Potts et al., 2013), and there is evidence for regional drying associated with the poleward retreat of rain-bearing westerlies (Stuut and Lamy, 2004). These significant changes in bioclimatic conditions likely had significant impacts on the regional semi-arid winter-rainfall flora, as revealed by phylogenetic evidence for the widespread species Elytropappus rhinocerotis (Bergh et al., 2007), palynological evidence (Scott, 1994; Scott et al., 2018) and are also thought to have driven a poleward expansion of woody vegetation in the summer rainfall sub-tropics (Scott, 1999). However, most palynological studies in the region have been point based and have not reconstructed individual species’ range shifts. A study that investigated the impact of palaeo-climatic changes on a southern African savanna tree species, S. Senegal, also suggests that the patterns of range contractions during the LGM, followed by range expansions through the mid-Holocene until current observed distributions seen for Northern Hemisphere species are likely to have also occurred for Southern Hemisphere species (Lyam et al., 2020).

There is strong bioclimatic modelling evidence that the current geographic range of A. dichotomum is limited by air temperature variables (Guo et al., 2016), and our modelling focussing on temperature as the primary driving variable showed the likely importance of low temperature limits to range edges in the southerly (poleward) and easterly reaches of A. dichotomum’s range (Figure 4). According to our SDM, areas of suitable climatic conditions for the persistence of A. dichotomum are projected to have contracted equatorward during the LGM when air temperatures would have decreased significantly in this region (Figure 2), thereafter, suitable climatic space is projected to have shifted southwards at the end of the LGM and into the mid-Holocene when temperatures became more favourable (Figure 2). Indirect evidence of this pattern can be seen in preliminary observations of geographic disparities in the genetic diversity of the northerly and southerly ranges of A. dichotomum (Figure 3). A general cline in variation characterises the species contemporary range, with northerly populations supporting, on average, higher levels of genetic diversity than populations in the southern and south-eastern regions.

When viewed together, these independent lines of evidence support the hypothesis of a glacial contraction to a northerly refugial state during the LGM (Stewart et al., 2010). Similar patterns have been found for Chinese pine (Pinus tabulaeformis), where ecological niche modelling and the presence of glacial lineages provide evidence for post-glacial spread from LGM refugia (Hao et al., 2018). Midgley et al. (2005) suggested that Pleistocene glacial cycles have played an important role in the creation of species clines, and in the current geographic patterns of high endemic richness of clades with low dispersal ability in the winter-rainfall western areas of southern Africa. Huntley et al. (2016) invoked a similar explanation for current geographic patterns of bird species richness in this region, and more broadly in southern Africa. Certainly, the impact of changing climate patterns through the LGM and mid-Holocene until the present day has been important for current geographic distributions and spatial patterns of genetic diversity for multiple long-lived tree species globally (Hao et al., 2018; Lyam et al., 2020; Yannic et al., 2020; Łabiszak et al., 2021).

What insights can be gained for conservation of A. dichotomum in the face of projected climate change? Modern plant taxa have persisted since at least the Pleistocene by both in situ adaptation and geographic range shifts in response to climate change related to glacial-interglacial cycles (Dawson et al., 2011), with evidence in the winter rainfall region of southern Africa for range shifts in E. rhinocerotis (Bergh et al., 2007). The past changes in climatic suitability projected by this study suggests that A. dichotomum may also have shifted its range since the LGM. There are four potential responses for plant species under climate change, namely persistence through migration to new ecological niches that are more favourable, persistence through in situ adaptation by genetic changes, persistence through in situ acclimation by phenotypic plasticity, or local extinction (Aitken et al., 2008; Feeley et al., 2012).

While species with shorter generation times are likely to be able to adapt via evolutionary processes (Lustenhouwer et al., 2018), it has been suggested that long-lived trees with extended juvenile phases are unlikely to be able to match the rate of anthropogenic climate change through in situ adaptation alone (Nathan et al., 2011; Kremer et al., 2012). Feeley et al. (2012) even suggest that, due to the high rates of projected anthropogenic climate change, adaptation to climate change by long-lived tree species may even be ruled out altogether. It is believed that migration is the most important autonomous response to climatic change for tree species (Feeley et al., 2012), allowing them to remain within their inferred climatic niche (Pitelka, 1997). However, the response of long-lived species to changing climates is projected to be characterised by significant lags, with species only reaching equilibrium with suitable climate hundreds of years after climate stabilisation (Kuparinen et al., 2010; Moran and Ormond, 2015). It is therefore likely that species will have to simultaneously shift range and adapt and acclimate to new local climatic conditions (Feeley et al., 2012; Moran and Ormond, 2015). Further, it has been suggested that adaptation may become important for long-lived tree species as climate-induced mortality of mature trees increases (Kuparinen et al., 2010). With the observed increase in mortality of mature A. dichotomum individuals in warmer populations (Foden et al., 2007; van der Merwe and Geldenhuys, 2017), and with its inherent long generation time, it is possible that evolutionary adaptation may be important to consider in future projections of changes for this species. Shifting range may therefore be only one part of the whole response for A. dichotomum to future projected changes in climate, with local adaptation and acclimation playing an important role as well. The ability of A. dichotomum to adapt to these changes in situ, will require genetic diversity within populations, where some individuals carry genes allowing them to better cope with extreme climatic conditions (Foden et al., 2019).

From the modelling conducted in this study, it can be inferred that a significant pre-Holocene range contraction likely occurred for A. dichotomum. For the LGM period, the ensemble models showed spatial agreement of the distribution of A. dichotomum during this time (Figures 2A,D and Table 3). Mid-Holocene projections showed congruent patterns across ensemble models, with the area of suitable climatic space for A. dichotomum projected to have been almost 6° south of that during the LGM (Figures 2B,E and Table 3). The models therefore suggest that, if A. dichotomum did indeed track these changes in climate suitability, the past range of A. dichotomum would have been found significantly equatorward (650 km north) of its current range and that the species was non-existent in South Africa at this time, with the only suitable habitat located in the higher latitudes of Namibia (North of 27°S).

Range expansions in this species are likely to be caused by rare long distance dispersal events of its wind dispersed seed, followed by exponential population growth (Hampe and Petit, 2005) and potential founder effects. Signatures of this process are predicted to include generally lower levels of genetic diversity in the poleward direction (Figure 3). The latitudinal trend in genetic diversity of A. dichotomum (Figure 3) supports our current conclusion of a paleo-historical poleward range shift from more spatially stable equatorward populations, into novel geographic locations when temperatures warmed significantly at the end of the Pleistocene (Potts et al., 2013). This pattern of genetic diversity may also allude to the importance of microclimatic refugia during the LGM, for example the Brandberg Massif in Namibia which currently supports higher genetic variation than the surrounding lowlands (Figure 3). Similar patterns of higher genetic variation at higher elevations have been found for Chinese pine (P. tabulaeformis) across the Loess Plateau in northern China as a result of persistence of the species in mountains surrounding the plateau during the LGM, and subsequent recolonisation of surrounding areas from local refugia (Hao et al., 2018).

A concern for slow-growing and long-lived tree species is whether they will be able to shift range fast enough to keep pace with the ever-increasing rates of anthropogenic climate change (Feeley et al., 2012), and the likelihood of significant lags of range shifts behind changing climates (Moran and Ormond, 2015). The projected post-glacial increase in climatic suitability in the poleward direction and the present pattern of populations isolated from each other by many tens of kilometres suggest that long-distance dispersal and establishment is possible for A. dichotomum at least on a century time scale. From the modelled climatically determined southern range limit at the LGM, it is possible to calculate a hypothetical migration rate of 0.4 km per decade to permit the current range extent, assuming range expansion in response to warming post-LGM. This rate is far lower than the average rate of ∼2 km/year suggested for other long-lived tree species (Settele et al., 2014), thus plausibly supporting a post-LGM range expansion of this magnitude for this species. However, this hypothetical rate is 15 times slower than the migration rate required for the species to track its optimal climatic space by 2070, as modelled in this study.

Models developed here using mainly temperature controls only were able to project the current range extent of A. dichotomum with acceptable accuracy (Figure 5), but it must be noted that the poleward range extent was slightly overestimated when compared to current observations of A. dichotomum populations (see Figure 1). It may therefore be that the more southerly reaches highlighted by the models may be suitable for A. dichotomum establishment, but that there is an inherent dispersal lag behind increased climate suitability (Moran and Ormond, 2015), and the species may currently be undergoing a range-filling process (Schurr et al., 2007) in these southern sites where high recruitment rates have been recorded. Understanding ongoing shifts in range at this southerly range edge will be vital in interpreting and monitoring range expansion and correctly attributing observed changes for further learning about autonomous climate change responses in this species. It will therefore be important to conduct empirical studies to pinpoint the role of low temperature relative to other controls on at least the southerly range edge of the species.

Current demographic trends across the species’ range have also been used to infer an incipient range shift (Foden et al., 2007). The southern and eastern populations show high population growth rates (Foden et al., 2007; Hoffman et al., 2010; van der Merwe and Geldenhuys, 2017) while northern and central populations (equatorward) tend to have higher mortality rates and low recruitment rates (Foden et al., 2007). These patterns have been attributed to warming and drying trends (Foden et al., 2007), but the attribution was subsequently challenged by an analysis of demographic patterns that split populations by arbitrarily defined climatic zones (Jack et al., 2016), thus making a direct test of the findings of the original attribution study unclear.

The modelled LGM to mid-Holocene (Figure 2) and then to current (Figure 5) range shifts of A. dichotomum align with initial observations of geographic disparities in the distribution of genetic variation in the more northerly and southerly populations, with a general decline in genetic diversity in the poleward direction (Figure 3). The palaeogeographic modelling, together with this preliminary genetic data, suggest that the northern populations are highly likely to be older and, as such, acted as the source region for the species’ southward range expansion (Hampe and Petit, 2005). A more detailed and densely sampled phylogeographic analysis of the species will allow a robust test of this hypothesis, providing an opportunity to reconstruct both the fine-scale and regional genetic structure of the species, and estimate historical migration rates across the distribution. The linkage between historical climatic events and the current pattern of genetic variation suggests that the phylogeographic history of A. dichotomum has likely played an integral part in shaping the current population demographics. A more nuanced understanding of this history will provide valuable insight for its importance to the species in a future of substantial environmental change.

What might be projected for the next few decades to centuries for this species? Surface air temperatures are projected to continue to increase across southern Africa over the next decades, with likely increases of 1–2°C along coastlines and 2–3°C in the interior (Department of Environmental Affairs, 2017), with extreme increases of up to 4–6°C over the western region under high emissions scenarios (Engelbrecht et al., 2011). Warmer temperatures will result in increased evaporation and therefore decreased water availability, especially if rainfall amounts are reduced (Schulze et al., 1993; MacKellar et al., 2014). Model projections to 2070 were consistent across ensemble models (Figures 2C,F and Table 3), and broadly consistent with previous distribution modelling for the species (Foden et al., 2007; Guo et al., 2016), in projecting significant loss of bioclimatic space in the equatorward parts of the range. However, this study suggests, in concurrence with Foden et al. (2007) and Guo et al. (2016), that the majority of the current poleward range will remain within the optimal climatic niche for A. dichotomum at least during this century, but in contrast to Foden et al. (2007), that additional bioclimatic space would become available towards the east (Figures 2C,F). The projected persistence of climatic suitability in some equatorward parts of its range into the future may be due to climatic buffering provided by the prominent elevational gradients in these regions (Foden et al., 2007) and could potentially conserve unique genetic variation harboured by these populations (Figure 3).

If successful conservation plans and strategies are to be developed to support natural adaptation to climate change, it will be important to firstly understand what capacity these organisms have to adapt on their own, and then subsequently to anticipate biological impacts on such organisms, in order to make the best use possible of natural processes to enhance species, community, and ecosystem resilience. While the models are able to predict the optimal niche space for A. dichotomum in the future, it is critical to take into account its life history traits, such as the fact that the species is very slow growing and is slow to reach sexual maturity (Kaleme, 2003). This will provide further insight as to whether the species will be able to migrate to this optimal niche space in time to track the changing climate, which is dependent on its dispersal capability (Foden et al., 2019).

The iconic desert tree succulent A. dichotomum provides a powerful test of the theory behind plant species responses to climate change, and its application to conservation action into the future. While poleward populations occupy a part of the potential range which we project to remain suitable in 2070, slow migration rates in potential expansion zones in the south and east (leading range edges), coupled with projected loss of climatic suitability in the north (trailing range edges), will likely narrow the geographic range of A. dichotomum (Figure 1). The generally lower levels of genetic variation of poleward populations may also be important in limiting the species’ potential to survive and adapt to a novel climate (Foden et al., 2019). With our preliminary genetic analyses suggesting higher genetic diversity in the trailing edge populations, they may become disproportionately important for the longer-term conservation of the phylogeographic history and adaptive capacity of the species, warranting a higher level of conservation priority in the near future (Hampe and Petit, 2005).

Further research will also expand our current understanding of both the large and fine scale drivers of genetic structuring of A. dichotomum across its geographic range. This could provide invaluable insights into the effects of historical and ongoing climate change on current genetic structuring of the species (Avise, 2000; Byrne, 2007) and strengthen a predictive framework for understanding the future response of A. dichotomum to anthropogenic climate change (Gavin et al., 2014). In particular, SDM-based efforts for reconstructing the paleo-range for this species could be improved through the inclusion of credible paleo-historical rainfall patterns that concur with empirically based reconstructions (Stuut and Lamy, 2004), at the local scales necessary to construct SDM projections. Finally, empirical observations and experimental work on this species holds promise for testing SDM assumptions and attributing local scale changes in recruitment and mortality of this species in response to ongoing observed changes in climate.

Conclusion

The link between past, present, and future may be extremely valuable when investigating climate change risks to species under anthropogenic climate change. While developing such information for all species, or even several species, is an impossible task, information for selected species in multiple areas may add significantly to an understanding of risks and potential adaptive responses. This longer-term perspective is especially important for understanding the range dynamics of a long-lived sessile organism. The preliminary population genetic data, combined with range modelling in this study suggests that the range of A. dichotomum was likely limited to northern (equatorward) reaches relative to its current distribution during the LGM and subsequently expanded polewards into the Holocene. Model projections to 2070 showed a predominantly eastward as opposed to the hypothesised poleward range shift, and diverged from those previously projected (e.g., Foden et al., 2007). Overall, the outlook for the persistence of this species in the wild is positive, as most of its current range appears to remain habitable into the future. However, the northerly (equatorward) populations appear to be at risk of local extirpation, as the calculated migration rate of the species is unlikely to permit them to track the projected warming rate over the next few decades. If this occurs, valuable adaptive genetic diversity of populations possibly selected under warmer and drier conditions may be lost. Beyond 2070, if emissions continue to rise and climate to warm, populations in easterly and poleward reaches of the range may be at risk that would be exacerbated due to their lower genetic diversity.

We would propose that proactive ex situ efforts could be made to ensure safeguarding of this genetic diversity, such as through carefully considered assisted colonisation informed by work that elaborates on this study. The findings presented here add to the growing literature on this iconic species that will inform conservation responses not only for this species, but potentially also for a wide range of endemics and near endemics of this species-rich region.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

LB and GM designed the study. LB led species distribution modelling. JB led the genetic analysis. K-AG led manuscript preparation. All authors contributed to writing.

Funding

GM acknowledges funding from the National Research Foundation Grant 118591 and Global Environment Facility grant GEF-5810 “Spatial Planning for Area Conservation in Response to Climate Change (SPARC).”

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

GM and K-AG gratefully acknowledge Hendrik van Zijl of Nieuwoudtville for his logistical support during field work expeditions. LB thanks Professor Wendy Foden for the discussion, comments, and support during the initial development of the script. LB would also like to thank Nicola Stevens for guidance and assistance with analyses conducted using geographic information systems (GIS) software. JB thanks Greg Hughes for assistance with the spatial representation of genetic variation. The authors thank Wendy B. Foden for collecting the leaf material used in genetic analyses. The authors thank the editor and two reviewers for their valuable contributions in shaping the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.715702/full#supplementary-material

Supplementary Figure 1 | Projections of the current distribution of Aloidendron dichotomum in southern Africa using the (A) generalised linear model, (B) generalised additive model, (C) generalised boosted model, and (D) random forest algorithms. Colour scale indicates density estimates (see Supplementary Data Sheet 1).

Supplementary Data Sheet 1 | Locality data (X:Y latitude:longitude decimal degrees) for Aloidendron dichotomum used for distribution modeling in this paper.

Footnotes

References

Aitken, S. N., Yeaman, S., Holliday, J. A., Wang, T., and Curtis-McLane, S. (2008). Adaptation, migration or extirpation: Climate change outcomes for tree populations. Evol. Appl. 1, 95–111. doi: 10.1111/j.1752-4571.2007.00013.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, C. D., Breshears, D. D., and Mcdowell, N. G. (2015). On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 6:129.

Google Scholar

Allouche, O., Tsoar, A., and Kadmon, R. (2006). Assessing the accuracy of species distribution models: Prevalence, kappa and the true skill statistic (TSS). J. Appl. Ecol. 43, 1223–1232. doi: 10.1111/j.1365-2664.2006.01214.x