Projecting Podocarpaceae response to climate change: we are not out of the woods yet

Abstract Under the changing climate, the persistence of Afrotemperate taxa may be threatened as suitable habitat availability decreases. The unique disjunct ranges of podocarps in southern Africa raise questions about the persistence of these species under climate change. Here, we identified likely environmental drivers of these distributions, characterized the current and future (2070) environmental niches, and projected distributions of four podocarp species in South Africa. Species distribution models were conducted using species locality data for Afrocarpus falcatus, Podocarpus latifolius, Pseudotropheus elongatus and Podocarpus henkelii and both historical climate data (1970–2000) and future climate scenarios (Representative Concentration Pathway [RCP] 4.5 and 8.5, 2061–2080) to estimate the current and future distributions. We also used this opportunity to identify the most important climatic variables that likely govern each species’ distribution. Using niche overlap estimates, a similarity test, and indices of niche expansion, stability and unfilling, we explored how niches change under different climate scenarios. The distribution of the study species was governed by the maximum temperature of the warmest month, temperature annual range, mean temperature of the wettest quarter, and precipitation of the wettest, driest and warmest quarters. The current distribution of A. falcatus was predicted to expand to higher elevations under RCP 4.5 and RCP 8.5. Podocarpus henkelii was predicted to lose most of its suitable habitat under RCP 4.5 and expand under RCP 8.5; however, this was the opposite for P. elongatus and P. latifolius. Interestingly, P. elongatus, which had the smallest geographic distribution, showed the most vulnerability to climate change in comparison to the other podocarps. Mapping the distribution of podocarps and understanding the differences in their current and future climate niches provide insight into potential climate drivers of podocarp persistence and the potential for adaptation of these species. Overall, these results suggest that P. elongatus and P. henkelii may expand to novel environmental niches.


Introduction
Conifers are a crucial component of global forests with both economic and ecological importance (Farjon 2018). Podocarpaceae (podocarp) is the second most species-rich conifer family and the largest clade in southern conifers with 19 genera and 187 species (Christenhusz and Byng 2016), and are one of the few gymnosperms that inhabit tropical forests in the Southern Hemisphere. Species in the Podocarpaceae were historically centred in Gondwana, subsequently expanding to Australasia, southernmost Africa, and currently also occur in Malaysia (Leslie et al. 2012;Quiroga et al. 2016). Podocarps were prominent in Gondwana during the Cretaceous period (Krassilov 1974) diversifying by the early Triassic period into the podocarps we know today (White 1981;Mill 1999). This diversification was possibly due to the onset of warmer and wetter climates because of the opening of the South Ocean (Atlantic and Indian Ocean; Morley 2011). During the mid-Cretaceous, podocarps went extinct in West Africa during the flattening of the uplands, reduced in geographic extent in India during the late Cretaceous when India drifted northwards (Morley 2000), and subsequently dispersed into its current Southeast Asian distribution following global cooling. In Africa, global cooling may have enabled the dispersal of podocarps from East Africa to the highlands of West Africa through northern and southern pathways (Migliore et al. 2022). Presently, podocarps have a largely pantropical distribution. Some taxa extend into subtropical and temperate latitudes, where they mainly occur in Australasia, Central and South America, and tropical montane Africa (Farjon 2010;Klaus and Matzke 2020). In Afrotemperate forests, podocarps persist within fynbos, heathland and grassland matrices, where fire and grass competition are the predominant factors that influence forest distribution Lawes 2009a, 2011;Adie et al. 2017).
Remarkably, most podocarps are restricted to humid, mountainous environments, including angiosperm-dominated forests (Kitayama et al. 2011). Many extant podocarps are slow growing, fire and drought intolerant, and these characteristics put them at a competitive disadvantage when co-occurring with angiosperms (Bond 1989;Brodribb and Hill 1998), which restricts them to these 'steppingstone' refugia.
Developing management strategies and practical measures to conserve podocarps is difficult without identifying key environmental variables to each species' climatic niche and predicted geographical distribution, particularly under future climate conditions. Moreover, climate change affects habitat requirements for many species (Bakkenes et al. 2002;Aitken et al. 2008;Essl et al. 2011;Moran 2020;Fricke et al. 2022). Consequently, determining whether climate change will affect the geographical distribution and environmental niches of podocarps presents another critical challenge linked to their ecological value and significance. The current distributions of A. falcatus, P. elongatus, P. henkelii and P. latifolius are relatively well known; however, little is known about the environmental drivers of podocarp distribution, their dispersal patterns and how climate change will affect their distribution and environmental niches.
Niche conservatism is used to test how species shift (expand and/or contract) their ranges in response to climate change (Wiens 2004;Wiens and Graham 2005;Wiens et al. 2010). This ecological and evolutionary process is an issue of concern particularly due to the expected adverse effects of climate change on species distributions and diversity (Araújo et al. 2013). Therefore, modelling the environmental niche and projecting it into geographic space allows us to test whether species exhibit environmental niche conservatism as ranges shift under climate change. Testing environmental niche conservatism is necessary to enable us to develop strategies to mitigate the negative effects of climate change and provide important ecological insights. For the purpose of this study, 'niche conservatism' is defined as when the predicted environmental niche tends to remain similar to the current environmental niche (Warren et al. 2008). Eeley et al. (1999) and Colyn et al. (2020) suggested that the ability of forest species to track suitable environmental conditions under climate change could be significantly limited by anthropogenic landscape change. Therefore, in this study, the ability of podocarps to persist under climate change will be determined by the ability of the species to track favourable niches and expand their geographic range.
Species distribution models (SDMs) are mathematical algorithms that estimate a species' climatic niche by characterizing a species' occurrence in relation to environmental factors (Elith et al. 2006). SDMs can be used to predict and analyse patterns of distribution and could estimate the risk these environmental drivers pose to species (Guisan and Thuiller 2005;Elith and Leathwick 2009;Guisan et al. 2013). Species distribution models quantify the relationship between plant distributions and environmental factors (Guisan and Thuiller 2005;Elith and Leathwick 2009;Wiens et al. 2010;Guisan et al. 2013); thus, demonstrating their vulnerability more accurately. Species distribution models often assume niche conservatism, which suggests that species can maintain their environmental niches (Peterson et al. 1999;Wiens 2004;Wiens et al. 2009;Peterson 2011). It is still controversial whether species' environmental niches are preserved in space and time (Peterson et al. 1999;Broennimann et al. 2007). As a result of the ongoing theoretical development and quantification of a species' environmental niche, researchers have advanced our understanding of how species fluctuate in their requirements for and tolerance of various factors (Soberón 2007). In this study, we used a combined environmental (niche overlap indices) and geographical approach (temporal transferability of SDMs) to characterize the distribution of A. falcatus, P. elongatus, P. henkelii and P. latifolius. The objectives of this study were to (i) identify the environmental variables that shape the distribution of A. falcatus, P. elongatus, P. henkelii and P. latifolius; (ii) model the current distribution of South African podocarps; (iii) project the current models onto two future climate emissions scenarios (Representative Concentration Pathway [RCP] 4.5 and RCP 8.5); and (iv) compare the future geographic distribution of podocarp environmental niches to their current environmental niche geographic distribution. Such comparisons will enable us to infer the vulnerability of these South African podocarps to the changing environment.

Study species
South Africa consists of four podocarp species from sister genera: Afrocarpus (A. falcatus) and Podocarpus (P. elongatus, P. henkelii and P. latifolius ;Farjon 2001). In Africa, podocarps are restricted to highland archipelagos in montane forests along the Afromontane Forest belt, where they persist in small forest patches within a grassland, fynbos or heathland matrix. The Afromontane Forest belt extends from Ethiopia along the eastern mountain range, all the way to the southern Cape in South Africa. Afrocarpus falcatus is present throughout Afromontane forests and is absent in West Africa. Afrocarpus falcatus and P. latifolius are widely distributed geographically across South Africa but are restricted to Afromontane Forest habitats. In South Africa, A. falcatus and P. latifolius occur along the southwest of South Africa and extend to the south and east coasts of the country up the Drakensberg mountains in the eastern parts of the country and extend northwards in temperate midland regions towards Zimbabwe. Podocarpus henkelii occurs in forest habitats and is restricted to summer rainfall areas in the mesic east of the country. Podocarpus elongatus only occurs in the extreme western and southern parts of the country (winter rainfall) in the forest and fynbos biome.
Generally, our study included three major steps, with relevant details reported below: (i) data were downloaded, compiled and variables were selected; (ii) A. falcatus, P. elongatus, P henkelii and P. latifolius occurrence and climate data were used to characterize niche and environmental spaces; and (iii) a variety of species distribution modelling algorithms were used to project the potential distribution of A. falcatus, P. elongatus, P henkelii and P. latifolius in South Africa.

Climate variables
Current climate data were downloaded from Chelsa version 2.1 (historical  and future [2061-80) climate data were downloaded under RCPs 4.5 and 8.5 emission scenarios from the Hadley Global Environment Model 2-Atmosphere Ocean (HADGEM2-AO) circulation model which provides adequate coverage for Africa (Martin et al. 2011). These data were downloaded using 'Chelsa.CMIP_5. download' function from the ClimDatDownloadR R package (Karger et al. 2017;Jentsch et al. 2021). All environmental variables had a spatial resolution of 30 arc seconds (~1 km 2 pixel area at the equator). To avoid obtaining spurious results and over-parameterization of the models due to variable multi-collinearity (Elith et al. 2010;Synes and Osborne 2011;Boria et al. 2014), variable correlations were minimized using dimension reduction techniques. A Pearson's correlation test based on all 19 climate variables was conducted for all species' presence points and used to exclude highly correlated variables (r > 0.65) from our models. To further avoid overparameterization of the model and obtaining spurious results due to variable multi-collinearity, we used the variance inflation factor on the remaining variables to remove significantly correlated variables. Seven predictor variables were used in subsequent processing: Max Temperature of the Warmest Month (BIO5), Temperature Annual Range (BIO7), Mean Temperature of Wettest Quarter (BIO8), Mean Temperature of Coldest Quarter (BIO11), Precipitation of Wettest Quarter (BIO16), Precipitation of Driest Quarter (BIO17) and Precipitation of Warmest Quarter (BIO18). These seven predictor variables were used to model the species distribution for all the species to allow for comparisons between species and are variables used by other studies predicting the species distribution of podocarps (e.g. Bernardi et al. 2020;Tesfamariam et al. 2022).

Species distribution modelling
Species distribution modelling uses occurrence data and environmental variables to simulate suitable/habitable environmental conditions for focus species. Here, we used 10 different algorithms available within the biomod2 package (Thuiller et al. 2021) in R version 4.1.1 to obtain an ensemble of predicted distributions. An ensemble approach was used because numerous studies have shown that ensemble SDMs improve the predictive performance of SDMs in comparison to using one algorithm (Grenouillet et al. 2011;Guo et al. 2015;Bernardi et al. 2020;Zurell et al. 2020). The following algorithms were used to develop the ensemble SDMs: generalized linear model, generalized additive models, generalized boosting model/boosted regression trees (GBM), surface range envelop/BIOCLIM (SRE), classification tree analysis, artificial neural network, flexible discriminant analysis, multiple adaptive regression splines, random forest and maximum entropy. The following parameters were set for the models: one run of 500 pseudo-absence points was used on one model evaluation, split into 80 % for training and 20 % for testing for each SDM, and variable importance was determined using three permutations for the full extent of South Africa. The remaining model values were set to default values. Area under the receiver-operating curve (AUC), true skills statistics (TSS) and the continuous Boyce index (CBI) were used to assess model predictive efficiency and performance (Hanley and McNeil, 1982;Allouche et al. 2006). The Boyce index uses presence data only to measure how different the model predictions are from a random distribution of observed presence across the prediction gradient. Boyce index values range between −1 and +1, where positive values indicate that the model predictions are consistent with the presence data, values close to zero suggest that the model is no different from a random model and negative values indicate no occurrence when there is. All statistical analyses were performed in R using scripts from Broennimann et al. (2012), which are now available in the ecospat R package (Di Cola et al. 2017;Broennimann et al. 2020). The potential distribution for each species was projected under the RCP 4.5 and RCP 8.5 future climate emission scenarios for the year 2070 (average for 2061-80). Representative Concentration Pathways are models of greenhouse emission scenarios which predict future climate under different greenhouse emissions scenarios. RCPs 4.5 and 8.5 were chosen because RCP 4.5 greenhouse gas emissions are projected to peak around 2040 and then decline, while in RCP 8.5, greenhouse gas emissions are projected to continue to increase throughout the 21st century, respectively, a 'best-case' and 'worst-case' scenario (Meinshausen et al. 2011). At the current rate of emissions and mitigation processes, these are the most likely scenarios. RCP 2.6 was not included because it was unrealistic, as it assumed peak emissions by 2020 and subsequent decline requiring negative emissions (Meinshausen et al. 2011;Van Vuuren et al. 2011).
To identify the change in the geographic area of A. falcatus, P. latifolius, P. henkelii and P. elongatus binary maps were generated from SDMs, and the current and future binary models were used to calculate the change in geographic area between the current and future climate emissions scenarios. Changes in geographic area were classified as 'loss', 'gain' and 'stable'. Loss referred to the area the species originally occupied, but no longer occupies after climate change. Gain referred to occupied areas that the species did not originally occupy and then occupy after climate change. Stable referred to areas that persisted after climate change, or also referred to as climatic microrefugia (Pang et al. 2021). All geographic information system analyses were done in ArcGIS v 10.8 (ESRI).

Calculating niche overlap in environmental space
Assessing climatic niche characteristics is a powerful approach to studying niche divergence and conservatism. Broennimann et al. (2012) presented a principal component analysis (PCA) method which places environmental variables into a two-dimensional space identified by their first and second principal components. When testing niche divergence and niche conservatism hypotheses, niche overlap can provide a relatively reliable measurement of overlap. We quantified the degree of shared environmental niche space between the current and the future climate emission scenarios of A. falcatus, P. elongatus, P. henkelii and P. latifolius. In this study, niche overlap between the current and future climate emissions scenarios for each podocarp species was computed using the Schoener's D statistic (Schoener 1968;Warren et al. 2008). Schoener's D ranges from 0 (no overlap) to 1 (complete overlap). Results were then interpreted using the classification scheme suggested by Rödder and Engler (2011) where D ranges indicated potential interpretations: 0.0-0.2 = no or very limited overlap, 0.2-0.4 = low overlap, 0.4-0.6 = moderate overlap, 0.6-0.8 = high overlap and 0.8-1.0 = very high overlap/identical niche. Subsequently, niche overlap was interpreted in the context of three general niche categories: niche expansion, stability/stasis and unfilling Glennon et al. 2014;Guisan et al. 2014;Di Cola et al. 2017). Niche stability (i.e. the proportion of the current range niche overlapping with the future range niche), niche expansion (i.e. the proportion of the future range niche not overlapping with the current range niche) and niche unfilling (i.e. the proportion of the current range niche not occupied in the future range niche) enable the ability to generate hypotheses about the potential drivers of podocarp niche dynamics between current and predicted climatic emission scenarios. The niche stability and expansion values always sum to 100 %. The niche unfilling value corresponds to the expansion value when shifting current and future ranges. Niche conservatism is defined as the tendency for species to retain their environmental niche in space and time and is represented by 'niche stability'.
We used an ordination technique that applies kernel density smoothers to species presences in environmental space  to determine the environmental niche occupied by each species under current and future climate emissions scenarios. The kernel density function is applied for smoothing the density of occurrences for each of the generated grid cells, which correspond to the vector of the available environmental conditions in the study area under each climate emission scenario. This smoothing approach removes any biases to obtain a better representation of the environmental conditions suitable for each species under current and future climate emissions scenarios. This approach was implemented using a PCA which is an ordination technique that is calibrated on the entire environmental space based on the focal variables of both current and future climate emission scenarios (hereafter referred to as PCA-env). Further details can be found in Broennimann et al. (2012).
We also performed the niche similarity test, which assesses if the environmental niches of each species are more similar or more dissimilar than expected by chance in comparison to each other under current and future climate scenarios, accounting for differences in the surrounding areas of the localities (background space) under current and future climate emissions scenarios (Warren et al. 2008;Warren et al. 2010Warren et al. , 2019. The niche similarity test compares the niche overlap of one species range randomly distributed over its background, while keeping the other unchanged (e.g. Current: A. falcatus → P. elongatus), and then carries out the reciprocal comparison (e.g. Current: P. elongatus → A. falcatus). For the niche similarity test, an α value of 0.05 was considered to indicate that niches were no more similar than expected by chance.

Model evaluation indices
We calculated the commonly used (AUC) and the less commonly used (TSS and CBI) model evaluation indices. Model performance was good for all the emission scenarios (Table  1). AUC varied from 0.990 to 0.992, and TSS varied from 0.681 to 0.978. The CBI suggested that all the model predictions were consistent with the distribution of the actual data as they all had positive CBI values. However, CBI values for P. elongatus under the current emission scenario were low.

Importance of environmental variables
Precipitation of the driest quarter (BIO17) and precipitation of the wettest quarter (BIO16) were the most important variables shaping the current distribution of A. falcatus and P. latifolius, respectively (Table 2). Temperature annual range (BIO07) was the most important variable constraining the future distribution of A. falcatus and P. latifolius. The current and future distributions of P. elongatus were predicted to be constrained by precipitation of the warmest quarter (BIO18). The most influential variable shaping the current distribution of P. henkelii was predicted to be BIO18. Interestingly, the maximum temperature of the warmest month (BIO5) was the most important variable constraining the distribution of P. henkelii under RCP 4.5, while BIO18 was the most important predictor of P. henkelii under RCP 8.5.

South African podocarp distribution under climate change
The ensemble SDM predicted high habitat suitability for A. falcatus, P. elongatus, P. henkelii and P. latifolius in the exterior part of the country, near the coastal limits and excluded the interior of the country (Fig. 2). Podocarpus latifolius occupied the largest current and future geographic area followed by A. falcatus. Afrocarpus falcatus and P. latifolius were distributed in the Limpopo, Mpumalanga, KwaZulu-Natal, Eastern Cape and Western Cape provinces, which span the summer and winter rainfall regions of the country. Podocarpus elongatus was restricted to the southern parts of the Western Cape, which is associated with winter rainfall. Podocarpus henkelii was constrained to the summer rainfall provinces of Mpumalanga, the midlands of KwaZulu-Natal and the Eastern Cape, with moderately suitable habitat in the Limpopo province. Although P. elongatus and P. henkelii had the most restricted distribution, P. elongatus occupied the smallest geographic area.
Under climate emission scenarios, A. falcatus was predicted to expand to higher altitudes into the interior of the country; this was the case under RCP 4.5 and RCP 8.5 (Fig. 3). The south-western part of the distribution of P. elongatus was predicted to expand to higher elevations to track favourable conditions; however, the south-eastern and the northern parts of the distribution were predicted to become less suitable under both RCPs. Interestingly, the suitable habitat of P. henkelii under current climate was predicted to lose most of its suitable habitat under RCP 4.5 and expand its habitat to higher altitude in its northern distribution. This was the opposite for P. latifolius, where it gained suitable habitat under RCP 4.5 and lost suitable habitat under RCP 8.5. A common finding was that all four podocarp species expanded their habitat by moving to higher altitudes.

Environmental niche overlap and similarity
Podocarpus elongatus occurs over the largest environmental niche, while P. henkelii occurs over the smallest environmental niche (Fig. 4). The current environmental niches of A. falcatus, P. elongatus and P. henkelii are predicted to contract in 2070; however, the environmental niche of P. latifolius is predicted to remain stable under RCPs 4.5 and 8.5. Niche overlap between the current and future environmental niches. The current and future environmental niches of A. falcatus, P. elongatus and P. latifolius showed very limited overlap (Table 4; Fig. 4). Interestingly, P. henkelii had the highest environmental niche overlap of the four podocarps with 59.17 % and 55.43 % overlaps between the current-RCP 4.5 and current-RCP 8.5 environmental niches, respectively. When the overlap between the current and future environmental niches was further broken down, we observed that the current and future environmental niches of all four podocarps were stable (Table 4). Afrocarpus falcatus showed greater niche expansion compared to niche unfilling. A greater extent of niche unfilling was observed for P. elongatus and P. henkelii in comparison to niche expansion for both niche comparisons (Table 4). Low niche unfilling and expansion were observed for P. latifolius between the current and future environmental niches.

Discussion
This study characterized the environmental niches of South African podocarps using ordination techniques to better understand how these species may respond to the changing climate. In addition, we used ensemble SDMs to identify the environmental constraints that shape the current and future geographic distributions of A. falcatus, P. elongatus, P. henkelii and P. latifolius. We found that the ensemble SDMs accurately predicted the known geographic distributions of each podocarp species, that podocarps were generally limited by rainfall, particularly around the seed germination timeframe, and that there are likely reproductive and dispersal mechanisms as well as physiological traits that potentially prevent species from dispersing outside their current distributions.
Identifying environmental variables that shape and maintain species' geographic distributions is an essential component of evolution and ecology. Among the seven environmental variables used for these models, total rainfall during the summer and winter months (irrespective of summer/ winter rainfall region) was predicted to shape the current and future distributions of podocarps in South Africa. These findings were expected as numerous studies have identified drought stress as a factor constraining podocarp distribution and persistence (Brodribb and Hill 1998;Tesfamariam et al. 2022;Twala et al. 2022). Consequently, the distribution of podocarps in high-rainfall areas may be due to traits that reduce their ability to survive in low-rainfall regions. For instance, podocarps have small single-veined leaves that limit hydraulic conductivity (Brodribb et al. 2007) and recalcitrant seeds that become less viable as they lose moisture (Noel and Van Staden 1976;Smith 1989a, 1989b;Negash 2003). Nearly all extant podocarps exhibit low to no drought tolerance (Brodribb and Hill 2004;Twala et al. 2022) except for the New Zealand species P. totara, which uses drought-avoidance strategies (Innes and Kelly 1992) and P. macrophyllus found across China and Japan (He et al. 2020). In addition, P. latifolius recruitment is significantly reduced under drought and deep shade where seedlings are outcompeted by grasses in the landscape (Adie and Lawes 2009b). Drought and shade are antagonistic selection pressures and together with rainfall appear to constrain podocarp distribution (Chiariello 1984;Givnish 1987;Brodribb and Hill 2000;Fiorucci and Fankhauser 2017;Lusk et al. 2019). Furthermore, the models indicated that South African podocarps tend to persist in moist forest patches where fires are an important driver of forest distribution, yet fires hardly occur in the patches themselves (Geldenhuys 1994;Kellman and Meave 1997).
In addition to rainfall, the future distribution of A. falcatus and P. latifolius was predicted to be limited by temperature  annual range (BIO7). This finding suggests that in the future, there will be variation in seasonal temperatures where A. falcatus and P. latifolius occur across the highlands and midlands in the north-eastern parts of South Africa compared to the south-western parts of the country. In the north-eastern parts of South Africa, the models predict that the future daily temperature fluctuations will be smaller in comparison to present higher daily temperature fluctuations in the south-western parts of the country. Notably, the future distribution of P. elongatus was also constrained by mean temperature during the wettest quarter, which coincides with the winter months for this species as it is currently restricted to winter rainfall areas of the country. Winter rainfall and temperature are particularly important for P. elongatus because it produces recalcitrant seeds during the winter months and minimum rainfall and temperature requirements need to be met for recruitment and persistence. Podocarpus henkelii may be limited to warm high-rainfall areas because it has foliar sclereids and broad leaves (Brodribb and Hill 1998;Adie and Lawes 2011;Brodribb 2011). Furthermore, P. henkelii produces seeds that are more drought sensitive than A. falcatus seeds (Dodd et al. 1989a;Negash 1995). The current and future potential distributions of podocarps obtained from our ensemble SDMs coincide with known geographic ranges (White 1981;Farjon 2001). The broad geographically distributed P. latifolius and A. falcatus co-occur with P. elongatus and P. henkelii in their respective geographic ranges; however, P. elongatus and P. henkelii do not co-occur. Podocarpus latifolius had the broadest geographic distribution followed by A. falcatus. Both A. falcatus and P. latifolius were predicted to co-occur in the north-eastern part of South Africa, along the southern coast, and through to the southwest of South Africa. This finding coincides with the actual distribution of each species. The geographic range extends across the savanna, grassland, forest, fynbos and Albany thicket biomes, and extends to the Indian Ocean coastal belt. Afrocarpus falcatus is the only one of these four podocarps that persists in the dry forest of the coastal lowlands (Adie and Lawes 2011;Coomes and Bellingham 2011), and this occurrence is likely as a relic from population expansion during the cool and moist conditions of the last glaciation (Finch and Hill 2008;Neumann et al. 2008). Podocarpus elongatus was predicted to occupy the smallest geographic area of all the examined podocarp species under current and future climate emissions scenarios. This taxon is presently distributed in the south and south-western parts of South Africa within the fynbos and succulent karoo biomes. Podocarpus elongatus can persist in this fire-prone habitat due to its ability to resprout after fire (Midgley et al. 1995). Unlike the other podocarps, P. elongatus produces epicormic sprouts that enable tolerance of fire and persistence in the fynbos biome that is prone to fire (Midgley et al. 1995), which makes it better adapted to the Mediterranean habitat. Podocarpus elongatus may be restricted to its Gondwanan origin and its disjunct distribution as a result of retreating to higher elevations to track favourable environmental conditions such as rising sea levels during the Last Glacial Maximum that lead to loss of coastal habitat (Geldenhuys 1992). Podocarpus henkelii is predicted to have a current distribution in summer rainfall regions with high rainfall in the eastern Limpopo province, the midlands of KwaZulu-Natal and the northern Eastern Cape province of South Africa.
The effects of climate change on the geographic distribution of a wide range of tree communities have been reported across the world (Bakkenes et al. 2002;Aitken et al. 2008;Essl et al. 2011;Aguilée et al. 2016;Burley et al. 2019;Kambach et al. 2019;Pavlović et al. 2019;Sedmáková et al. 2019;Fricke et al. 2022). Podocarpus henkelii was predicted to occupy the third largest geographic range under future climate emission scenarios (Fig. 2). Under future scenarios, P. henkelii may experience a contraction of its geographic distribution under relatively 'mild' climate change scenarios (e.g. RCP 4.5) yet will potentially expand its geographic distribution as temperatures rise (e.g. under RCP 8.5, Fig. 2). Its geographic range was predicted to contract in higher elevations and expand to lower elevations for both climate change emission scenarios. The southern expansion of its geographic area under both emission scenarios may suggest that the southern parts of South Africa will become hotter under climate change (Collier et al. 2008). Therefore, it is plausible that P. henkelii is a thermophilic species and can expand its distribution with increasing temperatures under climate change. This is also evident through the reduced geographic expansion of P. henkelii under RCP 8.5 in comparison to the cooler RCP 4.5 (Fig. 3). Twala et al. (2022) reported that P. henkelii was more heat tolerant than A. falcatus. Therefore, the ability of P. henkelii seeds to disperse across the grass matrix and/or establish in newly colonized areas could limit species expansion. In contrast, under both emission scenarios, A. falcatus and P. latifolius were predicted to expand to higher altitudes and lower longitudes to track favourable climatic conditions. The dispersal of A. falcatus and P. latifolius seeds by birds, bats, baboons and monkeys (Geldenhuys 1993;Negash 2003) will likely enable them to disperse to climatically favourable areas faster relative to P. elongatus and P. henkelii. The past responses of P. elongatus and P. henkelii to climate change indicate a tight synchronization with climate fluctuations, particularly in response to water availability. For P. latifolius, recruitment is less episodic than the other species thus making it less vulnerable to climatic instability. Podocarpus elongatus was predicted to lose suitable habitat in its southern distribution under both RCPs and expand northwards, particularly under RCP 4.5. However, shifting to higher elevations and longitudinal shifts seem to be the most common responses to climate change (Pavlović et al. 2019;Sedmáková et al. 2019).
Numerous studies have reported that plant species that occupy a large geographic range also occupy different environmental niches (Cardillo et al. 2019;Kambach et al. 2019;Quiroga and Souto 2022). We found that niche conservatism was evident between current and future RCPs for A. falcatus, P. latifolius and P. henkelii (Table 3). The current and future environmental niches under RCP 4.5 for P. elongatus also showed niche similarity; however, the RCP 4.5 → current, RCP 8.5 → Current and Current → RCP 8.5 comparisons showed niche divergence. The tendencies towards divergent environmental niches in podocarps cannot be explained by differences in the climate suitability alone. However, this divergence could be explained by differences in seed dispersal, establishment requirements, habitat use, competition and ecophysiological limits. There was no consistent trend in niche conservatism or divergence among the comparisons of current environmental niches and predicted future niches across the podocarp species examined. For example, P. henkelii may shift to new geographic areas under climate change/predicted emission scenarios, per niche conservatism, given the high niche overlap value and unfilling value (i.e. the future environmental niche space is available but currently unoccupied). Notably, there may be elevated propagule pressures that ultimately limit the dispersal of P. henkelii to sites with suitable environmental conditions, which will result in a reduced realized niche. It is also possible that P. henkelii may continue to spread in its geographic range even though the breadth of its environmental niche is predicted to decrease. Relative to the other three podocarps examined, P. henkelii had the lowest niche expansion value, which suggests that it will not expand into novel environmental niches when moving southwards geographically under climate change. On the other hand, the current environmental niche of P. elongatus diverged slightly from its future environment niche due to low niche overlap as a result of high niche expansion and unfilling (Table 4). This divergence may be attributed to differences in microhabitats and competition. In addition, the current and future environmental niches of P. latifolius showed low niche overlap, which was associated with low niche unfilling values. Podocarpus latifolius can occur within areas with variable gradients of soil moisture, temperature, rainfall and altitude along different latitudes, and it produces small seeds annually which are dispersed by birds (Geldenhuys 1993;Adie and Lawes 2009b). Moreover, frequent recruitment occurs due to lower seed susceptibility to desiccation and shade-tolerant seeds, which facilitates germination and reduces susceptibility to climate change. These characteristics indicate that P. latifolius is a vagile species. Despite this, results showed low niche expansion. This may be because P. latifolius currently occupies most of the suitable environmental niches available already, and there is no additional niche to expand into. The ability of a species to occupy different environmental niches may have important implications for the understanding of the vulnerability of the species under climate change (Wiens and Graham, 2005). Previous studies recognized niche shifts in a variety of plant taxa such as Picea sitchensis (Aguilée et al. 2016), Agave spp. (Gómez-Ruiz and Lacher 2019) and Betula utilis (Hamid et al. 2019) under climate change, but the causes of these shifts are not well known. However, several studies have shown that some species will not be able to track favourable climate fast enough to escape climate change (Zhu et al. 2012;Cang et al. 2016;Burley et al. 2019). The combination of some niche expansion and unfilling for A. falcatus under both RCPs resulted in moderate niche stability under global change. The inability of A. falcatus to occupy environmental niches different from those of its current environmental niche may be due to a combination of post-dispersal mortality and loss in seed viability after harvesting, coupled with infrequent fruit production (Geldenhuys 1993;Negash 2003). Consequently, dispersal may be the main factor that governs whether A. falcatus reaches suitable habitats under global change.
Overall, we found that podocarp distribution is determined primarily by seasonal drought, likely due to the recalcitrant seeds of podocarps which are produced generally when rainfall is limited. The difference in the distribution of these podocarps underscores the variation of environmental conditions that they are adapted to and may also be a reflection on the differences in their physiology. This is particularly important as the physiological traits of podocarps may limit their distribution across environmental gradients (Brodribb andHill, 1998, 1999;Brodribb 2011;Twala et al. 2022). Physiologically oriented distribution models should be useful to predict whether the physiological thresholds of podocarps are the mechanism responsible for podocarp distribution and not their reproductive means. Insights from this study are useful for ongoing conservation actions for South African podocarps by directing the re-introduction of these species in suitable habitats and guiding the establishment of effective networks of protected areas for future dispersal events.

Data Availability
The datasets generated during and/or analysed during the current study are available in the Open Society Foundation repository, DOI 10.17605/OSF.IO/RNG5K (https://osf.io/ rng5k/).

Conflicts of Interest Statement
The corresponding author confirms on behalf of all authors that there have been no involvements that might raise the question of bias in the work reported or in the conclusions, implications or opinions stated.