Habitat availability and environmental preference drive species range shifts in concordance with climate change

A progressive increase in air temperature is recognized as the most important mechanistic driver of species range shifts. However, only a few studies have simultaneously considered the influence of both extrinsic and intrinsic mechanistic drivers; there are still no studies on the roles of extrinsic and intrinsic drivers that regulate such species changes. We investigated how species will shift their geographical ranges to cope with future climate change and analysed the relative importance of the mechanistic drivers in governing species range shifts.

To infer the responses of species to climate change, several community metrics such as changes in distribution area and optimum colonization position of populations have been proposed (Lenoir & Svenning, 2015). Furthermore, changes in the distribution sizes of species in the cool (i.e., leading edge) and warm edges (i.e., trailing edge) have been reported recently (Chen et al., 2011;Parmesan et al., 1999;Rosalino et al., 2019). Although there are various directions of species shift, such as poleward, equatorward, westward, eastward, upward (or uphill) and downward (or downhill), in response to progressive climate change (Devictor et al., 2012;Domisch et al., 2013;Palmer et al., 2015;VanDerWal et al., 2013), here, we focused on the most significant direction, that is north-south direction, which has a clear thermal gradient as species distribution strongly depends on thermal conditions. Although progressive climate change is widely accepted to have major effects on global and regional biodiversity, only a few studies have attempted to quantify latitudinal range shifts of species across a large spatial scale, such as the continental scale (Bush, Nipperess, Duursma, et al., 2014;Jump, Matyas, & Penuelas, 2009).
Odonata is an ideal model taxon to investigate the effects of climate change across a large spatial scale by using various factors, such as tropical evolutionary origin of odonates that potentially limits distribution by temperature, a long history of scientific research, and extensive records on the ecology (distribution data) and biology (species traits) of odonates Collins & McIntyre, 2015;Hassall, 2015;Hassall & Thompson, 2008).
Because of their strong ability to fly long distances, most species of odonates tend to shift geographical ranges, rather than persist in original habitats through phenotypic plasticity or evolutionary adaptation, to cope with climate change (Dudaniec, Yong, Lancaster, Svensson, & Hansson, 2018;Parmesan et al., 1999), but a few species of damselflies have been shown to locally adapt (Dudaniec et al., 2018). This is of particular importance for studies incorporating predicted species distribution because the projection approaches, such as species distribution models (SDM), consider that a species can completely shift to new suitable habitats without dispersal constraints (Bush, Nipperess, Duursma, et al., 2014;Collins & McIntyre, 2015;Domisch et al., 2013;Li et al., 2015Li et al., , 2016Simaika et al., 2013;Thuiller et al., 2011).
Several potential driving mechanisms have been identified to interpret the changes in odonate distribution patterns, and temperature change, in most cases, is thought to be the primary driver (Chen et al., 2011;Corser, White, & Schlesinger, 2015;Li et al., 2014;Pateman, Hill, Roy, Fox, & Thomas, 2012;Sunday, Bates, & Dulvy, 2012). In addition, odonate species are typically restricted to the following three fundamental habitat types based on their adaptation to water velocity: lotic (running waters), lentic (standing waters) and mixed habitats (Hof, Brändle, & Brandl, 2006). Thus, temperature and water velocity are considered as potential driving factors that structure the community changes in odonates. Besides appropriate food, shelter and mates in suitable habitats for species survival, habitat availability and breadth have been considered important factors determining species range shifts (Mair et al., 2014;Morueta-Holme et al., 2013). Furthermore, trait-based (intrinsic) drivers such as species dispersal ability, body size, the length of wingspan, and the duration of flight have been reported as important factors (Angert et al., 2011;Auer & King, 2014;Fitt, Palmer, Hand, Travis, & Lancaster, 2018;Jaeschke, Bittner, Reineking, & Beierkuhnlein, 2013;Urban, 2015) because species traits have been documented to affect responses to climate change (Pacifici et al., 2017). These drivers do not act individually, but interact with each other (Morueta-Holme et al., 2013). For example, fragmented habitats may reduce the probability of species dispersal to new areas (Corser et al., 2015;Estrada et al., 2015). Therefore, both extrinsic and intrinsic drivers, including their synergistic effects, should be simultaneously considered to evaluate the changes in species range shifts.
Here, we used a comprehensive approach to simultaneously incorporate multiple community metrics and mechanistic drivers. We investigated how species will shift their geographical ranges to persist with future climate change and analysed the relative importance of the mechanistic drivers in governing species range shifts. We collected a group of commonly recorded odonate species (105 high-frequency species) across Europe and estimated the importance of the mechanistic drivers for each community metric. Our three specific predictions were as follows: (a) The overall distribution areas of odonate species will increase over time with the colonization of leading edges and extraction of trailing edges, as odonates are relatively warm dwellers that can gain more habitats under climate change conditions. (b) The increasing migration rates of community metrics will reduce over time as species range shifts are predominantly affected by temperature changes, and the increase in air temperature will reduce in the future.
(c) Besides environmental preferences (i.e., temperature and water velocity), habitat preferences will be considerably important. Overall, extrinsic drivers will override intrinsic drivers (i.e., species dispersal capacity) in interpreting spatio-temporal changes in odonates, as physical habitats with suitable thermal and water velocity conditions and food resources determine species survival.

| Data collection
We compiled the occurrence records of odonates across 16 countries in South, West and North Europe (Figure 1), where the ecology of odonates has been extensively studied. The following two databases were used to extract the occurrence records: (a) the occurrence records of 492 310 odonates (from 91 651 sites) were obtained from the Global Biodiversity Information Facility (GBIF; http://www.gbif.org; obtained on February 20, 2015). However, the number of records from France and Italy was relatively low, when compared with their areas, and (b) consequently, 5 137 and 5 538 records were extracted based on digital species distribution maps for France (http://www.libel lules.org) and Italy (http://www.odona ta.it), respectively. In total, 502 985 records from 102 326 sites and 207 F I G U R E 1 Geographical locations of (a) 105 high-frequency odonate species in Europe, (b) species richness in the present decade; species richness in (c) 2030s, (d) 2050s and (e) 2080s under RCP 2.6 (representative concentration pathway); species richness in (f) 2030s, (g) 2050s and (h) 2080s under RCP 4.5; and species richness in (i) 2030s, (j) 2050s and (k) 2080s under RCP 8.5 taxa were collected. Then, two criteria were used to filter the collected dataset. Firstly, all taxon lists were filtered to the taxonomic identification level, and only taxa identified to the species level were retained. Secondly, we ensured that a minimum number of occurrences (i.e., 100) was available for each species to maintain a high performance in SDM (Collins & McIntyre, 2015). After filtering with these two criteria, a new dataset was created with the selected 489 543 occurrence records from 100 278 sites ( Figure 1) and for 105 species (  and (c) emissions in RCP 8.5 continue to rise throughout the 21st century (IPCC, 2014). Prior to the analyses, highly correlated variables (Pearson's correlation coefficient |r| > .75) were excluded to reduce the redundancy of variables. Twenty-three environmental variables (elevation, slope, aspect, LC01-11, bio01 = annual mean temperature, bio02 = mean diurnal range, bio03 = isothermality, bio05 = maximum temperature of the warmest month, bio07 = annual temperature range, bio08 = mean temperature of the wettest quarter, bio09 = mean temperature of the driest quarter, bio12 = annual precipitation, and bio15 = precipitation seasonality) were finally selected for the following steps (Table S2).

| Species distribution models
With the present geographical records of each species and the selected 23 environmental variables in all four time periods, ensemble SDM was used to predict potential distributions of species in the present and future decades (2030s, 2050s, and 2080s) by using BIOMOD2 package (Thuiller, Georges, & Engler, 2015) in R 3.0.2 (R Core Team, 2018). All 10 algorithms described in BIOMOD2 were used in the modelling processes, namely artificial neural network, classification tree analysis, flexible discriminant analysis, generalized additive model, general boosting model, generalized linear model, maximum entropy, multiple adaptive regression splines, random forest and surface range envelop. For each algorithm, threefold cross-validation and 10 000 pseudo-absences were used, and the prevalence value was set as 0.5.
Prior to model projection, the occurrence records of a given species were randomly split into the following two sets: training (70%) and testing (30%) sets. The training set was used to calibrate the models, and the testing set was used to evaluate the models. True skill statistic (TSS) was used to evaluate predictive accuracy because TSS has been shown to be superior in quantifying the performance of SDMs with presence prediction data (Allouche, Tsoar, & Kadmon, 2006). TSS is an integrity measure that incorporates sensitivity (true presence) and specificity (true absence), which ranges from 0 (no discrimination) to 1 (perfect discrimination); only models with TSS ≥ 0.6 were selected for the subsequent analyses, in order to select only good algorithm models with a high predictive accuracy for the ensemble model. The selected models were then summarized in an ensemble model to produce more robust predictions (Araújo & New, 2007). Finally, distribution probability maps were converted into binary maps using a cut-off value that minimizes the differences between sensitivity and specificity (Liu, Berry, Dawson, & Pearson, 2005).
To account for variability and uncertainty in predicted species occurrence, instead of using all the 100,278 sampling sites, which were non-randomly distributed for the subsequent analyses, 50,000 random sites were selected across the study region; all analyses with the predicted values were performed based on these 50,000 random sites. The random sites were selected using the package sp (Pebesma et al., 2015) in R 3.0.2 (R Core Team, 2018).

| Potential driving factors
The species temperature preference index (STPI) was used as a probe to detect the most suitable temperature for a given species (Devictor et al., 2012;Li et al., 2016), and it was quantified as the annual mean air temperature conditions (i.e., bio01 in the WorldClim data) over the observed sites (Table S1). The preferred speed of present was not available for the selected 105 odonate species, and we therefore developed one water velocity index as an alternative for species flow preference in the larval stage. Based on the habitat list for a given species provided by the International Union for Conservation of Nature (IUCN; http://www.iucnr edlist.org), we assigned still habitats (e.g., lake, pool, pond, swamp and marsh) a value of 0, whereas the flowing habitats (e.g., river, stream, creek, spring, canal and channel) were assigned a value of 1. Then, the water velocity preference of a given species could be calculated as the arithmetic mean of the assigned values over the values of the observed habitats. Using this approach, the water velocity of the three typical habitat types was 0.88 ± 0.17 (lotic or running waters; mean ± standard error; n = 24), 0.12 ± 0.13 (lentic or standing waters; n = 44) and 0.29 ± 0.09 (mixed waters; n = 37), indicating the high reliability of this approach.
Species habitat availability was estimated as the proportion of species occurrence records in a given land-cover type among all land-cover types containing that species multiplied by the total area of the land-cover type in the study region. For species distributed in more than one land-cover type, habitat availability was summed across all species distribution land-cover types. Species habitat breadth was defined as the number of habitats for each species listed in the IUCN.
Information on direct dispersal ability of the odonate species is not available. Alternatively, body morphological parameters such as the length of wingspan, the body length of adults, and the duration of flight (DragonflyPix; http://www.drago nflyp ix.com) were used as proxies for species dispersal potential (Angert et al., 2011

| Statistical analyses
To test the first prediction, that is, the distribution areas of odonate species will increase over time with the colonization of leading edges and extraction of trailing edges, four community metrics were used.
The first metric, the change in distribution area, was estimated as the difference in proportional distribution area between pairwise time windows of the present decade and future, for example, between the present decade and 2030s, between 2030s and 2050s, and between 2050s and 2080s. Besides the distribution area, three colonization metrics were incorporated in our study. Among them, the change in optimum poleward colonization position was defined as the median migration distance (i.e., the position where the median of latitudinal values is located), the leading edge was defined as the area with the most northern community (i.e., sites with ≥99.9 percentile of latitudinal values), and the trailing edge was defined as the area with the most southern community (i.e., sites with ≤0.1 percentile of latitudinal values).
To test the second prediction, that is, the increasing migration rates of community metrics will reduce over time, the entire study period was split into three time windows, namely present decade-2030s, 2030s-2050s, and 2050s-2080s. The annual rates of changes in community metrics and air temperature were estimated for each time window. The analysis of variance (ANOVA) was then performed to evaluate the differences in the rates of these five metrics between each pairwise time windows. Tukey's multiple comparison test was performed when there was a significant difference (α = 0.05).
To test the third prediction, that is, species environmental preferences (i.e., temperature and water velocity) and habitat preferences (i.e., habitat availability and breadth) will override dispersal traits in governing species range shifts, general linear models were used to quantify changes in the rates of four community metrics (difference between 2080s and present decade) against six potential driving factors and their interactions. Particularly, nonlinear regressions were used for two remarkable nonlinear plots, for example, habitat availability−migration rate of the trailing edge (logarithm regression) and the duration of flight−migration rate of the trailing edge (exponential decay regression). Then, an information-theoretic approach was used to identify the best models for interpreting changes in each community metric. Model averaging was performed when the difference in delta AICc (corrected Akaike information criterion) between a given sub-model and the top-ranked sub-model was less than 2; otherwise, the top-ranked sub-model was selected as the best model using MuMIn package (Bartoń, 2015). All these statistical analyses were carried out using R 3.0.2 (R Core Team, 2018).
Furthermore, the species turnover rate between the present decade and 2080s was calculated according to the method reported by Peterson et al. (2002): where T is the turnover rate; E is local extirpation; C is local colonization; and S is species richness in the present decade.

| Model performance
The overall performance of models was good, with an average TSS of = 0.80, 0.82 and 0.81 (range: 0.63-1.00; Table S1) under RCP 2.6, RCP 4.5 and RCP 8.5, respectively. Thus, the ensemble projections for the 105 species were considered in the subsequent analyses.

| Patterns of species richness and turnover rate
The total species richness of the 50 000 randomly selected sites was predicted to increase until 2050s and then decrease until regions (logistic regression; Adj. R 2 = .31-.42, p < .0001; Figure 3b).
The changes in species composition between the present decade and 2080s presented an inverted bell curve along the latitudinal gradient for the species turnover rate under RCP 2.6 and RCP 4.5 (Gaussian regression; Adj. R 2 = .44-.57, p < .0001; Figure 3c). However, the species turnover rate did not present a clear trend under RCP 8.5 ( Figure 3c).

| Trends in community and temperature metrics
The  Table S3). The air temperature was projected to increase over time, at 0.27°C decade −1 under RCP 2.6 and 0.59°C decade −1 under RCP 8.5 (Figure 2f;  (Figure 2f; Table S3).

| Best models and driving factors
All the six candidate drivers were included in the best fitting models (Table 1, Table S4). Among them, STPI was the most important driver that structure the changes in all community metrics (importance value = 1.00); water velocity preference was selected by the models  Table S3. RCP refers to representative concentration pathway The linear regression analysis also showed that the wingspan was not a significant driver, as it did not present a significant correlation with any community metric (Table S5). STPI was significantly negatively correlated with the change rate of distribution area (Adj. R 2 = .13-.28, p < .05; Figure 4a, Table S5) and with the migration rate of the trailing edge (Adj. R 2 = .08-.26, p < .05; Figure 4i, Table S5) under all the three emission scenarios, but it positively correlated with the migration rate of the optimum position (Adj. R 2 = .10, p < .05; Figure 4e, Table S5) and with the leading edge under RCP 8.5 (Adj. R 2 = .03, p < .05; Figure 4h, Table   S5). Water velocity preference was also significantly negatively correlated with the changes in the distribution area under all the three emission scenarios (Adj. R 2 = .03-.09, p < .05; Figure 4b, Table S5) and positively correlated with the migration rate of the optimum position under RCP 8.5 (Adj. R 2 = .04, p < .05; Figure 4f, Table S5). Habitat availability was significantly positively correlated with the changes in the distribution area under RCP 2.6 (Adj. R 2 = .06, p < .05; Figure 4c, Table S5) and with the migration rate of the trailing edge under all the three emission scenarios (Adj. R 2 = .05-.11, p < .05; Figure 4j, Table S5). Habitat breadth was significantly negatively correlated with the migration rate of the optimum position under RCP 2.6 (Adj. R 2 = .03, p < .05; Figure 4g, Table S5) and with the migration rate of the trailing edge under RCP 4.5 (Adj. R 2 = .03, p < .05; Figure 4k, Table S5).
The duration of flight was significantly negatively correlated with changes in the distribution area under RCP 4.5 and RCP 8.5 (Adj. R 2 = .03-.05, p < .05; Figure 4d, Table S5) and with the migration rate of the trailing edge under all the three emission scenarios (Adj. R 2 = .05-.07, p < .05; Figure 4l, Table S5).

| Changes in the community metrics
Poleward shifts in the odonate species have been reported in various regions. For example, odonates in the UK (37 species) shifted by 6.8 km/year over the past four decades (Hickling, Roy, Hill, & Thomas, 2005), whereas a 6.4 km/year expansion in odonates was observed in South Europe between 1988 and 2006 (Grewe, Hof, Dehling, Brandl, & Brändle, 2013). On the basis of the projected species changes, six odonate species in Europe are expected to shift 5-14 km/year by 2035 (Jaeschke et al., 2013), whereas 191 stream invertebrate taxa are projected to shift 6.5-9.2 km/year across Europe by 2080s . When compared with these results, the selected 105 odonate species were predicted to exhibit relatively lower poleward migration rate, with 4.5 km/year in the optimum colonization position under the intermediate emission scenario, namely RCP 4.5 (Figure 2b). This migration rate was lower than the change rate of isotherms (8.0 km/year under RCP 4.5; Figure S1).
Thus, climatic debt is expected in our study, where the odonate species cannot catch up the poleward movement of isotherm (Devictor et al., 2012).

F I G U R E 3
Patterns of (a) species richness of odonates in the present decade (quadratic regression), (b) richness difference between 2080s and the present decade (logistic regression), and (c) species turnover rate between 2080s and the present decade (Gaussian regression) along the latitudinal gradient in Europe. RCP refers to representative concentration pathway TA B L E 1 Best models for explaining changes in odonate community metrics between the present decade and 2080s in relation to species temperature preference index (STPI; °C), water velocity preference (Velocity), habitat availability (Habitat; km 2 ), habitat breadth (Breadth), the duration of flight (Duration; month), and the length of wingspan (Wingspan; mm) by using general linear models Beyond the optimum colonization position, changes in the leading and trailing edges are particularly important to reveal the effects of climate change on community structure. A community's poleward range shifts are driven by either species expansion at the leading edge, species retraction at the trailing edge, or their combination effects (Hickling et al., 2005;Lenoir & Svenning, 2015;Thomas & Lennon, 1999). According to three principal categories (march, lean and crash) of species range shifts in response to climate change (Breshears, Huxman, Adams, Zou, & Davison, 2008), we believe that our case is the same as "march," where a community's range shift is influenced the combined effect of species colonization at the leading edge and species extraction at the trailing edge; this confirms our first prediction.
However, the total odonate species pool was predicted to unmatch  Figure   S1). Our findings also suggested that species at the leading edge were expected to expand slower than species at the trailing edge. This is consistent with the results of several studies; for example, Hickling et al. (2005) found that odonate species in southern Britain showed greater poleward shifts than those in northern Britain. In contrast, marine fish were found to have a high migration rate at the leading edge at a global scale (Poloczanska et al., 2013), whereas among 94 bird species in Finland the northern species (leading edge) had stronger shift than southern species (trailing edge) (Virkkala & Lehikoinen, 2014

| Mechanism underlying changes in the odonate species pool
Species range shifts primarily depend on thermal niche (Li et al., 2014;Palmer et al., 2015;Pateman et al., 2012), and our findings support this theory. Generally, the geographical ranges of species shift quicker in areas with intensive warming effects (Chen et al., 2011). Because of the expected reduction in the rate of increase in air temperature in the future, we observed a decrease in the rate of change in distribution areas, confirming our second prediction. In contrast, Urban (2015) performed a meta-analysis of 131 published studies and suggested that the extinction risk of species will be globally accelerated in the future, acknowledging the importance of temperature change. We believe the extinction risk of species is an accumulative process, and greater effects on community metrics may accumulate in response to future warming processes. As expected, the importance of water velocity was observed in the models of the leading and trailing edges.
The importance of environmental changes (i.e., temperature and water velocity) is likely due to the fact that environmental variables are the primary limiting factors that influence the decision-making process of the investigated species. When suitable temperature conditions and water velocities are available in potential habitats, odonates will undoubtedly move towards and occupy these potential habitats owing to their strong dispersal abilities. After colonizing the new habitats, the population will develop further when the supply of food resources (by means of habitat availability) is sufficient and the ecological niche (by means of habitat breadth) is broad enough (Hof et al., 2012;Mair et al., 2014). Therefore, despite the strong effects of environmental changes on the poleward movements of odonate species, we cannot ignore other potential drivers such as habitat preferences. Overall, habitat preferences were found to be important for migrations to the leading and trailing edges. This is because the expansion of habitat provides more routes and space for species expansion, whereas the fragmentation of habitat undoubtedly constrains species range shifts (Estrada et al., 2015). In our study, both environmental and habitat preferences explained species range shifts; this is consistent with the results of Morueta-Holme et al. (2013), who studied 85,000 plant species range changes across North and South America. However, we did not obtain strong evidence to indicate that the intrinsic drivers (i.e., the duration of flight and length of wingspan) are important, confirming our third prediction. This is consistent with the findings of several studies, which indicated that species traits are poor in predicting geographical changes in response to climate change (Angert et al., 2011;Corser et al., 2015;Grewe et al., 2013;Mair et al., 2014). Interestingly, we found that the length of wingspan was relatively important under RCP 8.5, indicating that once the temperature increases dramatically, the restrict of intrinsic driving mechanism will be uncovered.
Nevertheless, the importance of extrinsic drivers indicates that habitats with suitable temperature, water velocity and food conditions are a prerequisite for the survival of organisms; only when the constraint of extrinsic forces reduces, the intrinsic force may significantly affect species range shifts. In contrast, some studies have highlighted the importance of species dispersal traits in determining geographical range shifts (Auer & King, 2014;Estrada et al., 2015;Hof et al., 2012;Pacifici et al., 2017).

| Limitations of the study
Although we confirmed our three predictions, there are considerable unexplained variations in the best models, and this is consistent with the findings of several similar studies (Angert et al., 2011;Karssene, Chammem, Khorchani, Nouira, & Li, 2017;Li et al., 2016;Mair et al., 2014). This may be because we examined only six drivers.
Numerous factors were not considered in our present models, and they may be important, for example, some physical (e.g., precipitation, deforestation, habitat degradation and landscape dispersal barriers) and biological factors (e.g., other types of species traits, species interaction, the force of prey, colonization and extinction, and intraspecific trait variation) (Morueta-Holme et al., 2013;Parmesan et al., 1999). However, obtaining such variables is usually difficult in a broad-scale region. Parameter-based analyses of numerous species across a large spatial scale often interpret a small amount of species variations in response to environmental changes (Mair et al., 2014).
There were some limitations to our study, and these are similar to those of studies on climate change. Firstly, our study was conducted using data from specific European countries, where extensive odonate species records are available, but the complete geographical ranges of some species were likely not covered. Secondly, our modelled future distributions may be overoptimistic for some species because we did not account for all threats to odonates, such as dispersal barriers and habitat fragmentation, which may affect some species colonization routes and water velocity (Estrada et al., 2015). Thirdly, species with only a few presence records are typically excluded to maintain high model performance, which may limit the application of our approach for rare species, especially species of great conservation concern (Simaika et al., 2013). Fourth, our results may be also affected by the spatial resolution used (5 km × 5 km in our case), as climatic conditions vary greatly within a grid of such a resolution (Collins & McIntyre, 2015). Despite these potential limitations, our results are comparable between community metrics at our studied spatial extent, and we are confident that these limitations do not substantially alter our findings.
In our study, 105 records of odonate species were used, and such extensive data across a large spatial scale are vital for understanding biodiversity changes in response to different types of environmental changes. Importantly, we found that the drivers of species range shifts, such as climate change, vary with time, and extrinsic drivers (i.e., environmental and habitat preferences) are important to explain such responses. It is important to note that for conservation management, the rate of species range shifts is highly variable and likely to alter species compositions and ecological functions, and therefore, raise new demands for conservation planning processes. Nevertheless, the relatively low explanatory power of the selected best models implies that the choice of suitable potential drivers, albeit often neglected, is crucial for studies on large-scale range shifts to cope with climate change. We propose that follow-up studies for understanding species range shifts should consider incorporating more factors mentioned above for higher spatial resolutions, to improve the predicative accuracy of models and increase biological realism.

DATA AVA I L A B I L I T Y S TAT E M E N T
The environmental data and locations of odonates can be accessed from the websites mentioned in the Materials and Methods.