UC Merced Frontiers of Biogeography Climatic variables determining Rhododendron sister taxa distributions and distributional overlaps in the Himalayas

. Endemic species in mountains are vulnerable to rapid climate change. We elucidated distributional overlaps and related climatic variables for two endemic sister taxa of Rhododendron and a generalist with respect to current and future climate conditions. Our research questions are: (i) Which climate factors separate the distributions of Rhododendron cowanianum, R. lepidotum and R. lowndesii ? (ii) How large is the geographic overlap in current and future distributions? (iii) Is it likely that the species are able to track their niches in the future? To answer these questions, we performed species distribution modelling on binomial Rhododendron occurrences accompanied by random pseudo - absences and absences constrained by other Rhododendron taxa. We used Generalized Linear Modelling to select variables, and modelled the distribution of each species using Random Forest algorithms, predicting their potential distribution in current and future climates. We also examined range differences to identify the variables segregating the distribution of these sister taxa, and estimated current and future distributional overlap between and within species. Precipitation variables explained R. lowndesii distribution, whereas temperature variables explained distributions of the other two species. We found that sister taxa have similar climate niche and hence high overlap in geographic distribution in current climate (46–68%) and potentially in future climate (53 –77%). Under future climate conditions, the potential distribution area of R. lepidotum and R. cowanianum is predicted to be at a higher elevation, while the prediction for R. lowndesii is similar to its current geography. Our models suggest that there are more potential distribution areas for these narrowly distributed endemic species than are currently occupied, which illustrates that it is rather uncertain whether the Rhododendron species will be able to track the geographical location of their niches in the future.


Introduction
Climate change affects species distributions (Parmesan and Yohe 2003), and species with a restricted distribution may be more vulnerable to the changes in climatic factors that determine the boundaries of their distributions (Thuiller et al. 2005, Manish et al. 2016. Understanding the extent to which a geographical range shift is needed for species to be able to track their climate niche in response to climate change is currently a crucial scientific task. Several organisms have already shifted their range margins (Parmesan andYohe 2003, Telwala et al. 2013). These findings are underpinned by paleo-ecological data indicating that geographical range shifts were common during previous episodes of climate change (Willis and MacDonald 2011).
The magnitude of projected climatic change is large at high elevations and high latitudes (Pachauri et al. 2014). For species in the mountains, it is easier to track their climate niche than species in flat terrains because the speed of spatial climate change is slower in the mountains (Loarie et al. 2009). In addition, many species with a restricted distribution in the mountains have relatively wide climate niches compared to species with restricted distributions in flat terrains. However, both types of species have a high risk of extinction if climate change develops novel interactions between precipitation and temperature (Williams and Jackson 2007).
In general, most species are adapted to address variations in climate, but some species that are endemic to a mountain range or a mountain peak are more vulnerable to extinction (Thuiller et al. 2005, Manish et al. 2016 and are particularly vulnerable to climate change if they have dispersal limitations (Manish et al. 2016). Mountain features may increase dispersal limitations due to steeper environmental gradients, heterogeneous microhabitats and isolation mechanisms, such as sky islands, which support a large number of unique and endemic species (Steinbauer et al. 2016). In addition, the extinction of endemic species is a global process, rather than just the loss of one metapopulation (Malcolm et al. 2006).
Although the spatial distance between different climate types is short in the mountains, it is not certain that the current combination of different climatic variables will actually exist in the future. Therefore species' survival is not guaranteed if they cannot keep pace with the climate as they move (Pearson 2006). The steady increase in mean annual temperature interacts with precipitation and the timing of the growing season, which is rather short in high mountains. The species-specific responses to warming in the moun-tains will also pose new challenges, such as competition with new species (Williams and Jackson 2007) or a lack of essential vectors for pollination or seed dispersal (Hobbs et al. 2006, Abrol 2012).
In the Himalayan region, the disappearance of current climate conditions and the development of a novel climate is expected (Williams and Jackson 2007) including an increase in the total amount and intensity of the precipitations with a reduced number of rainy days (Pendergrass and Hartmann 2014). In such conditions, the dry season becomes drier and species are found to move downwards against the direction of warming to track their precipitation niche (Crimmins et al. 2011, Qiu 2015. The species-specific responses and pace of migration may promote novel species assemblages and interactions that can lead to uncertain consequences (Hobbs et al. 2006). In this context, the conventional conservation approaches that aim to conserve representative communities or vegetation types may be ineffective (Hannah et al. 2002). This is mainly because the idea of representative communities is rooted in plant phytosociology, which assumes that the plant community responds to climate change as a unit and not as each individual species (Gleason 1926). This view of nature will be challenged by climate change (Hobbs et al. 2006, Williams andJackson 2007), and future dynamic conservation approaches will have to focus on individual species because each species may respond to the ongoing changes differently (Parmesan andYohe 2003, Telwala et al. 2013). Breshears et al. (2008) describe three possible ways of species range shifting in response to climate change; they are 'march' (defined as, range shift by colonizing leading edge, a shift in optimum and retraction at tailing edge), 'lean' (a stable range with the optimum shifting within the existing range) and 'crash' (population decline with stable edges and optimum). As such, it is important to focus on species with narrow elevational ranges and restricted geographic distributions because these specialist species will have higher risks of extinction due to their small populations and narrow ranges.
Species Distribution Models (SDMs) are be-ing used to predict potential spatial and temporal distribution of species (Thuiller et al. 2005, Randin et al. 2006 although their relative success when transferred to future conditions is at stake (Araújo and Rahbek 2006). Species distribution shifts are mostly studied within a single taxon, between sister or descendent taxa and within communities (Thuiller et al. 2005, Mao andWang 2011). Sister taxa are assumed to have common ancestors and are therefore expected to show some degree of niche overlap because niches are, to some extent, conserved within a clade (Wiens andGraham 2005, Losos 2008) while maintaining some distinctions among themselves (Cavender- Bares et al. 2004). In the Himalayan region, studies on niches, distribution overlaps and shifts of sister taxa are rare (but see Vetaas 2002). We address this gap by studying Rhododendron sister taxa from the central Himalayas. The target sister taxa belong to the subgenus Hymenanthes, subsection Lepidota. One species has a wide distribution from the western and the eastern Himalayas to China (R. lepidotum Wall), whereas the two other species have restricted distributions in the central Himalayas (Nepal: R. cowanianum Davidian and R. lowndesii Davidian).
Here we seek to address: (i) Which climate factors separate the distributions of closely related R. cowanianum, R. lepidotum and R. lowndesii species? (ii) How large are the current geographic overlaps between them and what will their potential overlap under future climatic conditions? (iii) Is it likely that the species are able to track their niches in new geographic areas?

Study area
This study was carried out within the distribution range of the genus Rhododendron across Nepal in the central Himalayas.

Taxa
Members of the genus Rhododendron L. (Ericaceae) are phanerophytes, i.e., shrubs or medium-sized trees. Rhododendron has a wide temperature range, from warm temperate zones to alpine bioclimatic zones. There are 43 lower taxa of Rhododendron in Nepal between approximately 900 m and 5600 m above sea level (a.s.l.) (www.efloras.org). The Lepidota (Hutchinson) Sleumer subsection of the genus Rhododendron includes three sister taxa, R. lepidotum Wall, R. cowanianum Davidian, and R. lowndesii Davidian, which are distributed between approximately 2100-4700 m a.s.l. in Nepal 1 . Among these, the latter two are rare and endemic to Nepal (Rajbhandari et al. 2016).

Occurrence data
We compiled occurrence data from the National Herbarium and Plant Laboratories, Kathmandu, Nepal, the Royal Botanic Gardens Edinburgh, UK, the Tokyo Herbarium, Japan, the Global Biodiversity Information Facility 2 , and from field sampling. Initially, we recorded 25, 420 and 46 presence data for R. cowanianum, R. lepidotum and R. lowndesii respectively. Among the collected occurrence points, we filtered out some points with high uncertainty. First, we excluded points with a very crude accuracy of the location, i.e. latitudinal and longitudinal values with less than three digits after decimal place (number of points removed: 6 points for R. cowanianum, 37 for R. lepidotum, and 6 for R. lowndesii). Secondly, we omitted specimens with elevation below 900 m and above 5600 m a.s.l. as they were more than 1000 m below or above the lowest and highest record of the Rhododendron species concerned 1 . This yielded 19 presence points for R. cowanianum, 271 for R. lepidotum and 40 for R. lowndesii.

Pseudo-absence data
Most SDMs and niche models are based on presence-absence data. However, species data are mostly composed only of recorded presences. In such cases, absences are complemented by pseudo-absence data for environmental information (Elith et al. 2011). There is no consensus on how to generate the best pseudo-absence data, and most studies use the random pseudo-absence method, which is equal to or better than other methods (Barbet-Massin et al. 2012). We used two different methods to generate 'pseudoabsences' to test which one would perform better.
The first approach was to use the presence points for all Rhododendron species in Nepal except the target species as absence points combined with the presence of the target species (hereafter; "Rhododendron pseudo-absences" = "RhoPs"). The approach constrains the pseudo-absence points to be within the climatic envelope of the genus, thereby avoiding "naughty noughts" placed far outside the potential climate range (Austin and Meyers 1996). This kind of pseudo-absences has been used for Eucalyptus in Australia and Rhododendron in Nepal (e.g., Austin et al. 1990, Vetaas 2002. Among the collected occurrence points, we filtered out some points with high uncertainty using the two-step filter described in the previous section. With this method we obtained 890 Rhododendron pseudo-absences. The second approach was to use randomly generated equal numbers of pseudo-absences combined with presence data (hereafter; "Random pseudo-absences" = "RanPs") within the same elevational range.

Predictor variables
We used 22 water and energy related predictor variables, including 19 bioclimatic variables from the WorldClim 3 dataset (Hijmans et al. 2005), Annual BioTemperature (ABT; Holdridge 1947), Ellenberg's Climatic Quotient (EQ; Ellenberg 1963) and the Relative Radiation Index (RRI; Oke 1987). All climatic data required for preparing the ABT and EQ were taken from the WorldClim dataset (method details in Supplementary Material S1). All predictor variables were in a 30 arc-second resolution and the same coordinate system (WGS 1984), and can be made available upon request to the authors.
We prepared two groups of variables from the original set of 22. The first group was composed of all variables (hereafter the "set I" variables) and the second group was prepared by selecting a few effective variables from a Generalized Linear Model (GLM) using the bidirectional (forward and backward) selection method in R package stats (R Core Team 2016). Then, we dropped the non-significant variables. For RanPs this yielded 9, 10 and 6 variables for R. cowanianum, R. lepidotum and R. lowndesii, respectively, and for RhoPs 12, 10 and 7 (Supplementary Material S2), (hereafter the "set II" variables). The optimum GLM models (set II variables) were partitioned to obtain the deviance explained by temperature and precipitation related variables using the R-package ecospat (Broennimann et al. 2016).

Future climatic scenario
We used the Intergovernmental Panel on Climate Change's (IPCC) most extreme future prediction (worst-case scenario), Representative Concentration Pathway 8.5 (RCP8.5) for our future climatic scenario because when we look at last few years, it is hard to be optimistic that the world's countries will succeed in limiting the warming to 2°C by the end of the 21 st century (UNFCCC 2015), especially as recent monthly mean temperatures and annual mean temperatures have broken previous records (GISTEMP Team 2016). The RCP8.5 projects 2.6°C to 4.8°C warming by 2081 to 2100 compared to the 1986 to 2005 baseline (Collins et al. 2013). We took the average of five different downscaled General Circulation Models, namely the ACCESS1-0, BCC-CSM1-1, GISS -E2-R, MIROC-ESM-CHEM, and MPI-ESM-LR models, to reduce model-derived biases. We predicted our results for only one worst-case scenario and for a single future period in the 2070s (average of 2060 to 2080).
The values of the predictor variables that were in raster format were extracted to the presence, rhododendron pseudo-absence, random pseudo-absence and lattice files (regular grid points of 3 km resolution above 900 m a.s.l.) for current climate and future climate in ArcGIS 10.3 (ESRI).

Distribution modelling and variable range difference analyses
To answer the first research question, which climatic factors segregate the closely related three Rhododendron sister taxa, Tukey's Honesty Significant Difference (HSD) post hoc test was used to identify the difference in range for all 22 variables for each species using R package stats (R Core Team 2016). Species distribution models were prepared to predict the potential distribution of species in current and future climate using the Random Forest method (Breiman 2001). The predictions were portrayed into geographic space to analyse the overlaps between species. The Random Forest method was used among different techniques because it can handle multiple variables regardless of their eventual multicollinearity, low numbers of presence points and different prevalence ratios (Elith et al. 2011, Barbet-Massin et al. 2012). All analyses were performed in the R package RandomForest (Liaw and Wiener 2002).
In the Random Forest method, we fitted models on RanPs and RhoPs with the set I and set II variables. The datasets were partitioned at 3:7 ratios for test and training datasets. We grew 2000 trees, as growth appeared to stabilize by 1000 -1500 trees. The model was replicated five times. Each time, we evaluated the Area Under the Curve (AUC) of the Receiver Operating Characteristic (ROC) value. Important variables are listed based on their Mean Decrease Gini in Random Forest. Predictions of the relative index of occurrence (RIO) of species were made from each replicate of models on current climate and future climate lattice files. Then average RIO was calculated from five predictions. The predicted RIO value ranges from 0 to 1; where a higher value refers to more suitability of the location. At the end, we had a total of 24 different predictions. Then, RIO raster maps for each species were prepared by interpolating the average RIO using the Inverse Distance Weighted tool in ArcGIS 10.3 (ESRI). This raster was converted to ASCII format to be fed into the distributional overlap study.

Distribution overlap analysis
To answer the second research question, we studied the predicted distribution overlap between species using ENMTools (Warren et al. 2008) with three different available methods, including Schoener's D, I statistics and Relative Rank (RR), for both current and future climates. Then, we compared the predicted distributional overlaps based on the average of the three methods. The value ranges from 0 (no overlap at all) to 1 (complete overlap).

Geographic shift of climatic niche
To analyse the geographic shift of climatic niches of three Rhododendron species, i.e. the third research question, the predicted average RIO values of the lattice points were plotted against elevation for the current and future projected climate for each species, and the shift was analysed graphically as it could not be quantified because we did not convert the RIO into a binary value. We plotted the points with RIO above or equal to 0.02 for better illustration.

Effects of environmental dimension reduction analysis
The models with the set I and set II variables with RanPs and RhoPs were compared based on the AUC value and ROC curve plots in order to figure out the effect of dimension reduction in models. Then, the differences between their predictions were tested in ENMTools. In this analysis, the value ranges from 0 (complete dissimilarity) to 1 (total identical).

Variables segregating species distributions
All but two temperature variables had similar ranges for the three species. Out of 22 variables, R. lepidotum had five variables' ranges that were significantly different from R. lowndesii, while the ranges of six variables were significantly different between R. lepidotum and R. cowanianum, and the ranges of five variables were significantly different between R. lowndesii and R. cowanianum ( Fig. 2 and Supplementary Material S3). Mean temperature of the driest quarter (bio09) and precipitation of the driest quarter (bio17) were significantly different between R. lowndesii and both of the other species (Fig. 2). A twodimensional niche plot of these variables (Fig. 3), showed a higher overlap of the generalist R. lepidotum and both endemic sister taxa, while R. lowndesii and R. cowanianum had smaller overlaps. Climate variables that were significantly different between the realized distributions of the three species. Variable acronyms correspond to isothermality (bio03), mean temperature of the driest quarter (bio09), annual precipitation (bio12), precipitation of the wettest month (bio13), precipitation of the driest month (bio14), precipitation seasonality (coefficient of variation) (bio15), precipitation of the wettest quarter (bio16), precipitation of the driest quarter (bio17), precipitation of the warmest quarter (bio18), precipitation of the coldest quarter (bio19) and Ellenberg Climatic Quotient (EQ). Based on the Random Forest models of RanPs (results are not illustrated from RhoPs models as they were consistently poor, details below), the most important variables (based on mean decrease in Gini index) in both sets I and II that explained the distribution of R. lowndesii were precipitation of the wettest month (bio13) and precipitation of the warmest quarter (bio18). In the case of R. lepidotum, in set I the most important variables were isothermality (bio03) and precipitation of the coldest quarter (bio19) and in set II they were isothermality (bio03) and mean temperature of the wettest quarter (bio08). The distribution of R. cowanianum was mostly explained in set I by the Ellenberg Climatic Quotient and precipitation of the driest quarter (bio17) and in set II by precipitation seasonality (bio15) and Annual BioTemperature. Although there were differences in the most important variables among sister taxa (Table 1), there were only a few variable ranges that were significantly different between them (Fig. 2).
The variance partitioning analysis showed that the deviances explained by temperature and precipitation are 30% and 22.9% respectively for R. cowanianum in the optimal GLM (set II) for RanPs (the total explained deviance was 66.3%). For R. lepidotum, the total deviance explained was 51.9%, of which 37.1% was explained by temperature variables and only 6.3% by precipitation variables. For R. lowndesii, the deviance explained by precipitation variables was 66.5%, which was tenfold higher than the deviance explained by temperature-related variables (5.7%), and the total deviance explained was 64.0%.

Distribution overlaps under current and projected future climate
The distributional overlap analysis verified a high degree of distributional overlaps ( Variable acronyms stand for: annual mean temperature (bio01), isothermality (bio03), minimum temperature of the coldest month (bio06), mean temperature of the wettest quarter (bio08), mean temperature of the driest quarter (bio09), mean temperature of the warmest quarter (bio10), annual precipitation (bio12), precipitation of the wettest month (bio13), precipitation seasonality (bio15), precipitation of the wettest quarter (bio16), precipitation of the driest quarter (bio17), precipitation of the warmest quarter (bio18), precipitation of the coldest quarter (bio19)

Geographical shifts of climatic niche
The predictions of our Random Forest models predictions using RanPs suggested that climatic niches of R. cowanianum and R. lepidotum will move to higher elevations with projected warming. However, the climatic niche of R. lowndesii does not seem to move uphill in future climate projections (Fig. 5). The results were consistent across both sets of variables.

Effects of reducing environmental dimension on species distribution models
There were some differences in the predictions using set I and set II variables. The similarities are depicted in Table 3. The respective AUC values of the Random Forest models are also illustrated in the table. The ROC curves for the set I and set II variables were also close to each other (Supplementary Material S6).

Discussion
In this study, we verified that three closely related Rhododendron sister taxa have similar relationships to most climatic variables. As these taxa are phylogenetically highly related and geographically very close in the Himalayas, their distributions partially overlap. The distribution model suggests that the potential area of distribution of species adapted to arid environments will not move to higher elevations, whereas the potential area of distribution of the other two sister species will move to higher elevations in the future climate.
Here, the potential area of distribution shift in geography is based on the 'worst-case' climatic scenario.

Climatic factors segregating species distribution
Based on the Random Forest model with the set I and set II variables on RanPs, precipitation of the wettest month and the warmest quarter are the most influential variables for the distribution of R. lowndesii. This aligns with empirical data that this species is mainly observed in dry regions in Nepal, where any amount of precipitation is important. The R. lepidotum distribution is mostly related to isothermality, a measure of how variable is temperature within each cell derived from diurnal and annual temperature ranges. Distribution of R. cowanianum is related to both temperature and precipitation variables as the most important variables. In other words, for R. lepidotum and R. cowanianum, water is not the main limiting factor in their distributions, in contrast with R. lowndesii (Fig. 3). This finding agrees with Cavender-Bares et   Table 3. Similarity in the potential area of distribution (in percentage) between models with all variables (set I) and those with GLM-selected variables (set II) under current and future climate conditions and their respective model AUC values.
al. (2004), who found that phylogenetically close oak species share contrasting moisture preferences in North Central Florida. The most important variable lists differ between sister taxa ( Table 1) and most of the variables' ranges are similar between them (Supplementary Material S3), which supports previous findings that sister taxa possess similar climatic niches on a broad scale (Hof et al. 2010) and indicates the conservation of phylogenetic niches (Losos 2008).

Range shifts and distribution overlaps under current and projected future climate
The results of Tukey's HSD tests suggest that the highest distributional overlap is found between the generalist species R. lepidotum and the two endemic sister taxa. On average, the sister taxa have approximately 58% (set I) and 54% (set II) overlaps in their geographical distribution (Supplementary Material S5A and S5B). This over-lap is higher than the one found between descendent and parent species in the Tibetan Plateau estimated by Mao and Wang (2011). They found 32% to 36% overlap between Pinus densata, a descendent from the hybridization of its parent species P. tabuliformis and P. yunnanensis. However, the distributional overlap between the three Rhododendron species was smaller than the 80% of distributional overlap between sister taxa found on a study with 71 different species in the California Floristic Province (Anacker and Strauss 2014). We estimated potential geographic distribution overlaps between current and future climates, assuming that the species may be able to track the geographical location of their niche, but many factors such as soil conditions, vectors for pollination, and dispersal may hamper a potential shift in geographical location therefore projected changes are always rather uncertain (Parmesan and Yohe 2003, Rahbek 2006, Svenning et al. 2010). The degree of distributional overlap under future climate conditions is predicted to be almost the same between R. cowanianum and R. lepidotum, whereas it may increase between R. cowanianum and R. lowndesii, while the overlap between R. lepidotum and R. lowndesii is predicted to be slightly lower by set I variables and slightly higher by set II variables. This prediction agrees with the assumed niche conservatism within sister taxa (Wiens and Graham 2005). Within species, changes in the distribution of approximately 30% (set I) and 26-40% (set II) are predicted between current and future climate conditions (Supplementary Material S5A and S5B). Based on the predictions, to be able to track their current niches R. lepidotum will have to 'march', and R. cowanianum will have to 'lean' and 'march' (Fig. 5). These species may move upslope with predicted warming as seen in other Himalayan species (Telwala et al. 2013). However, it is not necessarily true that all species require shifting upslope with warming (Crimmins et al. 2011, Qiu 2015; for instance, the potential area for R. lowndesii in the future climate is predicted around its current elevation. This is explained by precipitation. In the Himalayan region, the amount of precipitation has an inverse relationship with elevation, moreover, the future precipitation is predicted to be less frequent (Pendergrass and Hartmann 2014), which means that dry areas will be drier. In this situation, species may tend to stay behind the temperature niche or move downhill to track their precipitation niche. Similar instances are reported by Crimmins et al. (2011) in California, USA andQiu (2015) in southern Tibet, China. This shows that geographical shifts along mountainsides are species-specific and more complex than just upward shifts (Gleason 1926, Halpin 1997. Climate change may be a real threat to some endemic species if they fail to migrate due to dispersal limitations or if lack of adequate soil conditions prevent them from establishing in a new geographical location even if it is within their climate niche (Thuiller et al. 2005, Pearson 2006, Manish et al. 2016). This will in essence create large challenges for contemporary strategic biodiversity conservation (Hannah et al. 2002). Moreover, species-specific geographic shift rates (Parmesan and Yohe 2003) may involve the emergence of new community assemblages leading to novel ecosystems under future climate conditions (Hobbs et al. 2006, Williams andJackson 2007). In this context, contemporary conservation practices may have to change from ecosystem and/or community oriented to individual species oriented (e.g., red-listed species) because conventional strategies for communities may not be suitable for rare and endemic species in a dynamic future context. Hence, conservation strategies should incorporate climate change and focus on mountains when selecting protected areas in the future (Araújo et al. 2004).

Effects of reducing environmental dimension on distribution models
Here, our strategy of dimensionality reduction provided good results. In general, the AUC has a positive relationship with the number of predictor variables (Synes and Osborne 2011). In contrast, we found a negative relationship in R. lepidotum. We found that the model performances with set I (including all the environmental variables) and set II (reduced set) variables are very close to each other when the prevalence ratio is higher, with low differences between the predictions. However, set I is better at low prevalence ratios. This suggests that the model can be simplified by reducing the number of predictor variables. Here, we separately selected variables for three species using GLM, which is a recognized method for selecting effective variables (Guisan et al. 2002), and generated different combinations of variables (Supplementary Material S2).
In this study, the prevalence ratio was not equal among species, as rare species had a low number of presence records. The lower number of occurrences for rare species can hinder statistical analysis. However, other studies have shown that such low occurrences of rare species data are acceptable and more accurate predictive models can be developed for rare and restricted range species (Franklin et al. 2009). The narrow environmental range and restricted geographic distribution may have enabled the SDM to predict with higher accuracy for endemic and rare species despite the low number of occurrences. Here, our results support the previous findings. We found a higher AUC value for both the rare and endemic species (R. cowanianum and R. lowndesii) compared to the generalist species R. lepidotum (Table 3).
The prediction accuracy and model performance measures do not only depend on the number of presences, but are also affected by the number of pseudo-absences (VanDerWal et al. 2009). Here, we tested models wherein the numbers of pseudo-absences were set equal to the number of presences (results not included here). There are many different ways to distribute the 19 pseudo-absence points for R. cowanianum in the study area. We found that when the pseudoabsences were at a distance from the presence locations, the AUC was higher and the prediction was better than when the pseudo-absence points were close to the presence locations, which agrees with VanDerWal et al. (2009). This is why the RanPs models always outperformed the RhoPs models. This finding reveals that sister taxaconstrained absence values are not better than randomly generated pseudo-absences. This result is consistent with a finding by Barbet-Massin et al. (2012). The reason behind the poor performance of the sister taxa-constrained absence value is because of a low discrimination power within the model between the targeted presences and the constrained pseudo-absences as they are both within close proximity.
In conclusion, our models suggest that there is high climate niche overlap and thereby high geographical overlap for the sister species, but there are also more potential geographical areas for the two endemic species not occupied, which may relate to dispersal limitation or other environmental factors. The modes indicate that R. lepidotum will have to 'march', and R. cowanianum will have to 'lean' and 'march' to track their future climate niche, whereas R. lowndesii may stay behind, because its distribution is determined by precipitation. This illustrates that responses to climate change are very individual and it is also rather uncertain whether the Rhododendron species are able to track the geographical location of their niches in the future.