Big data help to define climate change challenges for the typical Mediterranean species Cistus ladanifer L.

1. Introduction

The impact of the going on climate and land use changes on plants, species interactions, and ecosystem functioning was yet unpredictable. Climate impacts on species were estimated to be severe, considering changes predictable in near-surface temperature, precipitation, frequency and intensity of heat waves, forest fires, and pest outbreaks. In the Mediterranean region, an increase in seasonality is also expected, as well as in late frost, heat waves, intensive summer drought, punctual heavy rains, and fire risks (Alcamo et al., 2007). Therefore, bioclimatic envelopes are expected to shift considerably within one century (Loarie et al., 2009). The future climate change effect on plant ecosystems is predicted to affect the ability of plants to adapt to stresses, particularly in the Mediterranean region. The warming up is 20% faster than the global average, and land-use fragmentation might result in species’ collapse, considering that the plants’ adaptation to fast-changing climatic conditions is also unknown. How quickly plant species can adjust to global change, despite sizeable genetic variation detected in outcrossed species (Gomulkiewicz and Houle, 2009), is a question yet to be answered. The populations might go extinct before adapting to the new environment if genetic modifications are too slow (Jump et al., 2006). Climate change impact on well-adapted Mediterranean species was rarely studied considering the full range of their distribution (Serra-Varela et al., 2015; Viruel et al., 2020; Almeida et al., 2022), and including genetic and palynological data. The current study aimed to forecast this phenomenon in a species that is particularly adapted to low precipitation and high-temperature summers, the sclerophyllous shrub Cistus ladanifer L. (rock rose). Indeed, an increased risk of desertification and the need for biodiversity conservation in this region, a biodiversity hotspot, might be at stake, even for drought-resistant species.

The complex history of the Mediterranean flora has been influenced by three main factors: geological history, climate, and human activities (Thompson, 2005). This region is characterized by several mountain ranges taken as physical dispersal barriers for lowland species, but also as significant refugia during Pleistocene glaciations [starting at approximately 2.6 million years (Ma)] (Médail and Diadema, 2009). These extreme climatic events were decisive for the current floristic composition of the Mediterranean Basin, just as were the Messinian Salinity Crisis (6–5 Ma) or the intensification of climate seasonality during the Middle Miocene (16–14 Ma) (Thompson, 2005; Rundel et al., 2016). Through these periods, tropical and subtropical elements were progressively replaced by the drought-tolerant woody sclerophyllous plants that dominate the region at present (Suc et al., 2018).

The Mediterranean region is the third-world plant biodiversity hotspot (Buira et al., 2021, and references therein). Despite having less species density than other Mediterranean-climate regions (Cowling et al., 1996), it has the highest number of species compared to the other regions with similar climate, and harbors a global hotspot for recent rapid speciation (Valente and Vargas, 2013); being the adaptive radiation of Cistus (Cistaceae) a representative example (Guzmán et al., 2009). The Cistus diversification appears to have been prompted by the onset of the Mediterranean climate (Fernández-Mazuecos and Vargas, 2010), whose early stage took place, probably, in the western Mediterranean (Guzmán and Vargas, 2005) during the warmest periods of the Pleistocene (Thompson, 2005). Cistus ladanifer is one of the most recent species in the genus, which diverged into three subspecies (Demoly and Montserrat, 1993): africanus [the basal one, occurring in the north of Africa and the southern tip of the Iberian Peninsula (IP)], sulcatus (restricted to the southwestern coast of Portugal), and ladanifer (the most widespread subspecies, spanning the west and center of the IP, south of France and north of Africa).

Cistus ladanifer is a perennial sclerophyllous scrub with resin, due to labdanum production. The species occurs mainly in oak or pine forests understory and tends to dominate disturbed and low-biodiversity scrublands. The populations constitute early successional stages adapted to disturbances operating in Mediterranean ecosystems, particularly fire. Moreover, this species have high production of seeds with considerable longevity, leading to persistent seed bank formation, and seeds have high physical dormancy, broken by high temperatures generated by a fire in soil top layers (Bastida and Talavera, 2002, and references therein). Seed dispersion is mainly local due to barochorous dispersion (Bastida and Talavera, 2002), though long dispersal is possible by animals that feed on seeds and fruits (Malo and Suárez, 1998). Additionally, seed production with different degrees of dormancy might help the species’ fitness in harsh environments, with germination happening after a fire but harder seeds germinating afterward (Delgado et al., 2001). Cistus ladanifer is often perceived as a prejudicial plant, though with valuable applications in the pharmacological and aromatic industries due to its resin, the labdanum (Quintela-Sabarís et al., 2011; Frazão et al., 2018), and creating associations with other species of commercial interest (Alonso Ponce et al., 2011). Hence, increasing the scientific knowledge on this vast and valuable resource is crucial to improve current and future species management and to capture information about the predictable influence of climate change on very well dry-adjusted species.

Currently, the C. ladanifer native area occupies a large part of the IP, especially the occidental and central regions, but it is also present in southern France, northern Morocco, and Algeria, mostly in poor and acidic soils and, infrequently, in ultramafic and calcareous soils (Demoly and Montserrat, 1993). This species is well adapted to hot and dry landscapes and is expected to be highly resilient to climate change, which is starting to increase the aridity of the Mediterranean Basin (Guiot and Cramer, 2016). Nonetheless, doubts remain about the subspecies’ sulcatus future resilience since the distribution range will decrease considerably over the current century (Ferreira et al., 2021).

A phylogeographic analysis for C. ladanifer suggests the existence of glaciation refugia in southern Iberia and northern Africa, enclosed by the Rif and Betic mountain ranges, in North Africa and South of Spain (Quintela-Sabarís et al., 2011). This resulted, besides climatic oscillations, as a partial consequence of human activities, which have been historically prevalent in the Mediterranean Basin (Thompson, 2005), since this species is an avid colonizer of degraded lands. Wildfire intensification (Pérez-García, 1997) and land abandonment (Mendes et al., 2015) are both beneficial for the emergence of C. ladanifer scrublands, whose allelopathic behavior hinders the occurrence of other plant species (Sosa et al., 2010; Gallego et al., 2020). Apart from some indications of potential glacial refugia and dispersal routes (Guzmán and Vargas, 2009; Quintela-Sabarís et al., 2011), there was little information concerning C. ladanifer historical range. Improving the knowledge about this species’ past biogeography and how it will respond to climate change across the entire native range is needed. Species distribution models (SDMs) are increasingly employed to unveil those biogeographic patterns (Svenning et al., 2011; Wróblewska and Mirski, 2018).

The general objective was to investigate past species dynamics using molecular and palynological analyses and SDM modeling information, as the past data might help predict future range changes. However, future changes are foreseen to happen at a much faster pace. With the example of a well-adapted species, the climate change impact on the Mediterranean flora could be anticipated. A claim to policymakers for ecosystem management and maintenance could be targeted to minimize near-future predictable desertification in this region. Additionally, the rock rose has potential economic value due to the labdanum production, a resin from where valuable essential oils are extracted, with medical, cosmetic, and food industry putative applications (Frazão et al., 2018, and references therein) and improvement programs could be designed due to the gathered information.

Employing species distribution models methodology, this study aimed at (i) predicting the potential current C. ladanifer distribution in its native range area, (ii) identifying the most suitable areas and potential glacial refugia for this species in three past periods based on SDM modeling, phylogeography, and fossil data, and (iii) projecting the model into future periods, analyzing the effect of climate change on its current distribution.

2. Materials and methods

2.1. Study area

The study area included the western part of the Mediterranean Basin (latitude from 29.88°N to 49.05°N and longitude from 13.18°W to 13.00°E), from Portugal to southern Germany and western Italy, in Europe, and, from Morocco to northwestern Libya, in Africa (Figure 1). The Mediterranean climate prevails in the study area, with harsh summers, generally high temperatures and low levels of precipitation and climate heterogeneity, a biodiversity trigger; thus, this region includes one of the most remarkable non-tropical biodiversity hotspots (Myers et al., 2000).

FIGURE 1

Figure 1. Location of the study area and spatial representation of the occurrence dataset. Light green dots represent the species’ presence (see Supplementary Table S1 for details). Map of the pollen sites in the Mediterranean region with Cistus ladanifer and C. ladanifer Type records. The map was built with ArcGIS® Desktop v10.8.1.

2.2. Occurrence data

The species’ presences were collected from different data sources (Supplementary Table S1). The dataset was cleaned to retain high-quality distribution records: some were dubious or had coordinates with large uncertainty values. All the records located outside the study area or in regions where C. ladanifer is not autochthonous were eliminated. Thus, the data found in Macaronesia, the Mediterranean islands, the northern half of France, where the Mediterranean biome is not predominant, and the Moroccan Atlantic coast south of Larache (Mohammed Ater, personal communication) were removed.

The cleaned dataset, including 560,915 records, was further reduced to one occurrence per grid cell of ca. 5 km² to attenuate spatial autocorrelation, resulting in a final set of 2,833 records (Figure 1), to avoid the adverse effects of the spatial model bias. The selection of a single record per grid cell was sufficient, as already demonstrated elsewhere (Gorelick et al., 2017).

2.3. Environmental data

Nineteen bioclimatic variables (BIO1 to BIO19), with a spatial resolution of 30-s (~1 km²) (Supplementary Table S2), were downloaded from the Worldclim v1.4 (Hijmans et al., 2005), a database of high spatial resolution global weather and climate data. They were downloaded for the past, the Mid-Holocene [MH: ~6,000 BP (years before present)], the Last Interglacial (LIG: ~120,000–140,000 BP), the present, and the future. The future climate data comprised the years 2050 and 2070 and two Representative Concentration Pathways (RCP 4.5 and RCP 8.5). The RCP 4.5 is a moderate emission scenario for the stabilization of radiative forcing by 2,100 (Thomson et al., 2011), and the RCP 8.5 is the high-emissions scenario, called “business as usual” (van Vuuren et al., 2011), the worst-case scenario.

The CHELSA v1.2 (Karger et al., 2017) was chosen to download the bioclimatic variables for the Last Glacial Maximum (LGM ~22,000 BP), given the 30-s resolution unavailability for variables for this period in the Worldclim v1.4. The Global Climate Models (GCMs) CCSM4 (Community Climate System Model 4.0) and the MRI-CGCM3 (Meteorological Research Institute-Coupled Global Climate Model 3.0) were selected for all the periods under study, but the LIG, since the extracted information was only available, for this period, according to the displayed source (Otto-Bliesner et al., 2013).

The extracted variables were narrowed to an adequate subset of non-collinear variables using the function removeCollinearity from R package virtualspecies v1.5.1 (Leroy et al., 2016). The collinearity among the 19 bioclimatic variables was verified using Pearson’s correlation coefficient (Dormann et al., 2013), and the multicollinearity cutoff = 0.8 was used to select the subset of non-collinear variables.

The final set included six variables, half of them temperature-related (BIO1 – Annual Mean Temperature, BIO5 – Maximum Temperature of Warmest Month, BIO10 – Mean Temperature of Warmest Quarter), and another half precipitation-related (BIO14 – Precipitation of Driest Month, BIO17 – Precipitation of Driest Quarter, and BIO18 – Precipitation of Warmest Quarter). The summary statistics considered for each selected variable, the GCMs, and across all the considered periods, past, current, and future, were shown in Supplementary Figure S1.

2.4. Species distribution modeling

The current species distribution models, including the ensemble, were implemented in R (R-Core-Team, 2021) with the package biomod2 v3.5.1 (Thuiller et al., 2021), using the species occurrences, the six environmental variables previously selected, and four different algorithms. The ensemble model was further projected into the past (LIG, LGM, and MH) and future (RCP 4.5 2050 and 2070; RCP 8.5 2050 and 2070) scenarios. The algorithms involved in the study included two regression methods and two machine-learning methods. The regression-based methods were the Generalized Linear Model (GLM) (McCullagh and Nelder, 1983) and the Multivariate Adaptive Regression Splines (MARS) (Friedman, 1991), and the machine-learning ones were the MaxEnt (Phillips et al., 2006) and the Artificial Neural Networks (ANN) (Hopfield, 1982). The biomod2 default settings were applied for all the models, according to the typical usage of biomod2 (Hao et al., 2019). The resulting ensemble model based on different methods was suggested in the literature as a proxy for predicting species range due to the considerable variation in different model types’ performance and consensus inexistence on the best technique. Thus, individual models’ averaging results may increase the overall accuracy of predictions, and compromise among different algorithms could help constrain prediction error (e. g., Jones et al., 2010; Ochoa-Ochoa et al., 2016). Moreover, ensembles tend to perform better than most single models by quantifying model congruence, recognizing the highest uncertainty regions, and moderating inaccurate predictions risk (Stewart et al., 2022).

Three sets of 100,000 random pseudo-absences (or background) were randomly generated to guarantee that environmental variation in the study area was considered. Five-fold cross-validation was performed to test the algorithms’ accuracy by splitting presence records into training (75%) and test datasets (25%). Summing up, 15 models (3 plus 5) were executed in the biomod2 using the four algorithms, resulting in a total of 60 models (see Supplementary Figure S2 for a graphical overview of the methods).

The obtained models’ performance was evaluated using two criteria, (i) the Area Under the Curve (AUC) of the Receiving Operator Characteristics (ROC) and (ii) the True Skill Statistic (TSS), using the BIOMOD_Modeling function, and the argument – models.eval.meth. The AUC ranges from 0 to 1, with higher values corresponding to better models’ prediction. The AUC has the advantage of being both prevalence and threshold-independent. The TSS is also prevalence independent and ranges from −1 to 1, where 1 indicates perfect agreement and inferior to 0 indicates a performance no better than random (Swets, 1988; Fielding and Bell, 1997; Allouche et al., 2006; Supplementary Figure S2). Moreover, the importance of the variables for each model was calculated using the variables_importance function. This function returns the variable importance value for each variable in the model. The highest the value, the more influence the variable has on the model.

The ensemble model was obtained using the BIOMOD_EnsembleModeling function, and the algorithms produced models with a TSS value equal to or above 0.6 (Jones et al., 2010), considering the mean probabilities across predictions. This ensemble model was further projected into the future and the past using the BIOMOD_EnsembleForecasting function. For each period (except the LIG), the results of the two GCMs were averaged to reduce uncertainty.

The potential distribution maps obtained for all scenarios, with values between 0 and 1 for each pixel, were reclassified into five suitability categories using the equal interval approach: (0.0–0.2; 0.2–0.4; 0.4–0.6; 0.6–0.8; 0.8–1.0). To predict changes, expansions, and contractions in the species’ range from the past to the present and from the present to the future distributions, the ensemble model and each projection were converted to binary presence-absence maps of estimated presences (1 – suitable habitat) and absences (0 – unsuitable habitat) through the maximization of the TSS (Thuiller et al., 2021). The ArcGIS Combine tool was used to compare differences between pairs of binary maps to identify potential shifts and trends in suitable habitats. This function allowed us to determine if the C. ladanifer distribution area increased | habitat gained (0 → 1), decreased | habitat loss (1 → 0), remained suitable area | habitat maintained (1 → 1), and remained unsuitable area | unsuitable habitat (0 → 0).

2.5. Palaeobotanical and genetic data

The models projected in the past were integrated with the paleobotanical and genetic data to check their plausibility. Twenty-four sites corresponding to fossil pollen records were found in the literature and the Neotoma Paleoecology Database (Williams et al., 2018), further compiled in Supplementary Table S3 and mapped in Figure 1.

The Quintela-Sabarís et al. (2011, 2012) studies were used as a source of genetic data. In the first study, the authors gathered results from 38 C. ladanifer populations, covering almost the entire natural range of the species, including the three subspecies, screened for length variation of polymorphic chloroplast simple sequence repeats, (cpSSRs), to investigate the species phylogeographic structure to locate Quaternary refugia, to reconstruct its recolonization patterns and to assess the role of geographical features (mountain ranges, rivers and the Strait of Gibraltar) as barriers to seed flow and expansion through the Western Mediterranean. In the second study, nuclear markers were used (AFLPs) and 33 populations were screened. Four clusters were unveiled by using Bayesian analysis [Figure 3 in Quintela-Sabarís et al. (2012)].

3. Results

3.1. Variable importance and model performance

The precipitation variables BIO18 (1), BIO14 (2), and BIO17 (3) had decreasing importance in the ensemble model, followed by the temperature ones BIO10 (4), BIO5 (5), and BIO1 (6) (Table 1). Additionally, in the study area, all future scenarios predicted a temperature increase, considering the temperature-related variables (BIO1, BIO5, and BIO10) and a decrease in precipitation in the related ones (BIO14, BIO17, and BIO18). The temperature variables had a reduction from LIG to LGM and an increase from LGM to present. On the other hand, the precipitation increased from the past to the present (Supplementary Figure S1), and those results were consistent in both GCMs.

TABLE 1

Table 1. Variable importance for each modeling algorithm, Generalized Linear Model (GLM), Multivariate Adaptive Regression Splines (MARS), MaxEnt, and Artificial Neural Networks (ANN), and for the ensemble model (average values).

Figure 2 displays how each variable’s predicted probability of occurrence is affected by the ensemble model, considering all the other constants. A medium to high suitability (>0.6) for the species is observed when (i) the annual mean temperature (BIO1) ranged between 15 and 17°C, (ii) the maximum temperature of the warmest month (BIO5) ranged from 31 to 36°C, (iii) the mean temperature of warmest quarter (BIO10) around 23–24°C, (iv) the precipitation of driest month (BIO14) decreased bellow 5 mm, (v) the precipitation of driest quarter (BIO17) ranged from 29 to 300 mm, and (vi) the precipitation of warmest quarter (BIO18) ranged from 19 to 32 mm. The precipitation-related variables response curves suggested that the species tolerate relatively low precipitation content, and, in general, the probability of occurrence decreased when those variables increased.

FIGURE 2

Figure 2. The ensemble model probability of occurrence response curves for the selected environmental variables. BIO1 (Annual Mean Temperature, °C), BIO5 (Maximum Temperature of Warmest Month, °C), BIO10 (Mean Temperature of Warmest Quarter, °C), BIO14 (Precipitation of Driest Month, mm), BIO17 (Precipitation of Driest Quarter, mm), and BIO18 (Precipitation of Warmest Quarter, mm).

The GLM displayed the worst performance (AUC = 0.871 and TSS = 0.648), and the ANN the best one (AUC = 0.951 and TSS = 0.789), based on the evaluation statistics (Table 2). The ensemble model showed higher predictive accuracy since the AUC = 0.944 was close to the theoretical maximum AUC value (Phillips et al., 2006), and the TSS = 0.755, indicated that the model had a better performance than 0.6, the value that defined a good model (Jones et al., 2010) than the three algorithms (GLM, MaxEnt, and MARS), but not when compared with the Artificial Neural Networks.

TABLE 2

Table 2. Evaluation statistics for each algorithm and the ensemble model.

3.2. Current potential distribution range

In the ensemble model (Supplementary Figure S3E, green and mustard), the potential C. ladanifer distribution with medium to high suitability (suitability >0.6) included the south and south-western part of the Iberian Peninsula, the center and north-east of Portugal, and spotted areas in Morocco and Algeria. The center and southern part of the IP, south-east France, a few sites in Italy (including Sardinia), Tunisia, and the Atlantic coast of Morocco (Supplementary Figure S3E, orange and yellow) were characterized as low to generally suitable regions (0.2 > suitability >0.6). These patterns were also observed in the GLM, ANN, and MARS models (Supplementary Figures S3A,C,D, respectively), with some differences, particularly in northern Africa. These three models showed a high agreement in their predictions, though projecting more high suitability areas (Supplementary Figures S3A,C,D, green) for the species compared to the ensemble model. The MaxEnt (Supplementary Figure S3B) was the most disruptive algorithm with less suitable areas predicted, namely the medium-high suitability.

Figure 3 displays the past (LIG, LGM, and MH), the current, and the future potential C. ladanifer distribution, under the global climate change scenarios RCP 4.5 and RCP 8.5, for the years the 2050s and 2070s. Regarding the LIG (Figure 3A), a considerable part of the Mediterranean North Africa area and the southern and western IP were found suitable. In the LGM (Figure 3B), an enormous retraction in the species distribution occurred, which remained spotted nearby the sea and absent inland. The model predicted refugia, in this period, on the eastern coast of Spain, the Sardinian west coast, and along the Algerian and Tunisian Mediterranean regions. The MH environmental conditions were more suitable for the species development than during the LGM, with a consequent C. ladanifer range expansion from that period until the current conditions (Figures 3C,D).

FIGURE 3

Figure 3. Predicted habitat suitability for the different scenarios. (A) LIG, (B) LGM, (C) MH, (D) Current, (E) RCP 4.5 2050, (F) RCP 4.5 2070, (G) RCP 8.5 2050, and (H) RCP8.5 2070. Suitability class range: 0.0–0.2, non-suitable area; 0.2–0.4, low-suitability area; 0.4–0.6, general-suitability area; 0.6–0.8, medium-suitability area and 0.8–1.0, high-suitability area. The maps were built with ArcGIS® Desktop v10.8.1.

The future projections (Figures 3E–H) unveil a considerable habitat loss across the North of Africa and the central part of the IP, including the regions with the highest suitability values, where a higher number of presence points were recorded (Figure 1). The RCP 4.5 in the years 2050 and 2070 displayed similar scenarios (Figures 3E,F), with the species’ putative existence in the western part of the IP. Contrarily, in the RCP 8.5 scenario (Figures 3G,H), a large area of central and northeastern Portugal and western Spain will be unsuitable or low-suitable, particularly for the RCP 8.5–2070 scenario.

From the LIG to the current species’ presence, the species’ disappearance in the LGM was well discernible. Those losses in suitability across time can be seen in the binary presence-absence maps (Supplementary Figure S4), where the past, present, and future distributions of C. ladanifer are displayed. In the present, the species had the broadest distribution considering all the periods and scenarios (Supplementary Figures S4A–D).

3.3. Rock rose distribution changes

The current binary presence-absence map (Supplementary Figure S4D) was used as a reference to identify potential shifts and trends in suitable habitat areas between the current and the other periods and scenarios. Suitability changes in past and future scenarios relative to current were classified into: (i) unsuitable habitat, (ii) habitat loss (contraction), (iii) habitat maintained, and (iv) habitat gained (expansion) (Figure 4 and Supplementary Figure S5D). In general, two patterns were observed: (1) habitat gained from past scenarios to the present and (2) habitat loss from the present to the future scenarios (Figure 5).

FIGURE 4

Figure 4. Changes in the C. ladanifer distribution area between scenarios. (A) LIG-Current, (B) LGM-Current, (C) MH-Current, (D) Current-RCP 4.5 2050, (E) Current-RCP 4.5 2070, (F) Current-RCP 8.5 2050, (G) Current-RCP 8.5 2070. The maps were built with ArcGIS® Desktop v10.8.1.

FIGURE 5

Figure 5. Changes in the C. ladanifer distribution area between scenarios.

From the LIG period to the present (Figure 4A), habitat gain is superior in the area to habitat loss. As expected, the changes in suitable habitat across time from the LGM to the present (Figure 4B) were characterized by a significant increase in the area due to improved climatic conditions (Supplementary Figure S5D).

For the future scenarios (Figures 4D–G, 5), the results showed that habitat loss will increase according to the climate scenario severity, and, consequently, the species’ presence will decrease. Under RCP 4.5 (Figures 4D,E), the habitat loss is similar for the years 2050 and 2070; in contrast, the habitat gain is higher in 2070 due to the emergence of newly suitable areas. The RCP 8.5 (Figures 4F,G), the worst climate scenario, will provoke a sounding suitable area reduction, particularly in the 2070s.

Habitat loss will be more evident in North Africa and, in central and southern IP, persisting small isolated populations across North Africa and South Iberia (increase in fragmentation) and in a large area of Central Portugal. The species will remain in southern France in the RCP 4.5 scenario but will disappear there in the RCP 8.5 scenario, suggesting that this area will become unsuitable in the future (Figure 4). The model indicated the occurrence of some habitat shifts in the past, and it predicted that the shifts would have a shallow extent in the future.

3.4. Palaeobotanical and genetic data

The palynological record (Figure 1 and Supplementary Table S3) showed very little evidence from the Late Pleistocene to the Holocene, all of which are in the IP (24). The two oldest records, dating from the Pleistocene (ca.142 and 90 ka) in Bajondillo and Padul-Granada (numbers 1 and 2), confirmed the presence of C. ladanifer in the south of Spain before the LGM. The other two Pleistocene records dated ca. 13 ka, after the LGM, originated from sites nearby the sea by the southern Portuguese coast (number 3) and in Galicia, North of the IP (PRD-4, number 4). These presences are congruent with refugia on the western coast of Portugal, as indicated in the model produced for the LGM by the current study (Figure 3B). The pollen records from the Holocene were distributed along a strip across the middle of the IP, from the central Portuguese coast to the northeastern Spanish coast (Figure 1). This distribution was under the species expansion area during the MH and onwards. Interestingly, the most recent C. ladanifer fossil pieces of evidence (late Holocene, around 1.3 ka) were found in the northeastern part of Spain, indicating the species expansion, both in this area and on the southern French coast, were recent events in this shrub history. Nevertheless, inferences from some palynological pieces of evidence should be made with caution: the ‘C. ladanifer-Type.’ The North of Africa was under-represented in the palynological records, probably a less studied region than Europe.

Moreover, C. ladanifer pollen is mainly dispersed by gravity at short distances, leading to a reduced presence in the pollen records. The North of Africa’s lack of palynological data might be partially overcome by the information given by the genetic data. Chloroplast and nuclear DNA information indicate that rockrose populations from the IP belong to two main population clusters that originated from independent colonization events from the North of Africa (Figure 1). Therefore, according to the genetic data, the presence of C. ladanifer in Africa should preclude the oldest indications of rockrose presence in the IP (Guzmán and Vargas, 2009; Quintela-Sabarís et al., 2011). The present distribution of the central populations’ genetic lineages indicated separated putative glacial refugia in the southeastern and western IP. This information is congruent with the past changes in habitat suitability from LIG through LGM to the MH, inferred from the models produced in the current study (Figures 3A–C).

4. Discussion

4.1. Current and future scenarios

In the discussion, we considered the following suitability area classes: 0.0–0.2, non-suitable; 0.2–0.4, low-suitable; 0.4–0.6, generally suitable; 0.6–0.8, medium-suitable; and 0.8–1.0, high-suitable.

The rock rose suitable habitat was predicted to occur in a wide area of the western Mediterranean region, particularly in the Iberian Peninsula, southern France, and northern Morocco and Algeria. Those findings confirmed the species’ ability to adapt to various climatic conditions, though the Mediterranean climate is prevalent. Indeed, the species distribution foreseen by the ensemble model is a good representation of the C. ladanifer current distribution (Figure 1, light green dots). The species’ medium-high probability of occurrence (>0.6) happens mostly in its native range, in the center and southwest of the IP, northern Morocco, and Algeria. These medium-high probability of occurrence were more frequent under the current conditions compared to the other scenarios, thus confirming the current climatic conditions as the most suitable for the species occurrence. Additionally, there is a good match between two divides, a bio-climatological and a lithological one, and the potential distribution of the species. A significant lithological divide exists between western acidic and eastern basic substrates [see Figure 1 in Buira et al. (2021)]. The other divide is due to a temperate climate north and northwest of the IP (the Pyrenees and the Cantabrian Mountains). At the same time, the rest of the region is characterized by the Mediterranean climate, with a noticeable summer drought (Buira et al., 2021, and references therein). In the binary ensemble map (Supplementary Figure S4D), these two divides frame the C. ladanifer potential distribution, south and westward oriented, avoiding basic soils and temperate climate.

The results showed that the species distribution range would be affected by climate change, regardless of the severity of the scenario. Indeed, heat waves are expected shortly for the Mediterranean region, but also a decrease in precipitation, in particular in this region, compared to global change worldwide (Vescovi et al., 2010; Grímsson et al., 2020), and future climate conditions will become less suitable for the species occurrence. Nevertheless, it seems to cope further with temperature excess than with the precipitation decrease, since the species’ distribution was mainly determined by the precipitation-related variables. Future temperature increase impact on rock rose distribution will probably be moderated by the species’ temperature resistance ecophysiology. Indeed, the production of labdanum, a sticky resin, increasingly during summer, photoprotects the leaves from UV light and helps control leave temperature and evapotranspiration decrease (Chaves et al., 1993; Nunez-Olivera et al., 1996).

According to the results, under the moderate emission scenario, the potentially suitable area will decrease both in 2050 and 2070 to a similar degree, yet, under the higher emission scenario, the suitable area will have a significant reduction in 2070 compared to the present. This reduction will be principally due to medium-high suitability areas decrease, from the present to the future, especially for the worst-case scenario, with substantial habitat loss in the center of the IP and northern Africa, a retraction westward oriented, with a potential fragmentation increase. This species’ habitat suitability will change considerably, even enduring drought and temperature increase, and the geographic distribution patterns will be affected by future climate change. Other authors uncovered that global change would negatively affect the Mediterranean species (e.g., Quercus sp.) distribution areas (Nicol-Pichard and Dubar, 1998; Çoban et al., 2020). Likewise, a decrease in suitability area is predicted for many other Mediterranean species’ geographical distribution, leading to shifts in geographical ranges and habitat loss, like with strawberry tree and others (e.g., Peñuelas et al., 2017; Almeida et al., 2022). Nevertheless, the rock rose will not migrate northerly; the suitable areas will contract toward the Portuguese coast, particularly in the worst-case scenarios. The northward migration of the species was a consistent result obtained by several authors, including Arbutus unedo (Santos et al., 2002; Abel-Schaad and López-Sáez, 2013; Rocha et al., 2014; Ribeiro et al., 2019; Almeida et al., 2022).

The species’ response to climate change in the past could further help in understanding its response in the future, which will be further explored in the coming point.

4.2. Past scenarios

The projections into the past revealed a range contraction from the LIG to the LGM and a continuous range expansion from the LGM through MH to the present (see Figure 5). Thus, during the LIG, the model projected a wide suitable area for C. ladanifer in the North of Africa and two main suitable areas in the South of the IP separated by the Guadalquivir basin: one to the southeast, encircled by Betic ranges and other to the west, spreading from the Algarve in the SW to the Tejo valley (Figure 3A) The pollen record partially supports these projections, with confirmed presence of C. ladanifer in the Betic area around 100 ka (Pons and Reille, 1988; Cortés-Sánchez et al., 2008). Moreover, the genetic data indicate that the African populations of C. ladanifer were the origin of the Iberian ones, which further supports the model projections about habitat suitability in the North of Africa (Quintela-Sabarís et al., 2011, 2012).

The results indicated that the climatic conditions in the LGM were less suitable for the rockrose, with the disappearance of the areas with a medium-high probability of occurrence (dark green in Figure 3B) and a general contraction of the areas with a general probability of occurrence (orange in Figure 3B). Thus, three main areas emerge as potential refugia for C. ladanifer during LGM. First, a circle comprising East Spain (Valencia and Murcia regions) and the NE coast of Algeria and Tunisia. This circle is a typical hotspot of biodiversity and common putative refugia during this period (Médail and Diadema, 2009, and references therein; Feliner, 2011). Second, the central coastal region of Portugal has also been found to be a refugium for the Mediterranean species maritime pine and strawberry tree, based on the fossil record (Ribeiro et al., 2001; Almeida et al., 2022). Third, a cryptic refugium in the Betic and Rif coastal regions (South of the IP, North of Morocco, respectively). The fossil evidence for this period supports some of the model results, with the registered presence of C. ladanifer-Type pollen 13 ka in palynological sites in the South of Portugal and on the Galician coast (northwestern IP).

The IP is considered one of Europe’s most important Pleistocene glacial refugia. The ‘refugia within-refugia hypothesis’ (Gómez and Lunt, 2007) claims that this region is unlikely to be a unique continuous refugium area due to its high environmental variability. Instead, it is divided into multiple and isolated refugia of different sizes. The ‘refugia within refugia’ have been demonstrated to be a common pattern, which explains the current distribution of several species of animals and plants in intraspecific lineages in the IP (Gómez and Lunt, 2007, and references therein). Indeed, the LGM refugia identification should help understand the species’ postglacial migration during the Holocene, from those regions and its adaptation to local conditions (Gavin et al., 2014).

The scenario depicted by our model considering C. ladanifer distribution during the LGM also fitted this pattern and, despite the limited amount of pollen evidence from that period, supported by the current species’ genetic studies. Thus, four intraspecific lineages in C. ladanifer (whose distribution coincides with the refugia inferred by the model approach: Figure 3B) have been identified based on chloroplast and nuclear markers (Quintela-Sabarís et al., 2011, 2012): four lineages, (1) Rif region, (2) Betic area, (3) central coastal Portugal, and (4) East of the IP.

Since the LGM, the C. ladanifer suitable area has progressively increased, both in the North of Morocco and in the IP, and the species’ expansion reached its current distribution range. This model projection is also supported by the increasing number of pollen occurrences in the center of the IP and a more recent expansion to the northeastern IP region (around 1.3 ka, López-Sáez et al., 2009) and South of France. Viruel et al. (2020), in their study, found support for the carob tree persistence in Moroccan and Iberian (Betic region) refugia during the LGM, confirmed both with molecular and potential historical range dynamics by SDMs reconstruction. However, this species had multiple domestications events from native populations throughout the Mediterranean basin. Besides rock rose fossil pieces of evidence, the existence of a putative central coastal refuge in Portugal (4) was also supported by similar findings in other thermophilous Mediterranean taxa in this sheltered area due to the oceanic influence, such as olive tree, maritime pine, strawberry tree, and mastic gum (Ribeiro et al., 2001; Figueiral and Terral, 2002; Almeida et al., 2022).

Species’ postglacial expansion is dominated by the refugia placed in the ‘leading edge,’ whereas the populations’ expansion behind that front diffuses at a slower rate (Hewitt, 2001). This seems to be also the case in C. ladanifer, as lineages with refugia in Portugal and, especially East of the IP, expanded to the north and currently occupy most of the distribution area of this species. In contrast, the expansion from the Betic refugium was limited to the North by the presence of a physical barrier (the Betic mountain ranges) and to the West by the contact with the area occupied by other lineages. Indeed, Guzmán and Vargas (2009), in their phylogeography study using organellar sequences, claimed that the C. ladanifer origin and the subspecies divergence happened in the Upper Pleistocene, before the Strait of Gibraltar onset and the Mediterranean climate establishment. Their study indicated that the Cistus species arose after the Mediterranean Sea re-emergence (5.59–5.33 Ma) and the Mediterranean climate launch (3.2 Ma). The sea possibly sheltered the species during ice ages, and the current species distribution is the result of migration by animals eating the capsules and from the different detected putative refugia, seen both in the molecular and the past projected modeled suitability of occurrence (LGM).

Indeed, this species moved westwards from the eastern refugia (3), northerly from North Africa (1), and in a fanned way from the central Portugal refugia (4), during the Holocene warming up, in a relatively quick spread pace since the MH model predicted that the habitat suitability matched, in a certain extent, with the current species occurrence (Figures 3B–D). Those findings are in agreement with the LGM-modeled map information (Figure 3B), consistent with the hypothesis of the two refugia unveiled in eastern and southern IP, by the sea, another in North Africa, and the fourth one in central Portugal, also nearby the sea. The oceanic influence could have played a protective role against the last glaciation’s most arduous climate and was ultimately responsible for the species’ survival, as found in a study with the strawberry tree (Ribeiro et al., 2017).

As mentioned previously, the C. ladanifer seeds’ primary dispersal strategy is barochory, and this species’ gene flow is mainly local and that the seed dispersal is a significant component of gene flow, with the pollination dropping to zero at ca. 30 m (Quintela-Sabarís et al., 2012, and references therein). However, a secondary dispersal strategy made by animals (ants and deers) is an essential source of the species’ long-distance dispersal, with ready-to-germinate seeds spread by deer (Malo and Suarez, 1996; Bastida and Talavera, 2002). Wind and water are not to be discarded as additional seed spread vectors since the seed is very light and small. Therefore, this species’ high structure and genetic differentiation was consistent with low gene flow, compared to species with generally high gene flow — as in wind-pollinated species — such as pines (e.g., Ribeiro et al., 2001), which promotes genetic drift and the formation of clear lineages/genetic clusters (Quintela-Sabarís et al., 2011, 2012).

In summary, the distribution projections to the past offered by the models are congruent with the other two independent approaches (fossil records and genetics). This confirms the validity of our modeling approach to infer the effect of past changes in climate on C. ladanifer distribution and supported projections into the future under global change scenarios. This species will be dramatically affected by global warming, despite its resilience to drought and temperature. As commonly observed in other Mediterranean species, North migration is not the picture. Instead, the species will collapse westwards toward the Portuguese coast. The area of this species was predicted to be much reduced in the future, and the information given by this paper might help to understand the consequences for other equally drought-resistant species in this region. Policymakers should be put under pressure for urgent measures to mitigate the warming-up tendency and biodiversity depletion in this region, as well as expected for other thermophilous Mediterranean species. Certainly, habitat loss and spatial distribution conversion are predictable effects of climate change for Mediterranean species, and species acclimation, adaptation, and altitudinal and latitudinal migrations are expected (Carrión and Van Geel, 1999; Barredo et al., 2016). The Mediterranean region has many rare and locally endemic taxa that survived as small populations, many of which are threatened by habitat transformation and climate change (Peñuelas et al., 2017).

Additionally, the use of big data was also crucial, in this case for the species’ range-wide modeling, due to the immense amount of information requested: presence points, layers of bioclimatic variables, different periods in the past and future, different algorithms, and GCMs. Moreover, the current study is an example of how big data allows us to predict the future and reconstruct the past.

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

MR and PF conceived the study idea. MRF and AMA handled the species occurrence data preparation and bias correction, conducted the analyses using the GIS and the R package biomod2, and wrote the first draft of the manuscript. All authors significantly contributed to ideas and revision and editing and revising the manuscript.

Funding

This research was co-financed by the Regional Operational Program Centro 2020, Portugal 2020 and European Union, through European Fund for Regional Development (ERDF) under the program CULTIVAR CENTRO-01-0145-FEDER-000020 (MR, PF, MF, and AA) (https://icultivar.pt/), and further funded by FCT – Foundation for Science and Technology under the projects UIDB/00239/2020 (MR and NR), UIDB/00681/2020 (MR), and UIDB/ 05183/2020 (PF) (https://www.fct.pt/).

Acknowledgments

Palaeobotanic data were obtained from the Neotoma Paleoecology Database (http://www.neotomadb.org), and the work of data contributors, data stewards, and the Neotoma community is gratefully acknowledged. The proofreading made by Isabele Lavado is also acknowledged.

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.

Supplementary material

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

References

Abel-Schaad, D., and López-Sáez, J. A. (2013). Vegetation changes in relation to fire history and human activities at the Peña Negra mire (Bejar range, Iberian Central Mountain system, Spain) during the past 4,000 years. Veg. Hist. Archaeobot. 22, 199–214. doi: 10.1007/s00334-012-0368-9

CrossRef Full Text | Google Scholar

Alcamo, J., Moreno, J. M., Nováky, B., Bindi, M., Corobov, R., Devoy, R. J. N., et al. (Eds.). (2007). Europe. Climate Change. Cambridge: Cambridge University Press.

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

CrossRef Full Text | Google Scholar

Almeida, A. M., Martins, M. J., Campagnolo, M. L., Fernandez, P., Albuquerque, T., Gerassis, S., et al. (2022). Prediction scenarios of past, present, and future environmental suitability for the Mediterranean species Arbutus unedo L. Sci. Rep. 12:84. doi: 10.1038/s41598-021-03996-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Alonso Ponce, R., Águeda, B., Ágreda, T., Modrego, M. P., Aldea, J., Fernández-Toirán, L. M., et al. (2011). Rockroses and Boletus edulis ectomycorrhizal association: realized niche and climatic suitability in Spain. Fungal Ecol. 4, 224–232. doi: 10.1016/j.funeco.2010.10.002

CrossRef Full Text | Google Scholar

Barredo, J. I., Caudullo, G., and Dosio, A. (2016). Mediterranean habitat loss under future climate conditions: assessing impacts on the Natura 2000 protected area network. Appl. Geogr. 75, 83–92. doi: 10.1016/j.apgeog.2016.08.003

CrossRef Full Text | Google Scholar

Bastida, F., and Talavera, S. (2002). Temporal and spatial patterns of seed dispersal in two Cistus species (Cistaceae). Ann. Bot. 89, 427–434. doi: 10.1093/aob/mcf065

PubMed Abstract | CrossRef Full Text | Google Scholar

Buira, A., Fernández-Mazuecos, M., Aedo, C., and Molina-Venegas, R. (2021). The contribution of the edaphic factor as a driver of recent plant diversification in a Mediterranean biodiversity hotspot. J. Ecol. 109, 987–999. doi: 10.1111/1365-2745.13527

CrossRef Full Text | Google Scholar

Carrión, J. S., and Van Geel, B. (1999). Fine-resolution upper Weichselian and Holocene palynological record from Navarrés (Valencia, Spain) and a discussion about factors of Mediterranean forest succession. Rev. Palaeobot. Palyno. 106, 209–236. doi: 10.1016/S0034-6667(99)00009-3

CrossRef Full Text | Google Scholar

Chaves, N., Escudero, J. C., and Gutiérrez-Merino, C. (1993). Seasonal variation of exudate of Cistus ladanifer. J. Chem. Ecol. 19, 2577–2591. doi: 10.1007/BF00980692

PubMed Abstract | CrossRef Full Text | Google Scholar

Çoban, H. O., Örücü, Ö. K., and Arslan, E. S. (2020). Max Ent modeling for predicting the current and future potential geographical distribution of Quercus libani Olivier. Sustainability 12:2671. doi: 10.3390/su12072671

CrossRef Full Text | Google Scholar

Cortés-Sánchez, M., Morales-Muñiz, A., Simón-Vallejo, M. D., Bergadà-Zapata, M. M., Delgado-Huertas, A., López-García, P., et al. (2008). Palaeoenvironmental and cultural dynamics of the coast of Málaga (Andalusia, Spain) during the upper Pleistocene and early Holocene. Quat. Sci. Rev. 27, 2176–2193. doi: 10.1016/j.quascirev.2008.03.010

CrossRef Full Text | Google Scholar

Cowling, R. M., Rundel, P. W., Lamont, B. B., Kalin Arroyo, M., and Arianoutsou, M. (1996). Plant diversity in mediterranean-climate regions. TREE 11, 362–366. doi: 10.1016/0169-5347(96)10044-6

CrossRef Full Text | Google Scholar

Delgado, J. A., Serrano, J. M., López, F., and Acosta, F. J. (2001). Heat shock, mass-dependent germination, and seed yield as related components of fitness in Cistus ladanifer. Environ. Exp. Bot. 46, 11–20. doi: 10.1016/S0098-8472(01)00076-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Demoly, J. P., and Montserrat, P. (1993). “Cistus L” in Flora Iberica. eds. S. Castroviejo, C. Aedo, S. Cirujano, M. Laínz, P. Montserrat, and R. Morales, et al. (Madrid, Spain: Real Jardín Botánico, CSIC), 319–337.

Google Scholar

Dormann, C. F., Elith, J., Bacher, S., Buchmann, C., Carl, G., Carré, G., et al. (2013). Collinearity: a review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46. doi: 10.1111/j.1600-0587.2012.07348.x

CrossRef Full Text | Google Scholar

Feliner, G. N. (2011). Southern European glacial refugia: a tale of tales. Taxon 60, 365–372. doi: 10.1002/tax.602007

CrossRef Full Text | Google Scholar

Fernández-Mazuecos, M., and Vargas, P. (2010). Ecological rather than geographical isolation dominates quaternary formation of Mediterranean Cistus species. Mol. Ecol. 19, 1381–1395. doi: 10.1111/j.1365-294X.2010.04549.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira, M. R., Almeida, A. M., Quintela-Sabarís, C., Roque, N., Fernandez, P., and Ribeiro, M. M. (2021). The role of littoral cliffs in the niche delimitation on a microendemic plant facing climate change. PLoS One 16:e0258976. doi: 10.1371/journal.pone.0258976

PubMed Abstract | CrossRef Full Text | Google Scholar

Fielding, A. H., and Bell, J. F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49. doi: 10.1017/S0376892997000088

CrossRef Full Text | Google Scholar

Figueiral, I., and Terral, J. F. (2002). Late quaternary refugia of Mediterranean taxa in the Portuguese Estremadura: charcoal based palaeovegetation and climatic reconstruction. Quat. Sci. Rev. 21, 549–558. doi: 10.1016/S0277-3791(01)00022-1

CrossRef Full Text | Google Scholar

Frazão, D. F., Raimundo, J. R., Domingues, J. L., Quintela-Sabarís, C., Gonçalves, J. C., and Delgado, F. (2018). Cistus ladanifer (Cistaceae): a natural resource in Mediterranean-type ecosystems. Planta 247, 289–300. doi: 10.1007/s00425-017-2825-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, J. H. (1991). Multivariate adaptive regression splines. Ann. Stat. 19, 1–67. doi: 10.1214/aos/1176347963

CrossRef Full Text | Google Scholar

Gallego, J. C. A., Caro, J. G., Campos, V. H., and Lobón, N. C. (2020). Effect of leaf litter from Cistus ladanifer L. on the germination and growth of accompanying shrubland species. Plan. Theory 9:593. doi: 10.3390/plants9050593

PubMed Abstract | CrossRef Full Text | Google Scholar

Gavin, D. G., Fitzpatrick, M. C., Gugger, P. F., Heath, K. D., Rodríguez-Sánchez, F., Dobrowski, S. Z., et al. (2014). Climate refugia: joint inference from fossil records, species distribution models and phylogeography. New Phytol. 204, 37–54. doi: 10.1111/nph.12929

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez, A., and Lunt, D. H. (2007). “Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula” in Phylogeography of Southern European Refugia: Evolutionary Perspectives on the Origins and Conservation of European Biodiversity. eds. S. Weiss and N. Ferrand (Dordrecht: Springer Netherlands), 155–188.

Google Scholar

Gomulkiewicz, R., and Houle, D. (2009). Demographic and genetic constraints on evolution. Am. Nat. 174, E218–E229. doi: 10.1086/645086

CrossRef Full Text | Google Scholar

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., and Moore, R. (2017). Google earth engine: planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 202, 18–27. doi: 10.1016/j.rse.2017.06.031

CrossRef Full Text | Google Scholar

Grímsson, F., Bouchal, J. M., Xafis, A., and Zetter, R. (2020). Combined LM and SEM study of the middle Miocene (Sarmatian) palynoflora from the Lavanttal Basin, Austria: part V. Magnoliophyta 3 – Myrtales to Ericales. Grana 59, 127–193. doi: 10.1080/00173134.2019.1696400

PubMed Abstract | CrossRef Full Text | Google Scholar

Guiot, J., and Cramer, W. (2016). Climate change: the 2015 Paris agreement thresholds and Mediterranean basin ecosystems. Science 354, 465–468. doi: 10.1126/science.aah5015

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzmán, B., Lledó, M. D., and Vargas, P. (2009). Adaptive radiation in Mediterranean Cistus (Cistaceae). PLoS One 4:e6362. doi: 10.1371/journal.pone.0006362

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzmán, B., and Vargas, P. (2005). Systematics, character evolution, and biogeography of Cistus L. (Cistaceae) based on ITS, trnL-trnF, and matK sequences. Mol. Phylogen. Evol. 37, 644–660. doi: 10.1016/j.ympev.2005.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Guzmán, B., and Vargas, P. (2009). Long-distance colonization of the Western Mediterranean by Cistus ladanifer (Cistaceae) despite the absence of special dispersal mechanisms. J. Biogeogr. 36, 954–968. doi: 10.1111/j.1365-2699.2008.02040.x

CrossRef Full Text | Google Scholar

Hao, T., Elith, J., Guillera-Arroita, G., and Lahoz-Monfort, J. J. (2019). A review of evidence about use and performance of species distribution modelling ensembles like BIOMOD. Divers. Distrib. 25, 839–852. doi: 10.1111/ddi.12892

CrossRef Full Text | Google Scholar

Hewitt, G. M. (2001). Speciation, hybrid zones and phylogeography — or seeing genes in space and time. Mol. Ecol. 10, 537–549. doi: 10.1046/j.1365-294x.2001.01202.x

CrossRef Full Text | Google Scholar

Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G., and Jarvis, A. (2005). Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. doi: 10.1002/joc.1276

CrossRef Full Text | Google Scholar

Hopfield, J. J. (1982). Neural networks and physical systems with emergent collective computational abilities. Proc. Natl. Acad. Sci. U. S. A. 79, 2554–2558. doi: 10.1073/pnas.79.8.2554

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones, C. C., Acker, S. A., and Halpern, C. B. (2010). Combining local- and large-scale models to predict the distributions of invasive plant species. Ecol. Appl. 20, 311–326. doi: 10.1890/08-2261.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jump, A. S., Hunt, J. M., Martinez-Izquierdo, J. A., and Penuelas, J. (2006). Natural selection and climate change: temperature-linked spatial and temporal trends in gene frequency in Fagus sylvatica. Mol. Ecol. 15, 3469–3480. doi: 10.1111/j.1365-294X.2006.03027.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., et al. (2017). Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4:170122. doi: 10.1038/sdata.2017.122

PubMed Abstract | CrossRef Full Text | Google Scholar

Leroy, B., Meynard, C. N., Bellard, C., and Courchamp, F. (2016). Virtualspecies, an R package to generate virtual species distributions. Ecography 39, 599–607. doi: 10.1111/ecog.01388

CrossRef Full Text | Google Scholar

Loarie, S. R., Duffy, P. B., Hamilton, H., Asner, G. P., Field, C. B., and Ackerly, D. D. (2009). The velocity of climate change. Nature 462, 1052–1055. doi: 10.1038/nature08649

CrossRef Full Text | Google Scholar

López-Sáez, J. A., López-Merino, L., Mateo, M., Serrano, O., Pérez-Díaz, S., and Serrano, L. (2009). Palaeoecological potential of the marine organic deposits of Posidonia oceanica: a case study in the NE Iberian Peninsula. Palaeogeogr. Palaeoclimatol. Palaeoecol. 271, 215–224. doi: 10.1016/j.palaeo.2008.10.020

CrossRef Full Text | Google Scholar

Malo, J. E., and Suarez, F. (1996). Cistus ladanifer recruitment – not only fire, but also deer. Acta Oecol. 17, 55–60.

Google Scholar

Malo, J. E., and Suárez, F. (1998). The dispersal of a dry-fruited shrub by red deer in a Mediterranean ecosystem. Ecography 21, 204–211. doi: 10.1111/j.1600-0587.1998.tb00673.x

CrossRef Full Text | Google Scholar

McCullagh, P., and Nelder, J. A. (1983). Generalized Linear Models. New York: Routledge

Google Scholar

Médail, F., and Diadema, K. (2009). Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J. Biogeogr. 36, 1333–1345. doi: 10.1111/j.1365-2699.2008.02051.x

CrossRef Full Text | Google Scholar

Mendes, P., Meireles, C., Vila-Viçosa, C., Musarella, C., and Pinto-Gomes, C. (2015). Best management practices to face degraded territories occupied by Cistus ladanifer shrublands – Portugal case study. Plant Biosyst. 149, 494–502. doi: 10.1080/11263504.2015.1040483

CrossRef Full Text | Google Scholar

Myers, N., Mittermeier, R. A., Mittermeier, C. G., da Fonseca, G. A. B., and Kent, J. (2000). Biodiversity hotspots for conservation priorities. Nature 403, 853–858. doi: 10.2307/239995710.1038/35002501

CrossRef Full Text | Google Scholar

Nicol-Pichard, S., and Dubar, M. (1998). Reconstruction of late-glacial and holocene environments in Southeast France based on the study of a 66-m long core from biot, Alpes Maritimes. Veg. Hist. Archaeobot. 7, 11–15. doi: 10.1007/bf01393413

CrossRef Full Text | Google Scholar

Nunez-Olivera, E., Martinez-Abaigar, J., and Escudero, J. (1996). Adaptability of leaves of Cistus ladanifer to widely varying environmental conditions. Funct. Ecol. 10, 636–646. doi: 10.2307/2390174

CrossRef Full Text | Google Scholar

Ochoa-Ochoa, L. M., Flores-Villela, O. A., and Bezaury-Creel, J. E. (2016). Using one vs. many, sensitivity and uncertainty analyses of species distribution models with focus on conservation area networks. Ecol. Model. 320, 372–382. doi: 10.1016/j.ecolmodel.2015.10.031

CrossRef Full Text | Google Scholar

Otto-Bliesner, B. L., Rosenbloom, N., Stone, E. J., McKay, N. P., Lunt, D. J., Brady, E. C., et al. (2013). How warm was the last interglacial? New model–data comparisons. Phil. Trans. R. Soc. A 371:20130097. doi: 10.1098/rsta.2013.0097

PubMed Abstract | CrossRef Full Text | Google Scholar

Peñuelas, J., Sardans, J., Filella, I., Estiarte, M., Llusià, J., Ogaya, R., et al. (2017). Impacts of global change on Mediterranean forests and their services. Forests 8:463. doi: 10.3390/f8120463

CrossRef Full Text | Google Scholar

Pérez-García, F. (1997). Germination of Cistus ladanifer seeds in relation to parent material. Plant Ecol. 133, 57–62. doi: 10.1023/a:1009776910683

CrossRef Full Text | Google Scholar

Phillips, S. J., Anderson, R. P., and Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259. doi: 10.1016/j.ecolmodel.2005.03.026

CrossRef Full Text | Google Scholar

Pons, A., and Reille, M. (1988). The holocene- and upper pleistocene pollen record from Padul (Granada, Spain): a new study. Palaeogeogr. Palaeoclimatol. Palaeoecol. 66, 243–263. doi: 10.1016/0031-0182(88)90202-7

CrossRef Full Text | Google Scholar

Quintela-Sabarís, C., Ribeiro, M. M., Poncet, B., Costa, R., Castro-Fernández, D., and Fraga, M. I. (2012). AFLP analysis of the pseudometallophyte Cistus ladanifer: comparison with cpSSRs and exploratory genome scan to investigate loci associated to soil variables. Plant Soil 359, 397–413. doi: 10.1007/s11104-012-1221-8

CrossRef Full Text | Google Scholar

Quintela-Sabarís, C., Vendramin, G. G., Castro-Fernández, D., and Isabel Fraga, M. (2011). Chloroplast DNA phylogeography of the shrub Cistus ladanifer L. (Cistaceae) in the highly diverse Western Mediterranean region. Plant Biol. 13, 391–400. doi: 10.1111/j.1438-8677.2010.00371.x

PubMed Abstract | CrossRef Full Text | Google Scholar

R-Core-Team (2021). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna. Austria. Available at: http://www.R-project.org/

Google Scholar

Ribeiro, M. M., Piotti, A., Ricardo, A., Gaspar, D., Costa, R., Parducci, L., et al. (2017). Genetic diversity and divergence at the Arbutus unedo L. (Ericaceae) westernmost distribution limit. PLoS One 12:e0175239. doi: 10.1371/journal.pone.0175239

PubMed Abstract | CrossRef Full Text | Google Scholar

Ribeiro, M. M., Plomion, C., Petit, R., Vendramin, G. G., and Szmidt, A. E. (2001). Variation in chloroplast single-sequence repeats in Portuguese maritime pine (Pinus pinaster Ait.). Theor. Appl. Genet. 102, 97–103. doi: 10.1007/s001220051623

CrossRef Full Text | Google Scholar

Ribeiro, M. M., Roque, N., Ribeiro, S., Gavinhos, C., Castanheira, I., Quinta-Nova, L., et al. (2019). Bioclimatic modeling in the last glacial maximum, mid-Holocene and facing future climatic changes in the strawberry tree (Arbutus unedo L.). PLoS One 14:e0210062. doi: 10.1371/journal.pone.0210062

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocha, J., Almeida, R., Almeida da Silva, R., Amich, F., Martins, Á., Almeida, P., et al. (2014). Biogeographic trends of endemic and subendemic flora in the western Iberian Peninsula under scenarios of future climate change. Lazaroa 35, 19–35. doi: 10.5209/rev_LAZA.2014.v35.44196

CrossRef Full Text | Google Scholar

Rundel, P. W., Arroyo, M. T. K., Cowling, R. M., Keeley, J. E., Lamont, B. B., and Vargas, P. (2016). Mediterranean biomes: evolution of their vegetation, floras, and climate. Annu. Rev. Ecol. Evol. Syst. 47, 383–407. doi: 10.1146/annurev-ecolsys-121415-032330

CrossRef Full Text | Google Scholar

Santos, F. D., Forbes, K., and Moita, R. eds. (2002). Climate Change in Portugal. Scenarios, Impacts and Adaptation Measures - SIAM Project. Lisbon: Gradiva.

Google Scholar

Serra-Varela, M. J., Grivet, D., Vincenot, L., Broennimann, O., Gonzalo-Jiménez, J., and Zimmermann, N. E. (2015). Does phylogeographical structure relate to climatic niche divergence? A test using maritime pine (Pinus pinaster Ait.). Glob. Ecol. Biogeogr. 24, 1302–1313. doi: 10.1111/geb.12369

CrossRef Full Text | Google Scholar

Sosa, T., Valares, C., Alías, J. C., and Chaves Lobón, N. (2010). Persistence of flavonoids in Cistus ladanifer soils. Plant Soil 337, 51–63. doi: 10.1007/s11104-010-0504-1

CrossRef Full Text | Google Scholar

Stewart, S. B., Fedrigo, M., Kasel, S., Roxburgh, S. H., Choden, K., Tenzin, K., et al. (2022). Predicting plant species distributions using climate-based model ensembles with corresponding measures of congruence and uncertainty. Divers. Distrib. 28, 1105–1122. doi: 10.1111/ddi.13515

CrossRef Full Text | Google Scholar

Suc, J.-P., Popescu, S.-M., Fauquette, S., Bessedik, M., Jiménez-Moreno, G., Bachiri Taoufiq, N., et al. (2018). Reconstruction of Mediterranean flora, vegetation and climate for the last 23 million years based on an extensive pollen dataset. Ecol. Mediterr. 44, 53–85. doi: 10.3406/ecmed.2018.2044

CrossRef Full Text | Google Scholar

Svenning, J.-C., Fløjgaard, C., Marske, K. A., Nógues-Bravo, D., and Normand, S. (2011). Applications of species distribution modeling to paleobiology. Quat. Sci. Rev. 30, 2930–2947. doi: 10.1016/j.quascirev.2011.06.012

CrossRef Full Text | Google Scholar

Swets, J. A. (1988). Measuring the accuracy of diagnostic systems. Science 240, 1285–1293. doi: 10.1126/science.3287615

CrossRef Full Text | Google Scholar

Thompson, J. D. (2005). Plant Evolution in the Mediterranean. Oxford: Oxford Oxford University Press

Google Scholar

Thomson, A. M., Calvin, K. V., Smith, S. J., Kyle, G. P., Volke, A., Patel, P., et al. (2011). RCP4.5: a pathway for stabilization of radiative forcing by 2100. Clim. Chang. 109, 77–94. doi: 10.1007/s10584-011-0151-4

CrossRef Full Text | Google Scholar

Thuiller, W., Georges, D., Gueguen, M., Engler, R., and Breiner, F. (2021). biomod2: Ensemble platform for species distribution modeling. R package version 3.5.1. Available at: https://cran.r-project.org/package=biomod2

Google Scholar

Valente, L. M., and Vargas, P. (2013). Contrasting evolutionary hypotheses between two mediterranean-climate floristic hotspots: the cape of southern Africa and the Mediterranean Basin. J. Biogeogr. 40, 2032–2046. doi: 10.1111/jbi.12156

CrossRef Full Text | Google Scholar

van Vuuren, D. P., Edmonds, J., Kainuma, M., Riahi, K., Thomson, A., Hibbard, K., et al. (2011). The representative concentration pathways: an overview. Clim. Chang. 109, 5–31. doi: 10.1007/s10584-011-0148-z

CrossRef Full Text | Google Scholar

Vescovi, E., Ammann, B., Ravazzi, C., and Tinner, W. (2010). A New late-glacial and Holocene record of vegetation and fire history from Lago del Greppo, northern Apennines, Italy. Veg. Hist. Archaeobot. 19, 219–233. doi: 10.1007/s00334-010-0243-5

CrossRef Full Text | Google Scholar

Viruel, J., Le Galliot, N., Pironon, S., Nieto Feliner, G., Suc, J.-P., Lakhal-Mirleau, F., et al. (2020). A strong east–West Mediterranean divergence supports a new phylogeographic history of the carob tree (Ceratonia siliqua, Leguminosae) and multiple domestications from native populations. J. Biogeogr. 47, 460–471. doi: 10.1111/jbi.13726

CrossRef Full Text | Google Scholar

Williams, J. W., Grimm, E. C., Blois, J. L., Charles, D. F., Davis, E. B., Goring, S. J., et al. (2018). The Neotoma paleoecology database, a multiproxy, international, community-curated data resource. Quat. Res. 89, 156–177. doi: 10.1017/qua.2017.105

CrossRef Full Text | Google Scholar

Wróblewska, A., and Mirski, P. (2018). From past to future: impact of climate change on range shifts and genetic diversity patterns of circumboreal plants. Reg. Environ. Chang. 18, 409–424. doi: 10.1007/s10113-017-1208-3

CrossRef Full Text | Google Scholar