Fuzzy sets allow gaging the extent and rate of species range shift due to climate change

The recent modification of species distribution ranges in response to a warmer climate has constituted a major and generalized biogeographic change. The main driver of the shift in distribution is the disequilibrium of the species ranges with their climatic favourability. Most species distribution modelling approaches assume equilibrium of the distribution with the environment, which hinders their applicability to the analysis of this change. Using fuzzy set theory we assessed the response to climate change of a historically African species, the Atlas Long-legged Buzzard. With this approach we were able to quantify that the Buzzard’s distribution is in a latitudinal disequilibrium of the species distribution with the current climate of 4 km, which is driving the species range northwards at a speed of around 1.3 km/year, i.e., it takes 3 years for the species to occupy new climatically favourable areas. This speed is expected to decelerate to 0.5 km/year in 2060–2080.


Results
The climatically favourable areas for the Atlas Long-legged Buzzard include Northern Africa, Sicily, the Middle East, and the Southwest of the Iberian Peninsula. The latter area is the only continental European region where the climate is highly favourable (Fig. 2). In Africa, highly favourable areas are located in Morocco, the Mediterranean coasts of Algeria, Tunisia, Egypt, and two inland areas of Algeria (Tassili n' Ajjer national park and El Menia oasis) (Figs. 1b, 2). In the Middle East, favourable areas are in the margins of the Mediterranean and Caspian seas, the Arabian Peninsula, Cyprus, Iran, and Turkey. These areas are currently occupied by its sister subspecies, the Asian Long-legged Buzzard (Figs. 1a, 2). The climatic favourability model included six variables (Table 1). A small diurnal temperature range, cool temperatures during the wettest season, and scarce precipitation in the driest month are the best predictors for Atlas Long-legged Buzzard breeding because these variables were the first to be entered in the stepwise procedure. In addition, under the Wald test, these variables have the most weighting in the model ( Table 1).
The model had high classification power (Table 2), higher sensitivity than specificity, and a high over-prediction rate (OPR > 0.8). Discrimination capacity was high (AUC > 0.8) and highly significant (p = 6.516 × 10 -51 ). The Pch value of the model was remarkably high ( Table 2). The Hosmer and Lemeshow calibration test (Fig. 3) showed that differences between expected and predicted values were nonsignificant (p > 0.05) over the whole range of expected probability values. Good calibration is very difficult to obtain when there are many operational geographic units (OGUs) because small relative differences between observed and predicted probabilities generate significant values in the statistic 51 . www.nature.com/scientificreports/ Although forecast favourable areas did not differ substantially between the two periods of time, they differed from the present climatic favourability pattern (Fig. 4). Current breeding areas in North Africa are expected to have lower favourability values in the future, but other areas are expected to increase their favourability values in Sicily, Sardinia, and the Iberian and Arabian peninsulas (Figs. 2, 4).
There was a low mean uncertainty of less than 0.1 associated with climate scenarios (Fig. 5). For the period 2061-2080, there is an increase in uncertainty values in some areas, such as on the Atlantic coast of France and the northern coast of the Black Sea, where in the period 2041-2060 the uncertainty values are much lower. The highest uncertainty values are in the northern half of the Iberian Peninsula, the northern edge of the Caspian Sea, the coast of the Balkan Peninsula, and Turkey, particularly for the period 2061-2080. Low uncertainty values are found in the south of the Iberian Peninsula and the coasts of north-western Africa.
A northward disequilibrium of 0.039 degrees of latitude was detected between the latitudinal barycentre and the climatic favourability barycentre of the OGUs with reported breeding of the buzzard (Table 3). A latitudinal degree is roughly equivalent to 111.12 km, and so this disequilibrium was equivalent to 4.3 km. Applying the same equivalence, the FD 20-60 rate value was 1.3 km/year. If this rate is assumed to be similar to the current one, then  Table 1. The map was created using ArcMap software (ArcGIS 10.4.1) https ://deskt op.arcgi s.com/es/arcma p/. Table 1. Variables entered in the logistic regression model by the forward-backward step-wise selection process, ranked by their order of entrance. β are the coefficients in the logit function, Wald is the Wald's statistics value (representing the relative importance of the variable in the model) and p the significance of the coefficients. Codes of the variables are the same as in Table 5.  www.nature.com/scientificreports/ the latitudinal disequilibrium of 4.3 km between distribution and climatic favourability is equal to a temporal disequilibrium of a little more than 3 years. The FD  value was lower at 0.6 km/year. Table 4 shows the fuzzy logic indicators of other expected effects of climate change for each possible scenario and the ensemble forecast at each period of time. Increment values are low (I < 0.08) but generally positive: however, some negative values are expected in lower CO 2 emissions scenarios (2.6, 4.5, and 6.0). Overlap values of present and future favourable areas are high (O > 0.78). Maintenance values for currently favourable areas are not expected to be complete (0.92 at most) although they are generally high (M > 0.87). In all cases, shift values are lower than 0.11 but higher than 0.08, indicating that between 8 and 11% of currently favourable areas for breeding are expected to be lost but could be replaced by new opportunities for breeding in emergent favourable areas elsewhere.

Discussion
Five of the six significant predictive variables included in the favourability model were associated with climate (Table 1). In addition, the only non-climatic variable, altitude, was the last to be included in the model and had the lowest weight in the model under the Wald test criterion (Table 1) (Table 1). Especially in Europe, climate change is expected to increase the temperatures during day and night and to reduce the number of cold nights as well as increasing the frequency and duration of heat waves 7 , which will get the European region a warmer and drier environment. Extreme low temperatures, unlike extreme hot temperatures, are predicted to be less frequent during daily and seasonal timescales 7 , so despite the fact that stability in daily temperatures is considerably hard to foresee, climate change predictions suggest a low increment in the diurnal range of temperatures, except in the southern part of the study area. The annual mean precipitation will likely decrease in Europe, although in some humid regions it is predicted to slightly increase under some scenarios 7 . Under these conditions the possibility of establishment in Europe will likely increase in the coming years. The positive coefficient associated with altitude (Table 1) may be due to the fact that this species mainly nests in cliffs 52 . However, in the Iberian Peninsula it nests in trees 19,35 , which shows the plasticity of the species when choosing a nesting substrate. This behavioural change in nesting requirements could be mediated by a compensation of climatic favourability in areas where cliff-nesting areas are scarce or already occupied by other species. Such plasticity in nesting behaviour has also been demonstrated for other species 53 . This topographical variable could also be indicative of climatic conditions that are closely correlated with altitude 54 . As our models were based on fuzzy set theory, they were able to measure the spatial and temporal disequilibrium between the buzzard breeding area and the climatic favourability for breeding. The mean climatically favourable condition occurs 4.3 km northwards of the mean latitudinal range of the buzzard breeding areas. The mere existence of this spatial disequilibrium demonstrates that, in order to be sound, favourability distribution models do not have to assume that the species must be in equilibrium with their environment, as models based on niche theory have customarily assumed 55,56 . In fact, a spatial disequilibrium is needed for the climatic conditions of the favourable areas to attract the buzzard and lead to changes in its breeding range 49,50,57 . The FD  value suggests that the buzzard range is not in temporal equilibrium with the current climate by a margin of about 3 years. This result suggests that, if climatic conditions stop changing, the buzzard range could attain equilibrium with the climate in approximately 3 years. Otherwise, the effect of the climatic disequilibrium will continue to attract the breeding range to the favourable adjacent areas. This latitudinal range shift is expected to decrease to 0.6 km/year from 2060 to 2080, suggesting that the colonization of the Western Palearctic will decelerate according to current climatic projections.
Fuzzy logic operations (IOMS framework, Table 4) indicate that the climatic favourability for the Atlas Long-legged Buzzard is shifting rather than broadening. An initial indicator of this change is that the Increment values are low and even negative in some cases. This result would not be expected if the expansion of the range was northward, as this would include the maintenance of favourable areas in their southern range and the addition of new favourable areas in the north 13 . Maintenance values of only about 90% are mainly explained by the favourability loss in the southern limit of the current distribution, as in the Sahel region (Fig. 4). However, the main indicator is that the Shift values are more than zero, and practically identical to the part of the current  www.nature.com/scientificreports/ favourability area that is not expected to be maintained (Table 4). This result indicates that the loss of climatically favourable areas in the south are expected to be compensated for by the appearance of new favourable areas in the north, particularly in the coast of northern Morocco, the southern half of the Iberian Peninsula, as other study previously stated 2 , and the Mediterranean islands of Sardinia and Sicily (Fig. 4). This loss may effectively push the breeding distribution of the buzzard northward. High Overlap between present and future climatic favourability models indicates that the current and expected favourable areas are geographically close 23 , which could facilitate the movement of buzzards to the newly available zones that are currently unoccupied. There is a high degree of similarity between all predicted future scenarios. The forecasts for the 2061 to 2080 period have the highest uncertainty, which is normal when the forecast time period is further away from the present 7 .
It should be noted that the Asian Long-legged Buzzard is currently occupying areas in the Middle East that have been identified as highly favourable for the Atlas Long-legged Buzzard 58 . This situation is likely preventing the establishment of the Atlas Long-legged Buzzard because of possible competitive exclusion [59][60][61] . It may also be the main reason for the high over-prediction rate and high Pch of the climatic favourability model for the Atlas Long-legged Buzzard ( Table 2, Fig. 2), and that the range shift is more likely to be northwards rather than eastwards. This forecast is supported by empirical data, such as the increasing number of records of the Atlas Long-legged Buzzard in northern Spain [62][63][64] and Portugal 65 and its successful breeding since 2009 in the south of the Iberian Peninsula 19 .
However, other factors could interfere in the northwards colonization process, by delaying or even stopping the Atlas Long-legged Buzzard dispersion. For instance, hybridization with the Common Buzzard (B. buteo) has already been detected in Pantelleria Island (Italy) and southern Spain 35,66 . These areas represent new contact zones where these closely related species currently meet. Hybridization in these areas is favoured because the Atlas Long-legged Buzzard is currently rare and, hence, its choice of mate is restricted. If the offspring were infertile, hybridisation would be acting as a synecological barrier against the expansion of the Atlas Long-legged Buzzard into Europe 35,67 . Thus, some new climatically favourable areas could remain unoccupied by the Atlas Long-legged Buzzard in the future. It is therefore relevant to track its northward expansion into Europe by monitoring favourable areas on the distribution fronts, such as the Strait of Gibraltar, and in those areas detected as highly favourable in this study.
Many other African birds have already experienced or are currently experiencing a similar expansion pattern in the Iberian Peninsula. Some examples are the White-rumped Swift (Apus caffer), the Little Swift (Apus affinis), and the Common Bulbul (Pycnonotus barbatus), which is the most recent addition to the European avifauna and is already breeding 2,68-70 . Moreover, in the last decade, there has been an increasing number of reports of other African birds such as Rüppell's Vulture (Gyps rueppellii), Lanner Falcon (Falco biarmicus), Cream-coloured Courser (Cursorius cursor), Moussier's Redstart (Phoenicurus moussieri), and House Buntings (Emberiza sahari) 33,36,62,71 . Other species, such as the White-backed Vulture (Gyps africanus) and Bateleur (Terathopius ecaudatus) have also been recently reported in Europe for the first time 33,34 . These reports may indicate that European Mediterranean areas, and particularly the Iberian Peninsula, should be prepared for the establishment of new African fauna in the near future.
These species could use these areas as stepping stones to expand into the rest of Europe. For example, this pattern was followed by the Black-shouldered Kite (Elanus caeruleus), which started breeding in Tarifa (southern Spain), then moved northwards, and became established in France in 2013 69,[72][73][74] . Some typical Mediterranean birds are also moving northwards and reaching Central Europe, where they are breeding for the first time. This is happening, for example, in Switzerland, where typical Mediterranean species, such as the Short-toed Snakeeagle (Circaetus gallicus) and the European Bee-eater (Merops apiaster) have settled recently and experienced population growth 75 . These observations indicate that a change in species composition is already occurring in Western Europe due to climate change.
Favourability models have proven to be useful tools to identify new potential areas for the arrival and establishment of new individuals in the near future 76,77 . The practical usefulness of these tools must be coupled to a sound theoretical background 78 . Our results reflect the interaction between the Atlas Long-legged Buzzard and its environment within a dynamic biogeographic framework. The favourability function is well established in fuzzy set theory and it also represents the response function of the species to the environmental conditions 44,54,79,80 . This aspect may constitute a contact point between fuzzy set theory and niche theory. In our view, this contact point requires adding notions of graduality to the niche concept, such as those included in the Maguire 81 niche concept 57 . More theoretical work is needed on the relationships between fuzzy set theory and niche theory. This would allow forecasting that is simultaneously practical and theoretically comprehensive and, thus, useful for conservation management and the scientific understanding of the biogeographic processes underlying the species response to climate changes.

Methods
Study area. The study area was the European, Asian, and African zones from 20° 00′ W to 60° 00′ E and from 09° 30′ N to 70° 00′ N (Fig. 1b). This area, which comprises the Western Palearctic and surrounding areas, fully covers the current breeding territories of the Atlas Long-legged Buzzard, as well as the western distribution of the Asian Long-legged Buzzard. Climatic heterogeneity is high, covering sub-tropical, desert, Mediterranean, Atlantic, and tundra climates [82][83][84] . To analyse the distribution of the subspecies, the study area was converted into 1-degree latitude × 1-degree longitude grid cells (n = 3839) using the Create Fishnet and Intersect tools from ArcGIS. These cells were used as operational geographic units (OGUs) 85,86 .

Species distribution data.
There is such little information on the Atlas Long-legged Buzzard that it is considered a gap in the knowledge of Western Palearctic birds 42  World atlases were used to obtain the distribution data to determine the breeding area of the Atlas Longlegged Buzzard 38,42,58 . This information was updated using the IUCN red list 89 , e-Bird (https ://ebird .org/speci es/ lolbu z1), and some personal records from Morocco and Spain 2 . We identified the OGUs where the species had been reported as breeding at least once from 1980 to 2020. Thus, we obtained a binary target variable representing breeding/not breeding at each OGU (breeding: n 1 = 146, not breeding: n 0 = 3693, Fig. 1b). Table 5) related to topography and climate between 1950 and 2000 were used in the biogeographic modelling procedure. These variables were digitized in raster format at a resolution of 1-km 2 pixels. Values of these variables at each OGU were obtained by averaging the values of the 1-km 2 pixels within them using the zonal function of ArcGIS 10.4.1 software.

Predictor variables and future scenarios. A set of 21 environmental variables (list and sources in
Expected future values of the climatic variables were obtained for the periods 2041 to 2060 and 2061 to 2080 90 (https ://world clim.org/). Four different Representative Concentration Pathways (RCP) were used to represent the extreme (2.6 and 8.5) and intermediate (6.0 and 4.5) trajectories of future CO 2 emissions, with different effects on precipitation and temperature 7 . We also used two different Global Circulation Models (GCM; HadGEM2-ES and NorESM1-M) to consider other sources of uncertainty regarding the future climate of the study area 23,91,92 . This process resulted in eight sets of expected values of the climatic variables for each period of time.
Model for the present. We performed a logistic regression of the binary target variable on each environmental variable separately. This is a supervised machine learning approach that assesses the predictive power of each variable according to the significance (α) of the score test of the corresponding regression. Multicollinearity between these variables was reduced by calculating the Spearman correlation coefficients r between them. For each pair of variables with r > 0.8, only the one with the highest individual predictive power was retained 93 . On the basis of this subset of predictors, the False Discovery Rate 94 (FDR) was evaluated to control the increase in type I errors (i.e. familywise error rate), and therefore the likelihood of obtaining false significant results when a large number of variables are used in the modelling process 95 . Using the Benjamini and Yekutieli 96 procedure for all forms of dependency between test statistics, only variables significantly associated with the distribution of the subspecies with α < 0.05 under an FDR value of 0.05 were accepted in subsequent modelling procedures. Table 5. Variables selected to model the Atlas Long-legged Buzzard distribution grouped by environmental factor. Sources: (a) Ref. 116 ; (b) calculated from Alti with ArcGIS software; (c) Ref. 90 . Exc. is the procedure that excluded the variable, being a Spearman's correlation value, b FDR analysis and c step-wise selection process. www.nature.com/scientificreports/ A comprehensive model for the current probability of breeding at every OGU according to their climatic conditions was obtained by multivariate forward-backward stepwise logistic regression of the target variable on the remaining subset of predictors. This procedure starts with a model with no predictor variable (i.e. the null model), which yields a constant probability of breeding at each OGU equal to the prevalence of the OGUs where breeding was reported in the whole OGU dataset. Then, a significant combination of predictors (y or logit) is built by adding the variables that provide the most significant contribution to the model obtained in the previous step 51,97 . If no predictor variable significantly adds to the null model then no environmental model is produced. Variables with an overall, broad scale, predictive power are entered first in the modelling procedure, while those that add significant nuances to the previous model are entered in subsequent steps. To avoid redundancy, variables that do not significantly add to the predictive power of the final model are not included, as their effect, if any, is effectively included in the model via correlated variables.
The effect of prevalence on the resulting probability values was discounted to extract the pure response of the buzzard to environmental conditions 57,81 . This was done by obtaining favourability values (F) using the Eq. (1) 98 .
where n 1 and n 0 are the number of OGUs where breeding was reported or not reported, respectively, e is the Euler's number, and y is the logit resulting from the performed logistic regression.
Favourability values range from 0 (minimum favourability) to 1 (maximum favourability). A local favourability value of 0.5 indicates that the local probability of breeding of this subspecies is the same as its prevalence in the study area, i.e., is the probability expected by a null model unaffected by environmental predictors, where breeding is neither favoured nor unfavoured by the environment. Hence, favourability refers to the degree to which the environmental conditions favour the breeding of the subspecies 50,99 , being a favourability value of 0.5 the threshold separating favourable from unfavourable areas. The concept of favourability differs from that of suitability in: (1) there is only one way of obtaining favourability from probability and prevalence, while different algorithms yield idiosyncratic suitability values; (2) favourability values are interpretable in absolute terms, indicating the extent to which probability of local presence differs from that expected by chance in the whole sample, while suitability values are relative; and (3) favourability values represent the degree of membership of the localities to the fuzzy set of sites with conditions that are favourable for the species, which enables the application of fuzzy logic operations to distribution modelling 50 . However, given the continuous and fuzzy character of favourability 50 , the use of a favourability value of 0.5 as a cut-off point for crisply distinguishing favourable from unfavourable areas is not sufficiently informative 100  Model assessment. The relative weight of each variable in the final model was assessed using the Wald test 102 . The discrimination capacity of the resulting model was evaluated using the Area Under the Receiver Operating characteristic Curve (AUC) 103,104 , which has an associated significance value. We used Cohen's Kappa Index to measure the degree to which the favourability of the OGUs with reported breeding or no reported breeding in the dataset was higher or lower than 0.5, respectively (Kappa is described as the proportion of specific agreement, whose values range from − 1 to + 1) 105 . As described by Fielding and Bell 106 and Barbosa et al. 107 , we also applied a set of classification measures, whose values range from 0 to 1. These measures were sensitivity (the conditional probability of OGUs with reported breeding being classified as favourable), specificity (the conditional probability of OGUs with no reported breeding being classified as unfavourable), correct classification rate (CCR: the conditional probability of correctly classified OGUs), the over-prediction rate (OPR: the proportion of OGUs with no reported breeding in the area with favourability higher than 0.5), and the under-prediction rate (UPR: the proportion of OGUs with reported breeding in the area with favourability lower than 0.5). Good classification performance is shown by high Kappa, sensitivity, specificity, and CCR values and low over-and under-prediction rate values. Model calibration was assessed using the Hosmer-Lemeshow test 51 , where nonsignificant values (p > 0.05) indicate a good fit between predicted and observed probabilities (i.e. a well calibrated model). The test was performed by dividing the probability range of the model into 10 bins of equal range and checking that each bin included at least 15 OGUs and at least 5 OGUs in which reported breeding was expected in all of them, requirements that should be fulfilled for reliable model calibration 51 .
We evaluated the current potential for dispersion of the Atlas Long-Legged Buzzard using the potential change factor (Pch) in OGUs with breeding. This factor is the ratio between the OGUs with favourability value higher than 0.5 and OGUs with reported breeding. Values lower or higher than 1 indicate higher or lower potential regarding the actual breeding distribution of the buzzard, respectively 77 .
Projection to future climate scenarios. Future climatic favourability values (F f ) were obtained by replacing in the logit (y) of Eq. (1) the present values of the climatic variables with the expected future values according to each RCP and GCM and for each future period of time 23,108,109 . This process resulted in eight expected climatic favourability models for each period. An ensemble forecasting of the models was obtained for each period of time by calculating the mean values of the eight future climatic favourability models at each OGU. The uncertainty of the ensemble forecasting was computed using fuzzy set theory 43 , given that favourability values may be considered to be the degree of membership in the fuzzy set of areas favourable for buzzard breeding 98 . Thus, the favourability function is the membership function that assigns each OGU their degree of membership value 98 . The uncertainty of the ensemble forecasting at each OGU was computed as the difference at each OGU between

Assessment of latitudinal variation.
In biogeography, the barycentre (B) of a variable in the range of a species is the centre of gravity of a species distribution along the gradient of that variable [111][112][113] . When applied to geographic coordinates, the barycentre is indicative of the geographic centre of the distribution range 114 . The latitudinal barycentre of the OGUs with reported breeding of the Atlas Long-legged Buzzard (B bre ) was obtained as the arithmetic mean of the latitudinal values at the centre of the OGUs. For the present and future models, the latitudinal barycentres of climatic favourability (B F ) were obtained by weighting the latitude with the favourability using Eq. (2): where La and F are the latitude (in decimal degrees) and climatic favourability values, respectively, at each OGU. For the present model, the climatic favourability barycentre was applied to the following OGUs: (1)   where c(X) is the cardinality of the fuzzy set X: that is, the sum of the degrees of membership of all the OGUs (i.e. favourability values) in the fuzzy set X, which is a measure of the size of the fuzzy set; Fp and Ff are the fuzzy sets of present and future favourable areas for the buzzard, respectively; F f ∩ F p is the fuzzy intersection between future and present favourability defined by the minimum of the two favourability values in each OGU, representing the present climatic favourability conditions that are expected to persist in the future; F f ∪ F p is the fuzzy union between future and present favourability, defined by the maximum of the two favourability values in the OGU; and Min is the minimum of the two values in square brackets. Positive values in I indicate a gain in overall climatically favourable areas, whereas negative values indicate a net loss of climatic favourability. O, M, and S range from 0 to 1. O values closer to 1 are obtained when a large proportion of favourable areas are shared in the present and future models, whereas a value of zero indicates total separation between climatically favourable areas at present and in the future. An M value of 1 indicates that the present favourable areas will be completely maintained in the future projections. The lower the M values the less the currently favourable areas License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.