1. Introduction
During the century, biological invasions have gained attention as a paradigm of ecological disruption, historical variation, and ecological integration [1,2,3,4,5,6]. Undoubtedly, biological invasions and climate change are two key drivers behind current biodiversity patterns [7,8].
The biosphere has entered an era of rapid environmental changes resulting in significant negative impacts on ecosystem resilience [9,10,11,12]. However, scientific and societal unknowns make it challenging to predict how global environmental alterations, such as climate change and biological invasions, will affect species interactions and ecosystem processes in the future [4,13,14]. In the long term, these changes may have interacting effects and confound the consequences associated with each individual driver [14,15,16].
The expected global climate change is expected to further expand the risk of plant invasion through ecosystem disturbance and enhanced competition due to elevated CO2 atmospheric levels [2,17,18,19,20,21]. It is in combination with these pressures that global warming becomes potentially more problematic [22,23,24,25,26,27].
Expanded risk from invasive plants owing to climate changes has been identified for several species [28,29,30,31,32,33,34]. However, climate change might also decrease invasive plant competitiveness if conditions become climatically unsuitable [21]. In ecological research, a group of previous publications indicate that climate change is likely to favour particular invasive species (e.g., [2,20,35]), but few authors have identified specific concerns for climate change in relation with invasive species [14].
It is useful to scrutinize possible climate change responses of invasive species in order to (i) encourage discussion about the unique consequences of climate change for invasive species, (ii) identify key hypotheses that should be tested to improve the conceptual framework supporting our current understanding of both phenomena, and finally, (iii) advise strategic planning and the optimal allocation of the efforts devoted to control biological invasions and preserve biodiversity [14].
Many issues regarding the ecology of invasive alien species (IAS) in forest ecosystems cannot be answered directly by researchers using traditional hypothesis testing methods [36]. The early understanding of an invasion can gain relevance by using modelling tools that synthesize available information to outline potential management options [37]. Change scenarios are frequently used with species distribution models (SDMs) to assess shifts in species range induced by climate change [38,39,40]. Consequently, the relationships between species’ climatic and habitat preferences need to be assessed in order to recognize and predict changes in community composition [41,42,43,44]. If these relationships are assumed, then local stressors could be proactively identified to improve tolerance and strengthen resilience to global climate change [42,45]. Furthermore, innovations in predictive modelling techniques (e.g., [46]) and modelling software allow for applying and comparing several modelling methods, and selecting the best-performing SDM for ecology and conservation. SDMs correlate the observed distribution of a species to a set of environmental predictors, including climate data, and use this relationship to project the potential distribution into maps of future habitat suitability [47,48,49].
In their primordial state, oceanic islands typically supported a substantial proportion of endemic species, many of which have been lost as a direct consequence of anthropogenic effects, such as human habitation [50,51,52,53]. The biodiversity “crisis” is in need of urgent action since these remote islands are considered fragile ecosystems and highly vulnerable to biological invasions [50,51,54,55,56,57,58]. Thus, the biogeography and endemic biodiversity of oceanic islands makes them confined ecosystems that might serve as hot spots and laboratories for the study of ecological change [59]. It is notable they are a home to unique forms of endemic organisms that have evolved in isolation over millions of years [59,60,61]. Some of these ecosystems are the most sensitive and vulnerable to climate change and may therefore provide valuable monitoring sites for detecting the onset of human-induced global climate change [61].
According to Ferreira et al. [62], the Azorean archipelago provides a unique opportunity to study species–climate change relationships. Furthermore, the greatest expected impact of global climate change in the Azores might be a change in seasonal precipitation distribution, with wetter winters and drier conditions in the remaining seasons [62,63]. This change could have a potential effect on the islands’ water resources [62,63] and consequently on the vegetation.
Particularly, forest ecosystems are under pressure, as climate change may threaten forest-dependent species across a wide range of species groups [64,65]. The vulnerability of forest biodiversity, alongside with the dependency of many terrestrial species on forest ecosystems, highlights the leading role of forests in the global ecological sustainability [65,66,67]. To promote the sustainability of forests, one of the key issues for researchers is to recognize how site quality, species distribution, and eventually plant growth are affected by site conditions and how environmental variables may respond with changes in climate [68].
Pittosporum undulatum Vent. is the most widespread invasive woody plant in the Azores, where it is considered as one of the priority species for control actions. Another potential threat to native Azorean biota is Acacia melanoxylon R. Br. In contrast, Morella faya (Aiton) Wilbur is a native plant facing several threats, including plant invasions [69,70,71]. Current distribution of P. undulatum and A. melanoxylon are related to climate, topography, and some land-use effects. Furthermore, analysis of the areas under risk of invasion showed that protected areas are under potential threat [71,72].
The purpose of this study is to identify possible limitations when making SDM projections based on predicted climate change. A modelling application, BIOMOD 2 [73], which allows comparing results obtained using different modelling techniques, was used. We aimed to predict how habitat suitability might change in the future as a result of climate change, and to identify possible limitations of using solely climatic variables in modelling, as opposed to the inclusion of other data sources, namely topographic descriptors. In addition, we compared generalized linear models (GLMs), generalized additive models (GAMs), and a selection of machine learning techniques to model the potential change in the distribution range of P. undulatum, A. melanoxylon, and M. faya in Pico, São Miguel, and Terceira islands, according to the climatic changes predicted to 2100.
2. Materials and Methods
2.1. Study Area
The proposed study was applied in three distinct geographic areas: Pico, Terceira, and São Miguel islands, covering circa (ca.) 447 km2, 400.6 km2, and 744.6 km2, respectively. The use of distinct areas of an archipelago allows for the identification of the major factors underlying species distributions at different environmental scenarios and multiple spatial scales (e.g., see [62,71,72,74,75]).
The Azores archipelago (ca. 2346 km2), located in the North Atlantic Ocean and west of continental Portugal, is particularly well-studied and features rich data for invasive and native species occurrence, and environmental data [69,70,71]. The first human settlements were established in the Azores around A.D. 1440 [52]. More than 550 years of human presence has taken its toll on the local biodiversity, including the 420 endemic species [52,76]. The natural vegetation includes diverse communities, namely coastal vegetation, coastal and inland wetlands, meadows, peat bogs, and several types of native forests and scrubs. Nowadays, 81.1% of the vascular plant species are exotic, although most are not invasive [77]. In fact, before human settlement, natural forests could have covered around 75% of the land surface of the Azorean islands [78].
The replacement of the natural forest in the Azores has followed a clear temporal sequence. Human activities had restricted the native forest in most islands to areas above 500 meters above sea level (m a.s.l [78,79]). During the last decades of the 20th century, the reduction of native forest area was significant, with the clearing of large fragments, at mid and high altitude, for pasture [52,80]. Pittosporum undulatum, A. melanoxylon, or Eucalyptus globulus Labill. dominate most forest patches located in low-to mid-elevation areas. At higher altitudes, Cryptomeria japonica D. Don. dominates, along with the remaining stands of natural forests, particularly above 600 m a.s.l. [78,81,82].
2.2. Data source
2.2.1. Species Data
We assembled forest species data from the Forest Inventory of the Azores [83]. Occurrence records for studied species (P. undulatum, A. melanoxylon, and M. faya) were obtained by using all the polygons identified in the Forest Inventory. The data, provided in vector format (shapefile), were converted to raster (RST) and comma-separated values (CSV) formats for further analysis. All biological data were geocoded to a grid of 100 × 100 m in order to match with the environmental variables used (Figure 1, Figure 2 and Figure 3).
The three target species show different patterns of spatial distribution and range boundaries. Since, below 600 m a.s.l., most forest patches in Pico, Terceira and São Miguel islands are dominated by P. undulatum and A. melanoxylon (invasive plants), along with individuals remaining from former native forests, such as M. faya [71,72], these three species were selected. Furthermore, their distribution and ecology is relatively well known and stable compared with other plant taxa. Eucalyptus globulus Labill. and C. japonica are also important forest species but, in the Azores, there are mostly restricted to those areas where they were planted.
Pittosporum undulatum is a shrub or tree up to 15 m tall, and in its native range it is found in different types of Eucalyptus forest in Australia, being one of the most abundant species [84]. The canopy is pyramidal, 3–5 m in diameter, the leaves are evergreen, ovate, and with a wavy margin, and shows white-creamy perfumed flowers in loose umbellate cymes that bear one to four insect pollinated flowers [85,86,87,88,89]. It is widespread in the three islands but especially in Pico Island [71,74,90,91,92].
Acacia melanoxylon is an evergreen tree with a pyramidal or rounded dense crown, usually reaching 8–15 (20) m high [87]. In the Azores, this species has important settlements on Pico, in all the northern slopes of the island, and some scattered stands in Terceira and São Miguel islands, however its interest fell sharply with the introduction of C. japonica in the regional wood market [87].
Morella faya is an evergreen shrub or small tree, found on all the Azores islands and is frequently a dominant species in lowland and coastal native vegetation [78,79,93,94,95]. Nevertheless, it can also be found up to 600 m of altitude and even as high as 1000 m in Pico Island. In the Canary Islands and Madeira, it is found up to 1300 m [96,97]. According to the Forest Inventory of the Azores [83], M. faya is currently present in a total of 22% of the forested area in the Azores archipelago but in only 5% of that area it is the dominant species [71]. For this reason, it is fundamental to preserve M. faya and the associated flora and fauna [71,97].
2.2.2. Climate Data
The CIELO model [98,99,100] is a simple layer model based on the transformations experienced by an air mass crossing over the mountain, and simulates the evolution of an air parcel’s physical properties starting from the sea level up to the mountains [62].
Original climate data comprised 12 monthly time step variables developed for the project by the CIELO model [98,99,100] at a 100 × 100 m grid resolution. This model allowed us to obtain a high number of climatic variables associated with temperature (°C), relative humidity (%), and precipitation (mm) for the 20th and 21st centuries.
We related species distributions to climate using six variables derived from the original climate data set: Annual minimum temperature (°C), annual maximum temperature (°C), annual mean temperature (°C), annual minimum relative humidity (%), annual maximum relative humidity (%), and total annual precipitation (mm). The selection of these set of particular climatic variables was considered, since they were shown to be relevant for the study species [91,92].
Possible future climate scenarios in Azores were based on runs of the CIELO model forced with boundary conditions derived from RCPs suggested in the fifth Assessment Report of The Intergovernmental Panel on Climate Change [101]. The IPCC scenario RCP8.5 (more “severe” scenario) was chosen following a precautionary approach. This option is based on the historical and still actual trend on the greenhouse emissions, which follows the scenario that leads to radiative forcing of 8.5 W/m2 by the year 2100 [62,101,102].
2.2.3. Topographic Data
We also added topographic variables to the models, allowing the inclusion of geomorphological aspects, given their previous importance in former modelling exercises [71,74,91,92].
Topographic data were derived based on a digital elevation model (DEM) available in CIELO ([98,99,100]; see
The identification of the most suitable variables for modelling is decisive to maximize the accuracy of distribution models and their reliable projection in space and time. Table 1 shows a characterization of the ecogeographical variables used, and Table 2 shows the summary of present and future climatic variables for the three islands.
2.3. Modelling Approach
Potential changes in the spatial distribution of the target species were predicted under each climatic scenario by applying a modelling approach—BIOMOD (BIOdiversity MODelling; “biomod2”)—based in a computer platform for ensemble forecasting of species distributions, enabling the treatment of a range of methodological uncertainties in models, and the examination of species–environment relationships [103,104]. BIOMOD includes the ability to maximize the predictive accuracy of current species distributions and the reliability of future potential distributions, using different types of statistical modelling methods into different environmental conditions (e.g., climate, topographic, or land-use change scenarios) [73,104]. BIOMOD 2 ([105]; version 3.3-7.) is implemented in R ([106]; version 3.3.1; R Development Core Team, 2016) and is a freeware, open-source package.
We used a total of seven different SDM algorithms implemented in the BIOMOD 2 package, including modelling approaches with different characteristics and potentialities (Table 3): Two regression techniques, generalized linear models (GLM; [107]), generalized additive models (GAM; [108]); three machine learning methods, artificial neural networks (ANN; [109,110,111]); generalized boosted models (GBM; [112]); random forest (RF; [113]); one classification method, classification tree analysis (CTA; [114]); and rectilinear envelopes similar to BIOCLIM, namely surface range envelope (SRE; [115]). Regarding GLM and GAM, based on previous work [92], we have not included non-linear components (e.g., quadratic terms) in the modelling settings.
Most modelling techniques implemented in BIOMOD 2 require that species distribution data include presences and absences [73,104]. When data are presence-only, a possible solution is to generate random pseudo-absences [126,127]. This can be done in BIOMOD 2 using strategies of increasing complexity [104,128].
The modelling approach was divided into two parts based on pseudo-absence requirements: GLM-GAM-ANN-SRE and GBM-RF-CTA. Since those algorithms require presence and absence data, we generated pseudo-absences following the recommendations of Barbet-Massin et al. [129]. For GLM-GAM-ANN-SRE, we used 10000 random background absences across the whole study extent; and for GBM-RF-CTA we used a number of pseudo-absences equal to the number of occurrences. In order to minimize projection bias, we applied the “randomPoints” function in R package “dismo” [130,131].
As no independent data existed to evaluate the predictive performance of the models, species presences were randomly split into test data (30%) and training data (70%), and 10 replicate model runs per algorithm were performed to account for model variability [103,117].
A preliminary screening was performed to investigate which variables would be more suitable for modelling, based on model evaluation (described in Section 2.4). The models were tested in a sequential manner, starting with models developed using climatic variables only and, in a second step, calculating new models that also included topographic variables. We calculated a total of 843 models with different combinations of climatic or climatic and topographic variables—ecogeographic variables (EGVs), for each of the nine species/islands combinations. More specifically, we started with models incorporating all the climatic variables in Table 1, and progressed to calculate models including all the climatic and topographic variables in the same Table 4. During model evaluation, we defined the best variable set for each species/island combination.
2.4. Model Validation
Two evaluation measures were calculated: The true skill statistic (TSS; [132]) and the area under the curve (AUC) of the receiver-operating characteristic (ROC) plot [133]. Both evaluation measures were calculated via the proportion of the two prediction types “sensitivity” (correctly predicted presences) and “specificity” (correctly predicted absences) which take into account all four elements of the error matrix (true and false presences and absences), derived from the test data [117,134]. The TSS requires a specific threshold level and is a prevalence independent alternative to Cohen’s kappa [132]. It ranges from –1 to +1, with values above 0.4 indicating a statistically reliable model performance [132,135]. In terms of prediction accuracy, TSS values ranging from 0.2 to 0.5 were considered poor, from 0.6 to 0.8 useful, and values above 0.8 were good to excellent [136]. The ROC curve, proposed by Peterson et al. [137], was used to assess the predictive performance of the models. A good model is defined by a curve that maximizes sensitivity for low values of the false positive fraction [138]. To guarantee most accurate model predictions, AUC values lower than 0.5 were considered poor, in the range 0.5–0.7 considered fair, in the range 0.7–0.9 considered good, and excellent for values above 0.9 [139].
2.5. Model Projections
For the best models, including climatic only or both climatic and topographic variables, and for the nine species/islands combinations, we projected the respective potential distributions for the reference period (present) and for the future conditions, using the respective functions of BIOMOD 2. The two probabilistic ensemble forecasts were translated into two binary maps using the value that maximizes the AUC scores as the threshold for distinguishing presence and absence predictions. The two maps were then reclassified into four levels, according to the respective quartiles of habitat suitability.
2.6. Variable Importance
In BIOMOD 2, variable importance is calculated by shuffling a single variable of the given data set. The value returned is based on the correlation of the model predictions obtained (i) with the initial dataset and (ii) with the shuffled data set. The highest the value (i.e., closer to 1), the more influence the variable has on the model. A value of 0 assumes no influence of that variable on the model. Note that this technique does not account for interactions between the variables [105]. This calculation is accessible using the variables_importance() function.
3. Results
3.1. Summary of the Best Models
From the preliminary analysis of more than 800 models including different combinations of EGVs, the best models were selected based on TSS and AUC values. A summary of those models is provided in Table 4 and in Figure 4, Figure 5 and Figure 6.
For TSS and AUC, EGV set 1 (climatic variables only) gathered the least accurate predictions among the different environmental datasets. The models derived by using EGV set 2 (also including topographic variables) were among the best, according to the TSS and AUC mean values (Figure 4, Figure 5 and Figure 6).
Globally, model quality was highest for RF and GAM algorithms, with average TSS values ranging from 0.315 (RF) to 0.698 (RF) for EGV set 1, and from 0.399 (GAM) to 0.802 (RF) for EGV set 2; and average AUC values ranging from 0.715 (RF) to 0.915 (RF) for EGV set 1, and from 0.775 (GAM) to 0.953 (RF) for EGV set 2. Therefore, in general, models had good agreement between observed data and model predictions, especially when AUC was used. Regarding Pico Island and EGV set 2, the models for the three species showed excellent predictive performances (AM, RF, AUC = 0.928 ± 0.010, mean ± sd; PU, RF, AUC = 0.935 ± 0.003; MF, RF, 0.953 ± 0.006; Supplementary Material—Table S1). For Terceira Island and EGV set 2, the models for the three species showed good to excellent predictive performances (AM, GAM, AUC = 0.775 ± 0.013; PU, RF, AUC = 0.831 ± 0.016; MF, RF, AUC = 0.912 ± 0.016; Supplementary Material — Table S2). For São Miguel Island and EGV set 2, the models for the three species showed good to excellent predictive performances (AM, RF, AUC = 0.832 ± 0.006; PU, RF; AUC = 0.859 ± 0.007; MF, GAM, 0.903 ± 0.011; Supplementary Material—Table S3).
3.2. Relative Importance of EGVs
The importance of EGVs varied between species and islands (Supplementary Material—Figures S1, S2, S3). For P. undulatum, annual minimum temperature, total annual precipitation, and elevation were more important in Pico Island, while total annual precipitation and slope were more important in Terceira and São Miguel islands. For A. melanoxylon, annual maximum relative humidity, elevation, and slope were more important in Pico Island, annual minimum temperature and annual maximum temperature in Terceira Island, and total annual precipitation and slope in São Miguel Island. Regarding M. faya, total annual precipitation and winter hill shade were more important in Pico Island, and total annual precipitation and slope in Terceira and São Miguel islands.
3.3. Present and Future Projections
The probability data were projected into map format in order to compare the realized habitat distribution in the present conditions and the future expected distribution under the EGV set 2 (Figure 7, Figure 8 and Figure 9).
The present and future maps (Table 5) revealed changes in the species’ probability of occurrence (increase/decrease), thereby affecting species’ range sizes (gains/losses) at the study locations.
According to model projections, most species would decrease their geographic space when considering the change in the number of cells belonging to the third and fourth quartiles of probability (i.e., probability of occurrence above 0.5 that is above 500 in the maps).
Regarding Pico Island, in response to the future climate scenario, the three species showed a reduction in the probability of occurrence at the lowest elevation belt, and therefore a loss of the more adequate habitat, although with an increase in habitat quality at higher elevations, but still below the 0.5 probability value (Table 5; Figure 7). In Terceira Island, models predicted an extension of the more adequate habitat for A. melanoxylon and P. undulatum, although much lower for the latter. For M. faya, a pattern similar to what happened in Pico Island was found, with some degradation of habitat quality at the lowest elevation belt and a slight increase in habitat quality at higher elevations (Table 5; Figure 8). Still, the spatial distributions patterns of the species remained more favourable at relatively low elevations, except for A. melanoxylon, which would apparently gain habitat at the interior areas of the island. Finally, for São Miguel Island, M. faya, which presently has a very restricted distribution, is predicted to have a profound increase in its most suitable habitat, while for A. melanoxylon and P. undulatum, a pattern similar to that found in Pico Island is expected, with habitat degradation at the lowest elevation belt and some increase of habitat suitability at higher elevations, although still below a probability of 0.5 (Table 5; Figure 9).
4. Discussion
One important research topic is to understand the currently ongoing decline in the biodiversity of the ecosystems [140]. Changes in land use and climate have been identified as the main drivers of this decline [24,65]. Forest ecosystems are under pressure, as climate change may threaten species across a wide range of taxonomic groups [64,65] and the conservation of forest biodiversity depends on predicting the effects of many elements of global change. However, this is a daunting and complex task [2]. One of the greatest uncertainties in the available climate changes variables comes from the downscaling method used to refine projections at the regional scale [141].
An impressive number of new climate change scenarios have recently become available to assess the ecological impacts of climate change [40,41]. Among these impacts, shifts in species ranges analysed with SDM algorithms are the most widely studied [38,40,62].
Casajus et al. [40] suggest the use of clustering approach, k-mean clustering, which reduces the number of climate change scenarios needed to project species distributions, while retaining the coverage of uncertainty in future climate conditions. This methodology selects an objective subset of climate change scenarios and offers an efficient guidance to project species distribution through time.
Our results showed that, despite the good accuracy of the model’s predictions, there was considerable variation depending on the type of modelling approach, the species and the island. Therefore, modellers must be aware that their results will include several degrees of uncertainty, when projecting possible changes in species distribution attributed to climate change [40,131,142,143]. Moreover, in general, the use of climatic EGVs only led to a lower accuracy in model predictions, when compared with models also including topographic variables [144]. Furthermore, the association with topographic variables might have a connection with the history of land use in the target area. As an example, the present distribution of M. faya in São Miguel Island might be the result of forest clearing, leading to the limitation of species occurrence to sloped terrain or the ridges of volcanic craters, since the flatter areas were converted to agriculture [72,94,97]. A similar situation is seen with P. undulatum in São Miguel Island, where the plant has occupied native forest habitat and forest patches along water lines and sloped terrain [90], while in Pico Island, it has invaded large areas of abandoned vineyards at the lowest elevation belt. Therefore, we should stress that the previous history of land use might be reflected on the present distribution of the target species. Thus, their present distribution might not be at balance with the environment, thereby affecting future predictions largely based on the present relationship between climatic data and species distribution data, and leading to a possible underestimation of the extent of suitable habitat see [71,91].
In this context, a critical question arises. Which and how many climatic environmental/climate change scenarios are required to project species distributions, and carry out the impact analyses that cover the range of possible climate futures? Casajus et al. [40] propose an objective approach to select climate scenarios when projecting species distribution under climate change. Yet, there are no fixed guidelines when selecting the best type of climatic variables to map spatial projections, which is also contingent on the particular period of time to be considered and on the magnitude of environmental changes in the study area [145]. When environmental changes are weak, range dynamics might be determined predominantly by intrinsic population dynamics, as described by Lawton [146] and Pimm [147]. However, when environmental changes are strong, species distributional dynamics should be driven heavily by these changes [145].
Although a new climate might not directly favour invasive plant species over natives, many invasive plant species share traits that could increase their dominance in a transitioning climate scenario [2]. Moreover, new suitable areas might emerge in what were previously unsuitable or marginal habitats while previously suitable areas might be converted into unsuitable or marginal ones [145,148,149].
In our analysis, we found a clear effect of the type of modelling algorithm on model adjustment. RF and GAM showed the highest TSS and AUC values when compared to the other algorithms. One of the main advantages of GAM is that it offers a great flexibility in order to represent the relationship between the dependent variable and the explanatory variables [108]. GAM models are usually more adjusted due to their non-parametric component or semi-parametric component (i.e., “smooth functions”), and can fit nonlinear curves in a way that permits a great variety of functional forms [108,118]. On the other hand, RF is an ensemble classifier that consists of many decision trees, and generates an output class that is the mode of the classes generated by the individual trees [113]. It is one of the most accurate learning algorithms available; for many datasets, it produces a highly accurate classifier [150]. However, as this a single case study, caution is required before making specific recommendations about which algorithm to favour over others, as the algorithms that performed well here may not do so in other areas or under other circumstances [13]. In previous studies in the Azores, using MaxEnt, and targeting the same species and islands, we have found AUC values that were generally slightly below, or in particular cases, slightly above those presently obtained with RF or GAM [92]. It has been suggested that the use of multiple SDMs, which vary in complexity and statistical mechanisms, is a more robust way to assess species distributions [151,152]. One study reported a comparison of performance of five modelling approaches (GLM; GAM, RF; GBM; and MaxEnt) and concluded that none provided superior predictions in all performance criteria [153], in agreement with previous comparative studies [154,155].
Regarding potential model overfitting, an option to avoid this situation would be to use a metric such as AIC (Akaike’s information criterion) or other derived criteria (e.g., BIC, Bayesian information criterion; WAIC, Watanabe–Akaike information criterion), as reported in our previous research dealing with GLMs in a maximum likelihood or Bayesian context (see [92]). However, experts suggest that errors can occur in the estimation of the degrees of freedom (i.e., the number of parameters in the model is not the effective number of degrees of freedom in the models) when computing AIC for MaxEnt, GBM, RF, and other algorithms [156]. Therefore, we followed the guidelines suggested by Hauenstein et al. [156] and focused on cross-validation.
Although some authors consider that the predicted distributions obtained by different SDMs/algorithms are affected by a variety of factors, but still improved by combining them in an ensemble [103], we experienced some difficulties in following this guideline. In fact, we opted to base our projections in the algorithm that provided the best fit, since the projections obtained by the ensemble models often suggested extreme scenarios for species increase of decrease. This was even more accentuated for those islands/species where land use might have had a stronger effect on species distribution (e.g., M. faya in São Miguel Island). In addition, when ensemble model selection entirely relies on the AUC metric (which has been criticized; [157,158]), then the additional spatial error introduced may have a greater effect due to the inclusion model being sensitive to this error, while producing a high value of AUC [131]. In the case of our study, we further evaluated the predictive performance of the models, by randomly splitting data into test (30%) and training sets (70%), by replicating models 10 times, and also by using two evaluation metrics.
The importance of the EGVs varied between species and islands. However, elevation and slope were among the most important variables in the EGV sets. In terms of climatic variables, those related to annual precipitation values seemed to play an important part in habitat suitability. All these findings agree with previous modelling studies in Azores [71,72,74,91]. In general, the importance of EGV sets indicated that elevation and precipitation influence the ecological niche, suggesting that the studied species tend to avoid the high mountains, where temperatures are low and relative humidity and rainfall are high.
Regarding models projections and the expected changes in species distribution, our results revealed different outcomes for different species and islands. Namely, in some cases, species are predicted to clearly increase their distribution range, e.g., M. faya in São Miguel Island. This could be expected since the present distribution of this species was most likely reduced by land use changes that occurred in the past [71,97]. In contrast, for Pico Island, and particularly for P. undulatum and A. melanoxylon, a reduction in the number of cells within the third and fourth quartiles of habitat suitability was predicted. Apparently, habitat adequacy will decrease at low elevation, probably due to a temperature increase, while slightly increasing at higher elevations, but without a sufficient compensation for the habitat degradation that might occur at the lowest elevation belt.
In previous studies in Terceira and São Miguel islands, Ferreira et al. [62] found a tendency for upward shift in altitude in the suitable climate space of the analysed vascular plants, with the highest decrease occurring on coastal areas, a similar situation to what was found in the present study. However, regarding vascular plants, the mentioned researchers focused mainly on coastal species with a relatively narrow altitudinal distribution, whereas we analysed the distribution of more common and widespread species. Thus, differences in future predictions are expected when modelling species with very different ecological niches. Ferreira et al. [62] also found a possible increase in suitable climate space for the endemic Laurus azorica (Seub.) Franco, particularly in Terceira Island, which is agreement with our predictions of habitat gain for P. undulatum and A. melanoxylon in Terceira, and M. faya in São Miguel Island.
Meanwhile, prediction of range shifts under changing climatic conditions has important limitations, mainly related to the disregard of biotic interactions (e.g., [159]), dispersal limitations and intraspecific variation in niche breadth [160,161,162]. Moreover, land use might also play an important role in the present distribution of the target species [71,91], and will certainly affect their distribution in the future. Furthermore, the reliability of species distribution models depends on the quality of the used data [152], which results from historic factors that might have affected the present distribution (e.g., forest clearing, active plantation), besides the analysed climatic variables. Also, our results suggest that a potential extension of the range of an invasive species could not be the obligatory outcome of climate change, as seen for P. undulatum in Pico Island, with a possible degradation of its habitat at low elevation but a slight increase in habitat quality at high elevations. This still adds more uncertainty to the real outcome of climate change in the future distribution of native and invasive trees.
The future climate projections in Azores (PRAC 2017—Programa Regional para as Alterações Climáticas dos Açores; [163]) suggest an upward trend in global mean temperature, projections also verified by the IPCC (Intergovernmental Panel on Climate Change; [101]). However, due to the strong Atlantic influence, the future projections indicate that this increase will not be as marked as in the continental regions, namely Mainland Portugal [163]. This is mainly because of the geographical location of the islands, the greater oceanic thermal inertia, and the heat exchanges between this environment and the atmosphere [99,100,163]. Towards the end of the century, mean temperature will change with a projected increase between 1.4 and 1.9 °C for RCP4.5 and between 2.5 and 3.2 °C for RCP8.5. The increase of temperature will be more accentuated in the islands of the western group [163,164]. Changes in precipitation across the archipelago will not be uniform and some shift of precipitation from littoral to the interior will be expected [164]. Globally, for the entire surface of each island, a moderate increase in precipitation in the RCP8.5 scenario, for the period 2010–2039, is expected, especially for the central group. For the same scenario, a generalized decrease in precipitation is expected for all the islands by the end of the century, especially in the western group [164].
As seen above, the present study suggests that much attention has to be devoted to algorithm selection, differences within the same species for different habitats (i.e., islands or regions), and differences among species [131]. Peculiarities such as the previous historic land use of the target areas, might have constrained the present distribution of the target species, rendering predictions based on climate variables alone prone to potential error. Meanwhile, the inclusion of topographic or land use data in the modelling process should not be seen as way to replace climatic data, but integrated in a sequential manner, as proposed in our research. This would allow the identification of possible effects derived from previous land use, while avoiding the risk of biasing predictions under future climate, with unknown but potentially different land use. Therefore, our results emphasize the adoption of a careful study design and algorithm selection process. The uncertainties associated with climate change, as well as the unknown future of plants communities as a whole, including the response of both their indigenous and invasive components, highlight a pressing need for integrated modelling, monitoring, and experimental work to realize the ecological and geographical consequences of climate change, in order to ensure the resilience of forest ecosystems in a changing world.
5. Conclusions
Our results emphasize the need for the adoption of a careful study design and algorithm selection process. In fact, we found some relevant differences in model performance among different algorithms. In particular, RF was revealed to perform consistently well along the different islands and species. We also found that, in general, the EGVs associated with species distributions were similar to what was obtained in the previous modelling exercises. Also, models that including topographical variables performed better than models including climatic EGVs alone. This shows that other aspects, besides climate, most likely associated with historic land use changes, presently play an important role in forest species distribution, therefore potentially blurring our understanding of future range alterations associated with climate change. We also obtained some results that had not been previously anticipated, namely a global degradation of P. undulatum habitat in Pico Island, resulting from a potential loss of adequate habitat at the lowest elevation belt, which would not be fully compensated by the slight increase in habitat suitability that would be expected at high elevation. We should also mention that the future climate scenario used in this exercise (RCP8.5) mainly predicts an increase in temperature between 2.5 and 3.2 °C, with only a slight (and again uncertain) change in relative humidity and rainfall. Considering all this underlying uncertainty, the continuous but uncertain action of land use in the future, and the variation in the predictions that were obtained, depending on the island and on the target species, we conclude that the present predictions about changes in species distribution ranges associated with climate change are still profoundly uncertain.
Supplementary Materials
The following are available online at
Author Contributions
L.D.S. and L.S. conceived and designed the modelling questions, performed the modelling exercises, and analysed the data; E.B.d.A. and F.V.R. contributed with EGV information and R.B.E. revised the paper and discussion.
Funding
L Silva’s research on forests was supported by the project “FOREST-ECO2—Towards an Ecological and economic valorization of the Azorean Forest—ACORES-01-0145-FEDER-000014—Azores 2020 PO, 2016–2019”, supported by “Governo dos Açores” and FEDER. The climatic modelling performed by E B Azevedo and F V Reis was supported by the project “PRoAAcXXIs—Projecções das Alterações Climáticas nos Açores para o Século XXI, PO_ACORES-01-0145-FEDER-000037”, supported by “Governo dos Açores” and FEDER. This publication was also funded by national funds through FCT—Fundação para a Ciência e a Tecnologia in the frame of the project UID/BIA/00329/2019” and by the Azorean Government through DRCT-M1.1.a/005/Funcionamento-C-/2019 (CIBIO-A), Secretaria Regional do Mar, Ciência e Tecnologia.
Acknowledgments
We thank “Direção Regional dos Recursos Florestais” and “Governo dos Açores” for making the Azores Forest Inventory data publicly available.
Conflicts of Interest
The authors declare no conflict of interest.
Figures and Tables
Enlarge this image.Figure 1. Species distribution data for P. undulatum (PU), A. melanoxylon (AM), and M. faya (MF) in Pico Island.
Enlarge this image.Figure 2. Species distribution data for P. undulatum (PU), A. melanoxylon (AM), and M. faya (MF) in Terceira Island.
Enlarge this image.Figure 3. Species distribution data for P. undulatum (PU), A. melanoxylon (AM), and M. faya (MF) in São Miguel Island.
Enlarge this image.Figure 4. Evaluation of P. undulatum (three distinct islands) distribution models/algorithms (mean values) by the true skill statistic (TSS) and the receiver operating characteristics curve (AUC), using two different EGV sets: EGV set 1 (climatic variables) and EGV set 2 (climatic and topographic variables). The models/algorithms were: GLM (generalized linear models), GAM (generalized additive models), ANN (artificial neural networks), SRE (surface range envelope), GBM (generalized boosted regression models), RF (random forest), and CTA (classification tree analysis). TSS1: TSS of EGV set 1; TSS2: TSS of EGV set 2; AUC1: AUC of EGV set 1; AUC2: AUC of EGV set 2. Error bars correspond to the standard deviation. For more details, see Supplementary Material Tables S1–S3.
Enlarge this image.Figure 5. Evaluation of A. melanoxylon (three distinct islands) distribution models/algorithms (mean values) by the true skill statistic (TSS) and the receiver operating characteristics curve (AUC), using two different EGV sets: EGV set 1 (climatic variables) and EGV set 2 (climatic and topographic variables). The models/algorithms were: GLM (generalized linear models), GAM (generalized additive models), ANN (artificial neural networks), SRE (surface range envelope), GBM (generalized boosted regression models), RF (random forest), and CTA (classification tree analysis). TSS1: TSS of EGV set 1; TSS2: TSS of EGV set 2; AUC1: AUC of EGV set 1; AUC2: AUC of EGV set 2. Error bars correspond to the standard deviation. For more details, see Supplementary Material Tables S1–S3.
Enlarge this image.Figure 6. Evaluation of M. faya (three distinct islands) distribution models/algorithms (mean values) by the true skill statistic (TSS) and the receiver operating characteristics curve (AUC), using two different EGV sets: EGV set 1 (climatic variables) and EGV set 2 (climatic and topographic variables). The models/algorithms were: GLM (generalized linear models), GAM (generalized additive models), ANN (artificial neural networks), SRE (surface range envelope), GBM (generalized boosted regression models), RF (random forest), and CTA (classification tree analysis). TSS1: TSS of EGV set 1; TSS2: TSS of EGV set 2; AUC1: AUC of EGV set 1; AUC2: AUC of EGV set 2. Error bars correspond to the standard deviation. For more details, see Supplementary Material Tables S1–S3.
Enlarge this image.Figure 7. Present and future species distribution maps modelled for Pico Island: (a,b) P. undulatum, (c,d) A. melanoxylon, and (e,f) M. faya. Evaluation measure, AUC: Area under the receiver operating characteristic curve. EGV set 2. A probability scale from 0 to 1000 was used to facilitate map reading.
Enlarge this image.Figure 8. Present and future species distribution maps modelled for Terceira Island: (a,b) P. undulatum, (c,d) A. melanoxylon, and (e,f) M. faya. Evaluation measure, AUC: Area under the receiver operating characteristic curve. EGV set 2. A probability scale from 0 to 1000 was used to facilitate map reading.
Enlarge this image.Figure 9. Present and future species distribution maps modelled for São Miguel Island: (a,b) P. undulatum, (c,d) A. melanoxylon, and (e,f) M. faya. Evaluation measure, AUC: Area under the receiver operating characteristic curve. EGV set 2. A probability scale from 0 to 1000 was used to facilitate map reading.
General characterization of the ecogeographical variables used in the species distribution models.
| Environmental Variables | Code | Units | Data Type | Data Source |
|---|---|---|---|---|
| Annual minimum temperature | TMIN | °C | Climatic | CIELO model 1,2 |
| Annual maximum temperature | TMAX | |||
| Annual mean temperature | TM | |||
| Annual minimum relative humidity | RHMIN | % | ||
| Annual maximum relative humidity | RHMAX | |||
| Total annual precipitation | PT | mm | ||
| Elevation | DEM | m | Topographical | CIELO model |
| Aspect | ASP | ° | ||
| Slope | SLP | % | ||
| Curvature | CRV | |||
| Flow accumulation | FLA | |||
| Summer hill shade | SHS | |||
| Winter hill shade | WHS |
1 Present environmental data based on CIELO model; 2 Future environmental data based on the representative concentration pathways (RCPs) by screening and selection of ran CIELO model scenarios from the fifth Assessment Report (IPCC-AR5 2014).
Table 2Summary of present and future climatic variables for three distinct islands. TMIN: Annual minimum temperature; TMAX: Annual maximum temperature; TM: Annual mean temperature; RHMIN: Annual minimum relative humidity; RHMAX: Annual maximum relative humidity; PT: Total annual precipitation; Min.: Minimum; Max.: Maximum; Q1: First quartile; Q3: Third quartile. Present environmental data based on CIELO model; Future environmental data based on the representative concentration pathways (RCPs) by screening and selection of ran CIELO model scenarios from the fifth Assessment Report (IPCC-AR5 2014).
| Variables | Pico Island | Terceira Island | São Miguel Island | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Min. | Q1 | Median | Mean | Q3 | Max. | Min. | Q1 | Median | Mean | Q3 | Max. | Min. | Q1 | Median | Mean | Q3 | Max. | |
| Present | ||||||||||||||||||
| TMIN | 1.1 | 10.8 | 12.7 | 12.6 | 14.7 | 17.6 | 7.8 | 11.5 | 12.6 | 12.6 | 14.0 | 15.6 | 6.6 | 10.8 | 12.2 | 12.0 | 13.3 | 15.7 |
| TMAX | 7.4 | 16.1 | 17.9 | 17.9 | 19.8 | 22.9 | 13.5 | 16.9 | 18.0 | 18.0 | 19.3 | 20.9 | 13.3 | 17.3 | 18.6 | 18.4 | 19.7 | 22.2 |
| TM | 4.3 | 13.5 | 15.3 | 15.2 | 17.2 | 20.2 | 10.7 | 14.2 | 15.3 | 15.3 | 16.7 | 18.2 | 10.0 | 14.0 | 15.4 | 15.2 | 16.5 | 18.9 |
| RHMIN | 67.8 | 81.1 | 87.5 | 86.9 | 93.2 | 100.0 | 77.6 | 84.3 | 89.4 | 89.1 | 93.3 | 100.0 | 73.7 | 84.9 | 89.0 | 89.0 | 93.1 | 100.0 |
| RHMAX | 82.6 | 93.7 | 97.5 | 95.9 | 98.9 | 100.0 | 89.2 | 95.4 | 97.9 | 97.2 | 99.3 | 100.0 | 85.7 | 96.2 | 98.3 | 97.4 | 99.2 | 100.0 |
| PT | 974.0 | 1397.0 | 2301.0 | 2510.0 | 3346.0 | 6817.0 | 1126.0 | 1168.0 | 1621.0 | 1831.0 | 2280.0 | 4258.0 | 1027.0 | 1105.0 | 1382.0 | 1531.0 | 1864.0 | 2988.0 |
| Future | ||||||||||||||||||
| TMIN | 4.5 | 13.7 | 15.5 | 15.4 | 17.4 | 20.4 | 10.8 | 14.3 | 15.4 | 15.4 | 16.8 | 18.4 | 9.6 | 13.7 | 15.1 | 14.8 | 16.1 | 18.5 |
| TMAX | 10.7 | 18.9 | 20.6 | 20.6 | 22.5 | 25.7 | 16.3 | 19.6 | 20.7 | 20.7 | 22.0 | 23.6 | 16.3 | 20.1 | 21.4 | 21.2 | 22.4 | 25.0 |
| TM | 7.6 | 16.3 | 18.1 | 18.0 | 20.0 | 23.1 | 13.5 | 17.0 | 18.0 | 18.0 | 19.4 | 21.0 | 12.9 | 16.9 | 18.2 | 18.0 | 19.3 | 21.8 |
| RHMIN | 68.7 | 82.0 | 88.1 | 87.4 | 93.4 | 100.0 | 78.5 | 85.2 | 90.1 | 89.7 | 93.6 | 100.0 | 74.2 | 85.5 | 89.5 | 89.4 | 93.4 | 100.0 |
| RHMAX | 82.5 | 93.5 | 97.4 | 95.8 | 98.9 | 100.0 | 89.2 | 95.3 | 97.9 | 97.2 | 99.3 | 100.0 | 85.7 | 96.2 | 98.2 | 97.4 | 99.2 | 100.0 |
| PT | 934.5 | 1332.0 | 2280.0 | 2511.0 | 3397.0 | 7050.0 | 1086.0 | 1121.0 | 1579.0 | 1796.0 | 2243.0 | 4372.0 | 988.0 | 1064.0 | 1337.0 | 1488.0 | 1816.0 | 3021.0 |
Modelling algorithms used.
| Modelling Technique | Description | References |
|---|---|---|
| GLM: Generalized linear models | This is a generalization of ordinary linear regression, allowing for response variables that have normal (i.e., Gaussian) or non-normal error distribution models, such as the binomial and Poisson, through the use of a link function. Akaike’s or Bayesian information criteria (AIC or BIC) are used to select the most informative and parsimonious model. | [107,116,117] |
| GAM: Generalized additive models | A non-parametric extension of GLMs, where the linear predictors correspond to smooth nonlinear functions of the predictive variables. Typically use the same underlying distributions for the response variable as GLMs. | [108,118] |
| ANN: Artificial neural networks | Models are based on a stepwise progression of multivariate and univariate analyses, which are versatile in extracting information out of complex data, and which can be effectively applicable to classification and association. | [119,120] |
| GBM: Generalized boosted models | The individual models consist of classification or regression trees. In an iterative process, a final model is built by progressively adding trees while re-weighting the data poorly predicted by the previous tree. | [112,117] |
| RF: Random forest | Consists in an ensemble of unpruned classification or regression trees, created by using bootstrap samples of the training data and random feature selection in tree induction. | [113,121,122] |
| CTA: Classification tree analysis | Operates by recursively parsing the training observations on a binary splitting measure applied to explanatory variables, such as spectral responses. | [114,123,124] |
| SRE: Surface range envelope | It is a presence-only approach, based on the environmental conditions corresponding to the occurrence data, allowing the definition of the environment where a species can be found. Uses the extreme percentiles, as recommended by Nix or Busby. | [105,115,125] |
Sets of ecogeographic variables (EGV) that produced models with highest fit in Pico, Terceira, and São Miguel islands using BIOMOD 2. PU: P. undulatum; AM: A. melanoxylon; MF: M. faya. The modelling approach was initiated by calculating models with all the climatic or all the climatic and topographic variables described in Table 1. During model evaluation, we determined the best variable sets for each species/island combination, which are reported herein.
| Variables | Code | Pico Island | Terceira Island | São Miguel Island | |||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| PU | AM | MF | PU | AM | MF | PU | AM | MF | |||||||||||
| EGV set | EGV set | EGV set | EGV set | EGV set | EGV set | EGV set | EGV set | EGV set | |||||||||||
| Climatic | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | 1 | 2 | |
| Temperature | |||||||||||||||||||
| Annual minimum temperature | TMIN | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | ||
| Annual maximum temperature | TMAX | + | + | + | + | + | + | + | + | + | + | + | + | + | + | ||||
| Annual mean temperature | TM | + | + | + | + | + | + | + | + | ||||||||||
| Humidity | |||||||||||||||||||
| Annual minimum relative humidity | RHMIN | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | ||
| Annual maximum relative humidity | RHMAX | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | ||
| Precipitation | |||||||||||||||||||
| Total annual precipitation | PT | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | + | ||
| Topographic | |||||||||||||||||||
| Elevation | DEM | + | + | + | + | + | + | ||||||||||||
| Aspect | ASP | + | |||||||||||||||||
| Slope | SLP | + | + | + | + | + | + | + | + | ||||||||||
| Curvature | CRV | + | |||||||||||||||||
| Flow accumulation | FLA | + | + | + | |||||||||||||||
| Summer hill shade | SHS | + | + | + | |||||||||||||||
| Winter hill shade | WHS | + | + | + | + | + | |||||||||||||
Values of respective quartiles for habitat suitability maps: P. undulatum, A. melanoxylon, and M. faya based on the EGV set 2. Table shows the percentage of projected habitat space gain (+)/loss (-). P: Present; F: Future.
| Species | Island | Scenario | Q1 | Q2 | Q3 | Q4 | Gain (+)/Loss (−) |
|---|---|---|---|---|---|---|---|
| P. undulatum | Pico | P | 7599 | 2802 | 2891 | 31238 | −85 |
| F | 20246 | 19107 | 2994 | 2183 | |||
| Terceira | P | 23822 | 6035 | 4453 | 5506 | 7 | |
| F | 22073 | 7083 | 4971 | 5689 | |||
| São Miguel | P | 46414 | 7954 | 6687 | 13208 | −46 | |
| F | 40194 | 23234 | 9470 | 1365 | |||
| A. melanoxylon | Pico | P | 21668 | 3604 | 2735 | 16523 | −92 |
| F | 40939 | 2098 | 231 | 1262 | |||
| Terceira | P | 39604 | 150 | 38 | 24 | 174 | |
| F | 28501 | 11145 | 170 | 0 | |||
| São Miguel | P | 42701 | 8734 | 8332 | 14496 | −6 | |
| F | 35152 | 17637 | 16498 | 4976 | |||
| M. faya | Pico | P | 29802 | 1785 | 1740 | 11203 | −17 |
| F | 28708 | 5091 | 4952 | 5779 | |||
| Terceira | P | 25506 | 3628 | 2928 | 7754 | −32 | |
| F | 25024 | 7494 | 5483 | 1815 | |||
| São Miguel | P | 73236 | 413 | 270 | 344 | 5534 | |
| F | 28465 | 11204 | 10069 | 24525 |
© 2019 by the authors.