Biodiversity is heavily influenced by ongoing climate change, which often results in species undergoing range shifts, either poleward or uphill. Range shifts can occur provided suitable habitats exist within reach. However, poleward latitudinal shifts might be limited by additional abiotic or biotic constraints, such as increased seasonality, photoperiod patterns, and species interactions. To gain insight into the dynamics of insect range shifts at high latitudes, we constructed ecological niche models (ENMs) for 57 Odonata species occurring in northern Europe. We used citizen science data from Sweden and present-day climatic variables covering a latitudinal range of 1,575 km. Then, to measure changes in range and interactions among Odonata species, we projected the ENMs up to the year 2080. We also estimated potential changes in species interactions using niche overlap and co-occurrence patterns. We found that most Odonata species are predicted to expand their range northward. The average latitudinal shift is expected to reach 1.83 and 3.25 km y–1 under RCP4.5 and RCP8.5 scenarios, respectively, by 2061–2080. While the most warm-dwelling species may increase their range, our results indicate that cold-dwelling species will experience range contractions. The present-day niche overlap patterns among species will remain largely the same in the future. However, our results predict changes in co-occurrence patterns, with many species pairs showing increased co-occurrence, while others will no longer co-occur because of the range contractions. In sum, our ENM results suggest that species assemblages of Odonata—and perhaps insects in general—in northern latitudes will experience great compositional changes.
Human-influenced climate change has had major impacts on biodiversity (e.g., Walther et al. 2002, Parmesan 2006). The global average temperature has increased by 0.85°C since 1880 and is predicted to increase globally another 0.3 to 0.7°C by 2035 (IPCC 2014). Furthermore, owing to Arctic amplification, the largest temperature increases are observed at high latitudes in the northern hemisphere (Serreze and Barry 2011, IPCC 2014). Hence, at these latitudes, large decreases in seasonality are predicted (Xu et al. 2013) which in turn may affect species distributions (Root et al. 2003). Ultimately, climate change is predicted to increase species extinction risk (Urban 2015), but it will also affect phenology (Cohen et al. 2018), physiology (Pörtner and Farrell 2008), abundance (Carey and Alexander 2003, Johnston et al. 2013), and distribution (Velásquez-Tibatá et al. 2013, Vieira et al. 2018). Given that high latitudes in the northern hemisphere are predicted to be extra sensitive to climate change (IPCC 2014), knowledge about the effects of climate change on northern organisms is especially important.
To keep pace with climate change, species could adapt to new environmental conditions, or go through shifts in geographic distribution to track optimal conditions (Parmesan and Yohe 2003). Although the potential of species to tolerate change and to evolve may have been underestimated (Nadeau and Urban 2019), numerous studies have documented ongoing poleward or altitudinal range shifts (Hickling et al. 2005, Wilson et al. 2007, Barbet-Massin et al. 2012). For dispersal-limited species, most range shifts are predicted to occur in the form of contractions of species ranges, as high-latitude and montane habitats become reduced compared to mid-latitude and low-elevation regions, or because a range may contract faster at the trailing edge than it can expand at the leading edge (Lenoir and Svenning 2015). While assessments of vulnerability to climate change have been extensively documented in the birds and mammals (Pacifici et al. 2015), the most intense responses are expected in ectotherms, such as insects, which are particularly sensitive to ambient temperature change (Wilson et al. 2007, Sinclair et al. 2016). Insects represent the vast majority of animal species and biomass, and fill many niches within ecosystems (Prather et al. 2013, Scudder 2017). Given that this group of organisms is already undergoing dramatic population declines (Wagner 2020), it is of the utmost importance to understand the impact that future climate change may have on the survival of insect species.
Ecological niche models (ENMs) are increasingly used to infer how environmental variables are shaping species distributions (Elith and Leathwick 2009, Peterson et al. 2011). ENMs are spatially explicit analytical tools that relate species occurrence data with the bioclimatic variables to infer a range of suitable conditions for species and generate habitat suitability maps. Coupled with climate change scenarios, ENMs are often used for projecting future change in species distributions (Thuiller 2004), making it possible to identify areas of particular interest for conservation (Markovic et al. 2014), anticipate emerging threats from vector-borne diseases (Leta et al. 2018), and many other applications. Most ENMs include predictors which are ecologically relevant for the target species, such as temperature or precipitation data, topography, or land cover (Elith and Leathwick 2009). In addition, the use of ENMs to predict distributional changes comes with implicit assumptions regarding dispersal abilities. Both unlimited- and no-dispersal scenarios can be relevant and often run in parallel (Collins and McIntyre 2015), though attempts to model cost-weighted dispersal kernels have been performed (Bush et al. 2014). Another important assumption in these models is that a species' ecological niche does not change through time (niche conservatism). However, the phenotypic plasticity (Valladares et al. 2014) and the disruption of current or the emergence of future biotic interactions (Urban et al. 2013, Fitt and Lancaster 2017) are likely to alter ecological niches, which are hard to anticipate. Moreover, range shifts might be limited by additional non-climatic abiotic constraints increased at high latitudes, such as greater seasonality or lower solar insulation (Spence and Tingley 2020). Nevertheless, the assumption of the Eltonian noise hypothesis suggests that biotic interactions are averaged out at large spatial scales (Soberón and Nakamura 2009). Thus, ENMs integrating various climatic variables can still useful for predicting future species distributions under different climate change scenarios without explicitly considering the biotic interactions when they are largely unknown (Jenkins et al. 2020).
Odonates (dragonflies and damselflies) are valuable indicators that can be used to track climate change impacts on entire insect communities. This is partly because of the substantial dispersal abilities of the adult stage. While the adults are terrestrial, their larval stages occur in the freshwater environments. Thus, the niche of odonates might be driven by both aquatic and terrestrial conditions. They are widely distributed and the adult stage is conspicuous, thus, promoting studies and data collection (Bybee et al. 2016). However, extensive documentation of dispersal is still lacking in odonates but field and molecular studies indicate that some species have reduced dispersal abilities (a few hundred meters covered in an adult lifetime in Coenagrion mercuriale, Watts et al. 2004) while others such as long-distance migrants or arid-zone species, can span hundreds of kilometers (May 2013, Suhling et al. 2017). In the same population, dispersal distance can vary greatly with most individuals staying near their emergence site with few long-distance dispersers (Keller and Holderegger 2013). The colonization of new areas might stem from both recurrent small-scale movements and rare long-distance dispersal events. Currently, there is a growing interest in the application of ENMs to odonates (Patten et al. 2015, Pires et al. 2018, Rodríguez-Tapia et al. 2020). But, multispecies ENMs are still rare for odonates (Collins and McIntyre 2015) and represent a very low portion of arthropod ENMs published (Mammola et al. 2021). For instance, Hickling et al. (2005) found an average northward range shift of 74 km between 1960 and 1995 for 37 nonmigratory British odonates with an increase in range size for all but two species. More recently, Termaat et al. (2019) showed that temperature increase has benefited most odonate species regardless of their temperature preference in various European countries, namely, Sweden. However, neither of these two studies used an ENM approach. Moreover, another study, which was conducted in Germany using occupancy models, confirmed the positive effect of temperature preference as well as the key role of species habitat on distribution trends (Bowler et al. 2021). While ENMs have been constructed for two closely related northern European damselfly species (Wellenreuther et al. 2012), future projections are still lacking for high-latitude species in this region (Collins and McIntyre 2015, but see Li and Park 2020). Although regional-scale ENMs usually do not cover the entire distribution of species, they continue being useful in predicting distributional change and guiding future surveys (El-Gabbas and Dormann 2018), especially in the high latitudes neglected by the past ENM efforts (Collins and McIntyre 2015). Within odonates, a variety of interspecific interactions can be involved from intraguild predation and competition at both adult and larval stages (Moore 1964, Wissinger 1992) to territorial and mating interference (McEachin et al. 2022). Given the wealth of biotic interactions possible within odonates, it is not feasible to consider them a priori as potential factors shaping odonate distributions. Nonetheless, inferring from ENMs probabilities of species interactions and potential changes in these probabilities because of the ongoing climate change can be instrumental in anticipating more fundamental compositional changes. Indeed, ENM approaches may be essential in helping us understand how odonates, and insects in general, will respond to climate change in the northern latitudes.
While the range of boreal species is expected to contract dramatically, northern high-latitude areas are also expected to gain temperate species that will undergo northward range shifts (Langham et al. 2015). Despite the great dispersal abilities of odonates, their tropical evolutionary history and the associated physiological constraints may limit poleward colonization into areas of high seasonality and low-solar insulation (Hassall and Thompson 2008). Here, we used opportunistic occurrence data, generated by citizen science, spanning a latitudinal range of 1,575 km throughout Sweden (northern Europe), to model ecological niches for several odonate species, in order to address the following questions: 1) what is the present-day geographic distribution of odonate species richness in Sweden; 2) how will odonate distributions change by 2080, assuming niche conservatism and unlimited dispersal; 3) will niche overlap and co-occurrence patterns among species change over time as a result of climate change possibly leading to new interspecific competition or to the loss of existing interactions; and 4) whether traits, such as temperature preference, habitat, and phylogeny, affect niche overlap and vulnerability to climate change.
Methods
Occurrence Data
Opportunistic geo-referenced observations of 66 odonate species occurring in Sweden were extracted from the Swedish Species Observation System, Artportalen ( https://www.artportalen.se/), which is a freely accessible reporting system used by the citizen scientists from all around the country. Our analysis is thus limited to Sweden. Despite this limitation, this is a comprehensive data set, collected at a fine spatial resolution, and it covers a large geographic area, extending latitudinally over 1,575 km, with many data entries and extensive coverage. We extracted ca. 200,000 odonate records collected over a 30-yr period, from 1991 to 2020. However, we only considered data collected between 2006 and 2020, a period during which 2,500 or more records were reported annually. The retained records also represented a period when reporting was more geographically uniform, with the earlier records being more biased toward the south of the country. Citizen scientists are responsible for the accuracy of their records (date, location, and species identification) that are regularly validated by the experts. Records with low resolution (> 1 km) and outliers with doubtful or inaccurate species identification based on uploaded photos were discarded. Data were not filtered for life history stages or other indications of autochthony (oviposition, larvae, and exuviae). Geographic coordinates initially available in a country-specific coordinate reference system (SWEREF99 TM, EPSG:3006) were converted into WGS84 (EPSG:4326).
Environmental Predictors
To determine habitat suitability, we used a set of 16 ecologically and physiologically relevant climate predictors (ENVIREM, Title and Bemmels 2018) and also one topographic predictor (altitude). The climate predictors relate to various components of odonate life cycle (see Supp Table S1 [online only] (nvac056_suppl_supplementary_material.docx) for complete list) such as the length of the growing season for larval development, some levels of humidity or continentality that are important components of adult stage preferences and tolerance, and temperature seasonality that can constrain larval development and voltinism. The 16 ENVIREM climate predictors were generated from monthly temperature (maximum, minimum, and average) and precipitation climatologies, as well as solar radiation averaged over recent decades (1979–2013) available in the CHELSA v1.2 data set (Karger et al. 2017). This period was considered to reflect present-day climatic conditions and provide the finest spatial resolution available, with a rectangular grid of 0.0083 degree (ca. 1 km × 0.3 to 0.5 km depending on the latitude within the study area). The same variables were generated for future periods (2061–2080) under two greenhouse gas concentration trajectories (RCP4.5 and RCP8.5, intermediate and worst-case scenarios, respectively). The future climatologies we used for these two scenarios originate from seven different global circulation models (GCMs) from the CMIP5 generation (ACCESS1-0, CCSM4, CMCC-CM, CNRM-CM5, GFDL-ESM2G, HadGEM2-CC, and MPI-ESM-MR). These were chosen because they have been shown to yield satisfactory predictions in Europe (McSweeney et al. 2015) and low levels of interdependence (Sanderson et al. 2015). Altitude data were obtained from the EU-DEM v1.1 ( https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1) and the resolution was upscaled to fit that of the climate variables.
Environmental Niche Models
We used environmental niche models (ENMs) to predict changes in species distribution dependent on the future climate change. For accurate niche-modeling estimates, records should be evenly distributed geographically. However, despite removing data from early collection years, there was still unevenness present in the data, with more records from the south compared with the northern latitudes. In order to reduce spatial autocorrelation due to such uneven reporting efforts across the study area, we spatially thinned the occurrence data set for each species with the R package ‘spThin’ (Aiello-Lammens et al. 2015). This made it possible to keep a maximal number of records outside a neighborhood with a radius of 10 km around each data point. Nine of the 66 species with less than 25 records after the thinning step were not included for niche modeling. These species are either recent colonizers (Aeshna affinis, Anax ephippiger, Crocothemis erythraea, Lestes barbarus, and Lestes viridis), or are elusive restricted-range species (Nehalennia speciosa, Somatochlora sahlbergi, and Sympecma paedisca), or have probably been accidentally introduced (Sympetrum pedemontanum). Thus, the input data consisted of presence-only occurrences for 57 species (the number of records for each can be found in Supp Table S2 [online only] (nvac056_suppl_supplementary_material.docx)), with 17 environmental predictors. In order to assess potential multicollinearity among the environmental predictors, we ran the function removeCollinearity from the ‘virtualspecies’ R package with the default threshold set at 0.7. Despite the existence of multicollinearity among the 17 environmental predictors (seven groups of intercorrelated predictors), all were included in the ENM to increase model performance (i.e., limit the overestimation of range compared to ‘true’ range), and decrease the uncertainty associated with the choice of predictors (Beaumont et al. 2005) and the risk of missing relevant predictors for at least some species. However, we note that including all 17 predictors prevents any reliable analysis of predictor contribution to the model (Sillero and Barbosa 2021). Prior to performing ENM projections based on the future climatologies, we used the mobility-oriented parity (MOP) metric (Owens et al. 2013) on a random sample of 10% of all the geographic grid cells in the study area. We used MOP to estimate Euclidean distances in the multivariate environmental space between the present-day and future ranges. This metric helps to identify areas where future values of predictors are outside of their present-day range and thus models make extrapolations.
ENM algorithms were used to compute habitat suitability maps from occurrence data and environmental predictors. For each species, we ran four types of machine learning algorithms from the stacked species distribution model R package (‘SSDM’, Schmitt et al. 2017): multivariate adaptive regression splines (MARSs), generalized-boosted regression models (GBMs), random forests (RFs), and support vector machines (SVMs). The four algorithms have been shown to be robust to multicollinearity issues (Dormann et al. 2013) and yielded models with some of the best performances among the algorithms available in the package assessed over a selection of species ( Supp Table S3 [online only] (nvac056_suppl_supplementary_material.docx)). The transformation from habitat suitability to binary presence/absence was carried out using the threshold that maximized the true skill statistics (TSS; Allouche et al. 2006). TSS corresponds to the sum of sensitivity and specificity minus one (the sensitivity and specificity are the proportions of correctly predicted presences and absences, respectively). Model fit was assessed using the area under the curve (AUC) of the receiver operating characteristic (ROC) plot (Manel et al. 2001) and kappa statistic (Allouche et al. 2006). Each of the four algorithms was run twice on present-day variables and the eight models were then assembled (and weighed according AUC) into an ensemble ENM using the ensemble_modelling function from the ‘SSDM’ package (Schmitt et al. 2017). Model cross-validation was performed by splitting the occurrence data set into a training set (70% of the records) and an evaluation set (the remaining 30%). The present-day ensemble ENM were projected over the two future climate scenarios on separate GCMs, with the assumption of the unlimited dispersal. All the models and mapped predictions where visually assessed and compared with input occurrences to detect any obvious mismatches between data and models ( Supp Fig. S1 [online only] (nvac056_suppl_supplementary_material.docx)). From the three sets of 57 ensemble ENM for present day and future scenarios, a stacked SDM allowed to generate species richness maps by summing the probabilities of habitat suitability maps. This stacking method has been shown to overestimate species richness (Calabrese et al. 2014), but remains useful to compare different scenarios.
Species Traits and Range Shift Metrics
In order to assess potential correlations with range shift dynamics, three species traits (phylogeny, habitat, and temperature preference) were considered (see Supp Table S4a [online only] (nvac056_suppl_supplementary_material.docx)). Phylogeny was limited to suborder-level differences, Anisoptera and Zygoptera (i.e., dragonflies and damselflies). Larval habitat was either lentic (larval stage occurring in standing water), lotic (in flowing water), or generalist, following Hof et al. (2012). Temperature preference categorization was based on a European-scale Species Temperature Index (STI; Termaat et al. 2019), with cold- and warm-dwellers considered to have STI below 5.5°C and above 10°C, respectively. STI is defined as the annual average temperature a species experiences in Europe (excluding Russia). It is worth noting that temperature preference is not necessarily the result of an active choice by individual insects, but rather reflects the optimal range of temperatures to which a species is adapted.
The median latitude for each species' distribution was computed from the present-day ENM and future projections to estimate the expected speed of range shifts (assuming no dispersal limitation). We used median latitude (as in Rotenberry and Balasubramaniam 2020) rather than mean latitude as the median is less influenced by extreme values. For similar reasons, latitudinal amplitude was estimated with the inter-decile distance (the difference between the first and the last decile) to limit the influence of extreme values.
To characterize the nature of range shifts (i.e., identify expansion or contraction), from niche projections we estimated the net range change in square km (Scol – Sext; ‘col’ = colonization, ‘ext’ = extinction) between present day and the 2061–2080 period. We divided net range change by the present-day range area (Spres, in km2) to compute a range change index (following Buisson et al. 2010, see equation 1), exhibiting the magnitude of change accounting for more or less restricted species ranges. A positive value of this index means net range surface gain (range expansion) with a value of 100 meaning that the total range is expected to double. Conversely, negative values indicate the net range surface loss (range contraction) with a value of –50 meaning that the total range is expected to halve and –100 meaning total extinction of suitable conditions within the study area. Values near 0 indicate range stability or perfectly balanced range shift (surface gained ≈ surface lost).
Niche Over®lap Analyses
The probability of species interactions and potential changes in interaction due to range shifts can be captured by the niche overlap and co-occurrence analyses. From the binary maps produced by species ENMs, the suitable area or geographic space can be converted into a niche or environmental space (e-space) by sampling the environmental conditions of suitable versus unsuitable areas. Dealing with the e-space of different species allows for more relevant comparisons as e-space is relatively insensitive to spatial biases, even when distributions are out of equilibrium. We estimated the e-space for each species using the ‘humboldt’ R package (Brown and Carnaval 2019). The e-space was then used to determine niche overlap among species, based on both the present-day and future ENMs. We used Schoeners’ D as a measure for the niche overlap, which ranges from 0 (no overlap in environmental tolerance) to 1 (identical environmental preference, Warren et al. 2010). To limit overfitting, ENM locations were thinned so that locations were at least 50 km apart. This was performed only on a subset of 10,000 randomly selected ENM locations to decrease computation time. The accessible environment on the basis of which niche overlap calculations were performed was limited to a 100 km buffer area around each location. According to the definition given by Brown and Carnaval (2019), the niche overlap metric was calculated including non-analogous e-space, i.e., the accessible portion of e-space shared by both species in the pairwise comparison. E-space can be represented as a two-dimensional density grid for which the resolution was set to 50 instead of 100 to improve computation time and we kept the default kernel smoothing and kernel density threshold parameters of 1 and 0.001, respectively. E-space densities were also corrected by the abundance of different environments.
A multivariate approach was used to determine to what extent species traits would drive niche overlap patterns. We performed distance-based redundancy analysis (db-RDA; Legendre and Anderson 1999), using the capscale function from the ‘vegan’ R package, with three categorical variables (suborder, habitat, and temperature preference) as predictors. The contributions of these traits to niche overlap patterns were estimated using the varpart function in the same package. For each time period, we used a backward selection process, involving removal of non-significant traits (α = 0.05) until all the remaining traits were significant (Hyseni et al. 2021). Significance was determined using ANOVA, after adjusting for the main effects of the other predictors.
Co-occurrence Analyses
Species interaction probabilities were also assessed using the checkerboard analyses (Stone and Roberts 1990), performed with the ‘ecospat’ R package (Di Cola et al. 2017). We used the C-score, which has been shown to have good statistical power for detecting non-randomness (Gotelli 2000). The C-score quantifies the extent to which two species segregate (a C-score higher than expected by chance is indicative of a strongly antagonistic interaction such as competition or discrepancies in habitat requirements) or aggregate (a C-score lower than expected by chance is indicative of a positive interaction or similarities in habitat requirements). The magnitude of the significance is estimated with the standardized effect size (SES). In cases where there is no significant difference between observed and expected C-scores (i.e., a random pattern of co-occurrence; |SES| < 2), then no direct possible interaction is inferred. We simulated random co-occurrence patterns (i.e., null distributions of expected C-scores) using environmentally constrained matrices (i.e., habitat suitability from ENMs as site weights). We used habitat suitability scores derived from 10,000 randomly sampled cells of the present-day and future ENMs as the environmental constraints in order to rule out the effect of environmental conditions as much as possible (Peres-Neto et al. 2001), and used 10,000 permutations to generate matrices for the null models. The relationship between co-occurrence patterns and niche overlap for both the present-day and predicted future species range were explored, following Bar-Massada (2015). All the analyses were performed using R 3.6.2 software (R Core Team 2019).
Results
Quality of the Models
Of the 57 ensemble ENMs generated, 36 models (63%) had an AUC > 0.8 or Kappa > 0.6 ( Supp Table S2 [online only] (nvac056_suppl_supplementary_material.docx)), which is an indication of good-predictive ability (Thuiller et al. 2009). A negative correlation was found between AUC score and the latitudinal amplitude of ENM range (R2 = 0.70, Supp Fig. S2 [online only] (nvac056_suppl_supplementary_material.docx)), meaning that models for species with narrow distributions (occurring only in a limited fraction of the study area, e.g., only in the south) performed better than models for widespread species (occurring throughout Sweden). Despite this, we retained all the ENMs for subsequent analyses to avoid biased interpretations of patterns if we discarded widespread species by focusing only on ENMs with AUC > 0.8. The MOP analyses for future projections indicated reasonably low levels of extrapolation (values close to 1) in the northern portion of the study area for all the scenarios and GCMs ( Supp Fig. S3 [online only] (nvac056_suppl_supplementary_material.docx)). Thus, the present-day ENMs are fairly well transferable in such areas. However, the southern portion of Sweden is dominated by strict extrapolation, i.e., some of their climatic conditions exceed the range of the predictor values currently existing in Sweden. This suggests that range shifts occurring in the south should be interpreted with some caution.
To complement the reliability assessment of future projections, we projected the Sweden-trained ENMs over the present-day climate conditions in the central Europe and compared the suitability maps obtained with actual presence in this area ( Supp Table S5 [online only] (nvac056_suppl_supplementary_material.docx)). Moreover, the range change trends we predicted were also compared to the recent trends from other studies that used occupancy models. For most species, we found correct extrapolation to central Europe and consistent trends with other studies ( Supp Table S5 [online only] (nvac056_suppl_supplementary_material.docx)). Only Ophiogomphus cecilia showed an obvious discrepancy because of its odd distribution in the northern Sweden compared with occurrences at warmer temperatures in the rest of Europe.
Fig. 1.
Spatial distribution of present-day species richness and predicted change. (A) Results from stacked SDM showing present-day species richness in Sweden. (B–C) Projections of species richness change, in number of species, compared to present-day richness according to RCP4.5 (B) and RCP8.5 (C) climate change scenarios for 2061–2080. Species richness change is shown in shades of gray (in the online version gains are represented in shades of green and losses in shades of red), and for each shade the two numbers indicate the range of change. Only maps obtained with ACCESS1-0 GCM are shown, which are representative of the average GCMs (see Supp Fig. S5 [online only] (nvac056_suppl_supplementary_material.docx) for other GCMs).
Present-Day and Predicted Species Richness
Based on the stacked SDM output, we detected the highest species richness in the southern portion of the study area (ca. one third of Sweden; Fig. 1A). A low-reporting effort of odonate species in the north-western part of Sweden resulted in an underestimation of suitable habitat for widespread species. Despite this, we observed a strong negative latitudinal gradient in species richness, especially, along the Baltic coast, Mideast, and Northeast (Fig. 1A) where the reporting effort is greater than the further inland ( Supp Fig. S4 [online only] (nvac056_suppl_supplementary_material.docx)). ENM projections for the RCP4.5 and RCP8.5 trajectories suggest an increase in species richness in the northern two thirds of the country. This increase is more limited in high-altitude areas in the Northwest (Fig. 1B, Supp Fig. S5 [online only] (nvac056_suppl_supplementary_material.docx)). The level of inland colonization varied greatly according to the GCM used (see Supp Fig. S6 [online only] (nvac056_suppl_supplementary_material.docx) for some representative examples). The projections showed some species loss in areas of high present-day species richness, but this pattern needs to be interpreted with caution given model extrapolation results ( Supp Fig. S3 [online only] (nvac056_suppl_supplementary_material.docx)). Decreases in species richness predicted for some lowland and coastal areas of southern Sweden might in fact be compensated by the arrival of new species. For instance, four recent colonizers (originally excluded from niche modeling due to an insufficient number of records), plus four other species not yet present in Sweden, which are likely to reach the country within the next few decades. When included in analyses, these eight species compensated for decreases in richness projected for southern Sweden, except for the southernmost coastal areas, where these new arrivals would not compensate for local extinctions ( Supp Fig. S7 [online only] (nvac056_suppl_supplementary_material.docx)).
Winners and Losers
In total, 89% of the species modeled (51 out of 57, Fig. 2A) are expected to expand their range by 2061–2080 under the RCP8.5 scenario (88% under RCP4.5; Supp Fig. S8 [online only] (nvac056_suppl_supplementary_material.docx)). The average growth in suitable area is 133,260 ± 64,243 km2 (mean ± SD) and 78,789 ± 45,536 km2 according to RCP8.5 and RCP4.5 scenarios, respectively (Fig. 2B, Supp Fig. S8B [online only] (nvac056_suppl_supplementary_material.docx)). Suborder and habitat had no effect on the range change index (Kruskal–Wallis rank-sum test P = 0.55 for both). In contrast, temperature preference had a significant effect on range change (see blue and orange bars in Fig. 2A; Kruskal–Wallis rank-sum P < 0.001). A pairwise comparison using the Wilcoxon rank-sum test indicated that only the cold-dwelling species group (STI < 5.5°C) were significantly different from the intermediate (P < 0.001) and warm-dwelling species (P < 0.001) groups. The range change in the terms of absolute surface showed a similar pattern with regard to cold- and warm-dwelling species (Fig. 2B). Only a few species are expected to undergo genuine range shifts involving both contractions in the south and expansions in the north. The majority of these are cold-dwelling species such as Aeshna caerulea and Coenagrion johanssoni but also, interestingly, a few intermediate-temperature species such as Leucorrhinia caudalis and Epitheca bimaculata (Fig. 2B). These last two species are absent in the north, which suggests that it is not only species with ranges restricted to the north that show this pattern of predicted retreat in the south and expansion in the north. On the basis of the ensemble ENMs and projections, the median latitude of species range in Sweden is expected to increase from 58.9 ± 1.81 to 61.1 ± 2.06 degree north by 2061–2080 under RCP8.5 scenario. It represents a northward shift of 241 ± 132 km between 1996 and 2070 (average years of each period) thus an average median latitude shift of 3.25 ± 1.78 y–1 (Fig. 2C, Supp Table S4a [online only] (nvac056_suppl_supplementary_material.docx)). On average, shifts in the first and last latitude decile (assumed to reflect lower and upper range limits) are 1.11 ± 1.04 km y–1 and 5.67 ± 2.54 km y–1, respectively ( Supp Table S4b [online only] (nvac056_suppl_supplementary_material.docx)). Similarly, under RCP4.5 scenario a northward shift in median latitude of 135 ± 97 km is expected between 1996 and 2070, which represents an average shift of 1.82 ± 1.31 km y–1. Under this scenario, lower and upper range limits are predicted to shift at a speed of 0.67 ± 0.69 km y–1 and 3.71 ± 2.15 km y–1, respectively ( Supp Table S4b [online only] (nvac056_suppl_supplementary_material.docx)). Of the three species traits, only temperature preference had a significant effect on median latitude shift (the Kruskal–Wallis rank-sum test P < 0.001). In contrast, only warm-dwelling and intermediate species groups differed significantly (Wilcoxon rank-sum test P < 0.001), with an average median latitude shift of 2.39 ± 0.79 and 3.66 ± 0.88 km y–1, respectively. Median latitude shift of cold-dwelling species was more dispersed and reached on average 2.98 ± 2.32 km y–1.
Fig. 2.
Predicted range change between the present day and the 2061–2080 RCP8.5 scenario. (A) Range change index corresponds to net range change surface relative to the present-day range surface. Cold-dwelling (species temperature index, STI < 5.5°C) and warm-dwelling (STI > 10°C) species are shown in dark and light gray, respectively (blue and orange in the online version). Intermediate-temperature species (5.5°C<STI < 10°C) are represented by intermediate shades of gray. Species temperature preference categories are specified in Supp Table S4a [online only] (nvac056_suppl_supplementary_material.docx). (B) Range change surface. Surface area predicted to become unsuitable or suitable for species is represented as negative and positive values, respectively. (C) Average annual shift of range median latitude between the present day and 2061–2080. The values correspond to the average of seven different GCMs with standard deviation. Striped bars indicate species ensemble ENM with AUC and Kappa scores below 0.8 and 0.6, respectively. Range change for the 2061–2080 RCP4.5 scenario is given in Supp Fig. S8 [online only] (nvac056_suppl_supplementary_material.docx).
E-Space Niche Overlap Analyses
The level of niche overlap greatly varies among species pairs (Fig. 3). Five species tend to display very low-niche current overlap with other species, e.g., Ophiogomphus cecilia–Coenagrion johanssoni. These five species can be considered as the northern species in Sweden. Other species show high overlap, e.g., Erythromma najas and Somatochlora metallica. After the backward selection process of removal of non-significant traits, temperature preference and habitat were retained as significant in present-day and 2061–2080 RCP4.5 scenarios and only temperature preference in 2061–2080 RCP8.5 scenario (Table 1). Using db-RDA variance partitioning, we determined that temperature preference explained a larger portion of variance in niche overlap patterns than habitat and suborder for present-day and the two 2061–2080 scenarios considered ( Supp Fig. S9 [online only] (nvac056_suppl_supplementary_material.docx)). In accordance with the assumption of niche conservatism, the observed patterns seem to remain in future projections, but some differences emerged (Fig. 3). For example, Erythromma najas and Orthetrum coerulescens showed cases of higher overlap in the future than present. These differences may come from portions of non-overlapping e-space niches that may no longer exist in the future within the study area resulting in increased niche overlap or from some extrapolation of the model.
Fig. 3.
Pairwise niche overlap between the 57 species performed on the e-space derived from present-day (A) and future (B–C) ENMs. Niche overlap is estimated with Schoener's D index. Degree of niche overlap is indicated by shades of gray, with darker colors indicating high overlap (in the online version, the color scale is blue to red, with the latter indicating high overlap). Species traits (suborder, habitat, and temperature preference) are indicated by different shades on the left of the heatmaps. Future projections correspond to 2061–2080 RCP4.5 (B) and 2061–2080 RCP8.5 (C) scenarios based on a single GCM (CMCC-CM). Species are ordered according to a hierarchical clustering analysis performed using the complete-linkage method. Niche overlap computation parameters are detailed in the methods.
Table 1.
ANOVA table of db-RDA with species traits that have a significant effect on niche overlap patterns for present day and two 2061–2080 scenarios
Co-occurrence and Niche Overlap
While the e-space niche overlap analysis aims to determine species similarity in terms of habitat preference only, the co-occurrence analysis focuses on species interaction probability alone, without the habitat suitability component. A total of 94.8% of species pairs had co-occurrence patterns for both present day and the future (2061–2080 RCP8.5 scenario) significantly different from random (|SES| > 2). This high proportion of significant results might stem from the high number of ‘sites’ (10,000) considered. For present-day co-occurrence patterns, 90.2% of pairs were aggregating and 9.8% were segregating more strongly than expected at random (Fig. 4). The proportions remained similar for the future scenario (89.2% aggregating and 10.8% segregating).
On average, absolute values of future SES increased compared with the present day (mean ± SD, 17.9 ± 6.1 for present day, 31.1 ± 11.8 for future scenario) with 1,217 species pairs showing a decrease in SES (increased aggregation pattern) and 296 pairs showing an SES increase (decreased aggregation). Most changes occurred without a change in the sign of SES. However, 3% of pairs (n = 46) showed a change in the sign of SES between the present-day and future scenarios. For most of them (n = 37), the change was significant from negative to positive values of SES, meaning a change from aggregated to segregated co-occurrence patterns. As species interactions were not considered explicitly in stacked ENMs this change might reflect disjunct species distributions. This is supported by some markedly shifting species such as Coenagrion johanssoni and Leucorrhinia caudalis ( Supp Fig. S6 [online only] (nvac056_suppl_supplementary_material.docx)) involved in 16 and 10 of those segregating pairs, respectively.
For both present day and future, we found a negative correlation between co-occurrence (SES) and niche overlap (R2 = 0.34, Fig. 4) which is in line with the expected relationship between high-segregation pattern (more positive SES) and low-niche overlap. However, niche overlap and C-score SES did not show perfect linear correlation suggesting that some additional information was present in the e-space niche overlap compared with the purely geographic metric of co-occurrence.
Fig. 4.
Relationship between species co-occurrence and niche overlap for the present day and the 2061–2080 RCP8.5 scenario (dark circles and light squares, respectively, red and blue in the online version). Gray lines denote differences from present day to future for the same species pair. Only species pairs with significant standardized effect size (|SES|> 2) for both the present day and future are displayed (N=1,513 pairs). A positive/negative SES means that two species segregate/aggregate more strongly than what is expected by chance. The dashed lines represent the thresholds of significance. Niche overlap and co-occurrence for future scenario stem from one GCM only (CMCC-CM). Density plots show distribution spread and shift in both niche overlap and co-occurrence in the future scenario compared with the present day.
Discussion
Using an ENM framework, we inferred that northern European odonate assemblages will experience great compositional changes during the next few decades. Climate change will likely result in northward range shifts, northward range shifts of many species, thus changing species richness patterns and influencing species interactions. Furthermore, our results suggest that cold-dwelling species are more likely to suffer from range shifts as well as changes in species co-occurrence, which could ultimately reshape interaction networks.
Based on our ENM results, the majority of species in the study area is expected to undergo northward range expansion by 2080. Considering all the species, the average median latitude shift is expected to occur at a speed of 1.83 and 3.25 km y–1 under RCP4.5 and RCP8.5 scenarios by 2061–2080. Empirical northward expansion has already been documented for British odonates in a study comparing range expansion from 1960 to 1995 (Hickling et al. 2005). This study found an average range margin shift of 2.96 km y–1, which is close to what we found for the worst-case future climate change scenario. Interestingly, in Britain, a more recent analysis (Platts et al. 2019) confirmed range shifts at a similar rate (2 km y–1) for multiple invertebrate groups, with most of the variation within rather than among groups. The similarity in rates of range margin shifts across studies, including our study, suggests that climate change will push many insect taxa northward and also that the unlimited dispersal assumption in projections is reasonable. Odonates, in particular in Scandinavia, already have recent examples of remarkable northward range expansions with Anax imperator or Sympecma fusca that moved about ten-fold or more the average rate (Flenner and Sahlén 2008). Another recent study, covering the whole latitudinal range of Europe, but without a full coverage of the area, found that the majority of the 99 species studied had increased their range in the 1990–2015 period (Termaat et al. 2019).
For some cold-dwelling species, i.e., those with an STI < 5.5°C, our models predicted a range contraction since the habitat at their current southern range will become unsuitable. Interestingly, the study by Termaat et al. (2019) did not find that cold-dwelling species, as defined in their study, differed in range change compared with the warm-dwelling species between 1995 and 2015. However, that study defined cold-dwellers as species with STI below the mean STI of the European species (STI < 9.8°C) which corresponds to both our cold-dwelling and intermediate species groups. Our projections are more in line with temperature preference found to be positively correlated with the long-term trends in Germany, where cold-adapted species using standing waters showed the highest decreases (Bowler et al. 2021). It will be of interest to monitor if range contractions of cold-dwelling species predicted in our study will occur during the upcoming decades while we still are adapting conservation measures to limit local extinction and promote colonization of the most vulnerable species.
We predicted that species richness will increase in the northern areas and may decrease in the south due to local extinctions. Qualitatively, our models retrieved the species richness pattern known in Sweden based on observations with a rather high-species richness in the south-east and a much lower species richness in the northern areas (Kalkman et al. 2018). Yet, quantitatively our ENMs tend to underestimate the extent of some species ranges, especially in the northeast, thus, resulting in an underestimation of species richness. We are aware of two other studies that have used ENMs to study range shifts at the community level in the European odonates. Markovic et al. (2014) modeled aquatic biodiversity and found that odonate species richness were very little affected by climate changes until 2050 and explained this by their relatively high-aerial dispersal ability compared with strictly aquatic species groups. Li and Park (2020) modeled odonate distribution covering the western and the northern Europe and found an increase in species richness toward the north until 2050, and thereafter, a decline until 2080. The results from these two studies differ somewhat from ours, in that our models predicted clear range expansion and increased species richness toward the north, which was not found or was less evident in these other studies. The difference is probably due to our study being restricted to the northern part of Europe. We note that there was a slight decrease in species richness in some parts of the southern Sweden, which gives further support for the scenario of a south-north difference with regard to climate change. Nevertheless, this decrease is probably overestimated, because recent colonizers and other species that are yet to arrive and were thus not included in ENMs could counteract local losses in species richness. Hence, we conclude that the majority of species will experience range expansion in the north, which will lead to a less pronounced latitudinal gradient in Sweden. This also suggests that the current biogeographical and geological boundary known as Limes Norrlandicus (a boundary between the more oceanic and continental climatic regions), coinciding with the range limits of many species in southern Scandinavia, might be crossed by at least some odonate species during the next few decades. However, although we included a broad variety of ecologically relevant variables in our models, potential effects of high-latitude environments are almost impossible to anticipate in the absence of the experimental work (Spence and Tingley 2020). Since odonates are insects, we suggest that these change in richness patterns might hold for insects in general, because studies comparing different insect taxa show very similar responses to climate change (Maes et al. 2010, Domisch et al. 2011). In general, extinction risk of species will increase with future climate change (Urban 2015). But, at northern latitudes, only cold-dwelling species may be negatively affected, while warm-dwelling species will likely benefit from climate change.
We found interesting disparities in the niche overlap between species pairs. The general pattern we observed was that cold-dwelling species exhibited less niche overlap compared with the warm-dwelling species, a pattern which persisted until 2061–2080 for both climate change scenarios. Furthermore, the C-score analysis resulted in an overall increase in co-occurrence among species by 2061–2080, for both the climate change scenarios (Fig. 4). Such an increase in niche overlap was also found by van Beest et al. (2021) in a study on vertebrates at high northern latitudes (Greenland). More studies investigating changes in the niche overlap at the northern latitudes are needed to confirm that this pattern is generalizable in the northern latitudes and to understand its biological implications. These changes in the niche overlap were mostly seen in warm-dwelling species and although the majority of changes was toward higher overlap, there were notable exceptions. For instance, Erythromma najas, Somatochlora arctica, and Coenagrion puella showed a rather high present-day niche overlap with other species, followed by much lower niche overlap in the future. Conversely, Brachytron pratense, Sympecma fusca, Sympetrum sanguineum, and Leucorrhinia pectoralis experienced increased niche overlap in the future, compared with the present-day niche overlap with other species. More work is needed to understand why these different groups of species respond differently, how dependent on scale are these results, and whether this is not a result of the assumption of niche conservatism in niche modeling. Indeed, several studies have found high-niche conservatism in insects including odonates (Peterson et al. 1999, Wellenreuther et al. 2012, McCreadie et al. 2017, but see Hill et al. 2017, Torres et al. 2018). However, temperature-induced changes in interspecific competition have the potential to influence distributional patterns (Yang and Rudolf 2010, Suhling and Suhling 2013).
Besides intrinsic factors such as STI and competition (Yang and Rudolf 2010, Suhling and Suhling 2013), extrinsic factors might affect climate-induced range shifts and changes in species richness of odonates in the future. For instance, Li and Park (2020) found that, while STI was the main factor, extrinsic factors such as land cover type and water velocity also had a significant effect in predicting future distributional changes. However, the effect of intrinsic and extrinsic factors on niche overlap is not well studied. We included STI (temperature preference), suborder, and habitat in our analysis on niche overlap and found that only STI affected the niche overlap patterns significantly. The absence or limited effect of habitat as a trait driving niche overlap patterns could be explained by the coarse spatial resolution we used for computing niche overlap, which may have missed fine-scale niche differentiation that can occur within habitats (e.g. Hyseni and Garrick 2019).
The predictive abilities of our models were affected by a number of limitations and uncertainties. The level of thinning we applied to occurrence records served to limit the spatial autocorrelation in areas where the reporting effort was high, but for the species occurring in both the northern and southern Sweden the density of occurrence data included in the model was still higher in the south than in the north. This likely led to an underestimation of the present-day ranges and thus an overestimation of their predicted range expansion. Although AUC is assumed to be independent of prevalence (Manel et al. 2001), we found the poorest ENM predictions corresponded with the most widely distributed species. Moreover, a higher level of thinning would have degraded the performance of ENM without improving model fit ( Supp Table S6 [online only] (nvac056_suppl_supplementary_material.docx)). Some studies have managed to correct this bias by including sampling bias predictors in the models, such as distances to the main roads, cities, and protected areas (El-Gabbas and Dormann 2018). Also, the imperfect overlap between the date of occurrence records (2006–2020) and climatic variables (1979–2013) may lead to some small overestimations of range change given that range changes have certainly already occurred in the meantime. However, the magnitude of climate change experienced within the 1979–2020 period is likely to be minor compared to that between present-day and the 2061–2080 period. This makes our ENMs still relevant for projecting over this timescale. In addition, the use of several GCMs helps in obtaining reasonable projections. At last, our ENM projections involve three main assumptions. First, we assumed species distributions to be only constrained by abiotic conditions, which are assumed to be the dominant constraint at high latitudes (Sirén and Morelli 2019). Second, we assumed current distributions are not substantially affected by interspecific interactions at this scale (Eltonian noise hypothesis) and dispersal limitations. Third, ENM projections assume the availability of habitats in areas where colonization is predicted, which may be reasonable given that a substantial portion of Sweden—about 10%—consists of wetlands, especially in the northern areas (Jeglum et al. 2011). An interesting prospect would be to refine niche models using such wetland data.
Despite the great potential of citizen science to collect valuable biological data, the present study also underlines the disparities in reporting effort from easily accessible versus inaccessible areas (Mair and Ruete 2016). We like to note that enhanced data quality could be reached with a little extra effort to the citizen scientist. For instance, more detailed information with respect to the breeding status in odonate occurrence data (only available for 2% of the records) would allow to refine ENM so that areas suitable for breeding are inferred effectively (Patten et al. 2015). All this reaffirms the importance of reliable and extensive biodiversity data to help us evaluate range shifts in a changing world (Mihoub et al. 2017), and especially so at the northern latitudes.
Acknowledgments
We thank the Swedish Species Observation System and all the citizen scientists who submitted their sightings and made this study possible. The computations were enabled by resources in project SNIC 2021/22-171 provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
Author Contributions
M.P.: Formal analysis, validation, visualization, writing—original draft, writing—review and editing. F.J: Conceptualization, project administration, supervision, writing—review and editing. C.H.: Conceptualization, methodology, supervision, validation, writing—review and editing.
