Distribution trends of European dragonflies under climate change

Poleward range shifts of species are among the most obvious effects of climate change on biodiversity. As a consequence of these range shifts, species communities are predicted to become increasingly composed of warm‐dwelling species, but this has only been studied for a limited number of taxa, mainly birds, butterflies and plants. As species groups may vary considerably in their adaptation to climate change, it is desirable to expand these studies to other groups, from different ecosystems. Freshwater macroinvertebrates, such as dragonflies (Odonata), have been ranked among the species groups with highest priority. In this paper, we investigate how the occurrence of dragonflies in Europe has changed in recent decades, and if these changes are in parallel with climate change.


| INTRODUC TI ON
Climate change has a profound impact on the occurrence of many species of plants and animals (Parmesan & Yohe, 2003;Root et al., 2003;Walther et al., 2002). One of the most distinctive consequences is the poleward shift of species distribution ranges as a result of increasing temperatures, resulting in changes in the composition of species communities (Chen, Hill, Ohlemüller, Roy, & Thomas, 2011;Hickling, Roy, Hill, Fox, & Thomas, 2006;Kampichler, Van Turnhout, Devictor, & Van der Jeugd, 2012;Lindström, Green, Paulson, Smith, & Devictor, 2013;Mason et al., 2015). Species vary in their response to climate warming, due to different temperature requirements and different dispersal and colonization capacities. In general, warm-dwelling species and species with good dispersal capacity are more likely to be "winners" than cold-dwelling species and species with poor dispersal capacity (Franco et al., 2006;Rosset & Oertli, 2011;Virkkala & Lehikoinen, 2014). As a consequence, communities are predicted to become increasingly composed of warmdwelling, mobile species. This may seem straightforward, but the effects of climate change on species' trends and community compositions have only been studied for a limited number of taxa (but see Hickling et al., 2006;Mason et al., 2015), mainly birds, butterflies and plants (Bertrand et al., 2011;Britton, Beale, Towers, & Hewison, 2009;Clavero, Villero, & Brotons, 2011;Davey, Devictor, Jonzén, Lindström, & Smith, 2013;Devictor et al., 2012a;Jiguet et al., 2010;Roth, Plattner, & Amrhein, 2014;Virkkala & Lehikoinen, 2014). To gain a better understanding of how climate change affects total diversity, more taxa need to be covered, including taxa from different habitats. Freshwater macroinvertebrates should be ranked among the faunal groups with highest priority, as they have very different life histories from birds and butterflies and occupy very different ecosystems. They are known to react quickly to a wide range of changes in their habitats (Rosenberg & Resh, 1993). Furthermore, fresh water covers only 0.8% of the Earth's surface, while supporting almost 6% of all described species, most of which are insects (Dijkstra, Monaghan, & Pauls, 2014;Dudgeon et al., 2006). At the same time, they are among the most severely threatened ecosystems in the world, with aquatic species being more threatened than terrestrial species (Collen et al., 2014;Darwall et al., 2018;Dudgeon et al., 2006). For these reasons, freshwater invertebrates have been indicated as an essential future addition to Europe's biodiversity monitoring programme (Feest, 2013;Thomas, 2005).
Unfortunately, monitoring freshwater invertebrates comes with drawbacks. Most groups are so species-rich that collecting, sorting and identifying samples to species level require much effort and experience. Therefore, the number of specialists studying these groups is, in most countries, limited, which results in an incomplete picture of species' distributions. Dragonflies (Odonata) present an exception to this rule. Adult dragonflies are large, occurrence of dragonflies in Europe has changed in recent decades, and if these changes are in parallel with climate change.

Location: Europe.
Methods: We use data from 10 European geographical regions to calculate occupancy indices and trends for 99 (69%) of the European species. Next, we combine these regional indices to calculate European indices. To determine if changes in regional dragonfly communities in Europe reflect climatic warming, we calculate Species Temperature Indices (STI), Multi-species Indicators (MSI) and Community Temperature Indices (CTI).
Results: 55 of 99 considered species increased in occupancy at European level, 32 species remained stable, and none declined. Trends for 12 species are uncertain. MSI of cold-dwelling and warm-dwelling species differ in some of the regions, but increased at a similar rate at European level. CTI increased in all regions, except Cyprus.
The European CTI increased slightly.
Main conclusions: European dragonflies, in general, have expanded their distribution in response to climate change, even though their CTI lags behind the increase in temperature. Furthermore, dragonflies proved to be a suitable species group for monitoring changes in communities, both at regional and continental level.

K E Y W O R D S
citizen science data, climate change, Community Temperature Index, Multi-species Indicator, Odonata, Species Temperature Index colourful insects, which are easy to spot and relatively easy to identify at species level, making them attractive to a large public. With a manageable 143 species recorded in Europe (Kalkman et al., 2018), they constitute a suitable group for citizen science projects. Furthermore, dragonflies are well established as useful organisms to assess and monitor aquatic and wetland ecosystem quality (Oertli, 2008), and they are known to react quickly to climate change (Bush, Theischinger, Nipperess, Turak, & Hughes, 2013;Hassall, 2015).
In most European countries, dragonfly recording has increased in recent decades, mediated by the publication of several good field guides and national distribution atlases. This has resulted in a steep increase in available distribution data from citizen science projects and the publication of a European distribution atlas in 2015 (Boudot & Kalkman, 2015). The majority of these distribution data refer to records collected without standardization, which are unsuitable for straightforward calculation of distributions trends. However, previous studies have shown that these "opportunistic" records can be used to derive reliable trend estimates of dragonflies on a national scale, if occupancy models are applied. These models take the imperfect detection of species into account, and thereby, they may simultaneously correct for observation and reporting bias as well (Isaac, Van Strien, August, De Zeeuw, & Roy, 2014;Van Strien, Termaat, Groenendijk, Mensing, & Kéry, 2010;Van Strien, Van Swaay, & Termaat, 2013). Moreover, Van Strien,  showed in a pilot study, using records of a single species from five western European regions, that occupancy indices from multiple regions can be combined to calculate supraregional indices and trends.
In this paper, we investigate how the occurrence of dragonflies in Europe has changed in recent decades, and if these changes are in parallel with climate change. We use distribution data from 10 European geographical regions-ranging from Sweden to Cyprusto calculate occupancy indices and trends for as many dragonfly species as possible. Next, we combine these regional indices to calculate European indices. To determine if changes in regional dragonfly communities in Europe reflect climatic warming, we calculate Species Temperature Indices (STI), Multi-species Indicators (MSI) and Community Temperature Indices (CTI). We hypothesize that (a) warm-dwelling species have more positive trends than cold-dwelling species, that, as a consequence, (b) warm-dwelling species have increased their share in regional communities and (c) that these effects increase on a south-north gradient through Europe, as the ratio of warm-and cold-dwelling species decreases with increasing latitude.
F I G U R E 1 Participating European geographical regions, here considered as countries or lower administrative divisions 2 | ME THODS

| Species records
We gathered dragonfly distribution records, from 1990 onwards, from the following European geographical regions (countries or lower administrative divisions, hereafter referred to as "regions"): Sweden, Britain (United Kingdom excluding Northern Ireland), the Netherlands, North Rhine-Westphalia (a German state), Bavaria (a German state), Flanders (a Belgian region, including Brussels region), Wallonia (a Belgian region), France, Andalusia (a Spanish autonomous community) and Cyprus (Figure 1). The resulting data set included records of 99 species, which equals 69% of all dragonfly species recorded in Europe (Kalkman et al., 2018).
All records used in this study cover adult dragonflies only. The majority of these records are "opportunistic," that is, collected without a standardized field protocol and without a design ensuring the geographical representativeness of sampled sites. The period of data coverage and the number of records per unit area vary considerably among regions, depending on data availability (Supporting Information Table S1). All data in each region were validated by experts to prevent false-positive records. To standardize the geographical reference system, all observations were mapped in the ETRS89/ETRS-LAEA (EPSG:3035) reference system. Because we used 1 × 1 km grid squares as the definition of a site in our analyses, all observations were referenced to 1 × 1 km ETRS-LAEA squares.

| Generating non-detection data
Occupancy models require detection/non-detection data col- Almost all data obtained were records of species presence. The non-detection records of a given species were generated from the information of sightings of other dragonfly species, following Van Strien et al. (2010) and Van Strien, Van Swaay et al. (2013). Any observation of a given species was taken as 1 (detection), whereas we rated 0 (non-detection) if any species other than the given species had been reported by an observer at a particular 1 × 1 km site and on a particular date within the species' closure period.

| Annual occupancy estimates and trends: regional level
First, we calculated annual occupancy per species, for each region separately. We applied the same dynamic occupancy model as Van Strien et al. (2010) and Van Strien, Van Swaay et al. (2013) to estimate annual occupancy ψ, adjusted for detection probability p.
Because all parameters in the model may differ between regions, the analyses were performed separately for each region and the regional results were combined in a second step. The description of the model is derived from Royle and Kéry (2007) and Royle and Dorazio (2008). Here, ψ is the proportion of suitable 1 × 1 km squares that is occupied. A square is defined as suitable if the species had been recorded there at least once in 1990-2008. The occupancy model consists of two hierarchically coupled submodels, one for occupancy and one for detection, the latter being conditional on the occupancy submodel. The occupancy submodel estimates annual probability of persistence φ t and of colonization γ t and computes the annual occupancy probability per site recursively through: Thus, whether site i occupied in year t−1 is still occupied in year t is determined by the persistence probability, and whether site i unoccupied in year t−1 is occupied in year t depends on the colonization probability. All occupancy probabilities per site together yield the estimated annual number of occupied 1 × 1 km sites per region.
The same sites were included in the analysis for all years; estimates for sites not surveyed during some years were derived from sites that were surveyed in those years.
The detection submodel estimates the yearly detection p, but in addition, p is made as a function of covariates. We used the Julian date as a covariate for p because the detection of the species is expected to vary over the season, due to changing population size during the course of the flight period. Detection is also reduced if observers do not report all their sightings. Hence, we include the incompleteness of recording as a covariate for detection. We distinguished: (a) single records of any species on one site and date without records of other species, (b) short day-lists, that is, records of two or three species made by a single observer on one site and date and (c) comprehensive day-lists, that is, records of more than three species per observer, site and date. These lists may or may not include the species in question. These category thresholds are sufficiently low not to be confounded by real differences in species number between sites. In most 1 × 1 km sites in the regions, there are more than three species to be found and often many more.
Effects of both covariates were included in the detection submodel via a logit link: where p ijt is the probability to detect the species at site i during visit j in year t, α t is the annual intercept implemented as a random effect, β 1 and β 2 are the linear and quadratic effects of the date of visit j and δ 1 and δ 2 are the effects of short day-lists and comprehensive day-lists, relative to single records.
We fitted the models in a Bayesian mode of inference using JAGS (Plummer, 2017) on the computer cluster LISA (https://subtrac.sara.nl). We chose uninformative priors for all parameters, using uniform distributions with values between 0 and 1 for all parameters except δ 1 and δ 2 (values between −10 and 10), β 1 , β 2 (values between −10 and 10) and α t (values between 0 and 5) for the standard deviation of the normal distribution used as prior for the random year effect.
For each analysis, we ran three Markov chains with sufficient iterations to ensure convergence as judged from the Gelman-Rubin Rhat statistic and saved the last 93 iterations for use at supraregional level. This number of iterations is an empirically obtained compromise between the reliability of the estimates and data handling capacity. The model produced annual estimates of occupancy per region, which were converted into annual indices with first year = 100.
The trend in occupancy was considered significant if its confidence interval did not include zero.

| Annual occupancy estimates and trends: European level
In the next step, the annual occupancy estimates per region were aggregated to obtain European occupancy indices and trends for the period 1990-2015. Missing yearly values from a particular region were estimated ("imputed") from averaged year-to-year occupancy ratios in all other regions. For example, 1990 was missing in the Swedish data set. To impute occupancy estimates of Swedish species, we applied the 1991/1990 ratios from all other regions with data from both years. As a consequence of these imputations, confidence intervals increased for years with lacking data from one or more regions, especially when these were large regions (e.g., France).
Regions differ in the number of sites surveyed, so a naive aggregation has the risk of biased European trends. Hence, we developed a procedure to weigh regions according to the sampling intensity in relation to the range of species in each region. This procedure is an adaptation of procedures applied by Van Swaay, Plate, and Van Strien (2002) and Gregory et al. (2005). Weights were calculated as the quotient of relative range and relative sampling intensity to compensate for oversampling and undersampling. Relative range was defined as the range of a species in a region, as a percentage of its total range in all regions for which an occupancy index could be obtained. Relative sampling intensity was defined as the number of 1 × 1 km squares surveyed at least once in this period within the regional range of the

| Multi-species Indicators
To determine whether warm-dwelling species have more positive trends than cold-dwelling species, we calculated Multi-species Indicators (MSI), by combining the trends in occupancy indices of cold-dwelling and of warm-dwelling species, respectively. We did this for each region separately and for Europe as a whole. Colddwelling species were defined as species with STI lower than 9.8°C, which is the median STI of all species included in our study.

| Occupancy indices and trends
The number of species for which a regional trend could be calculated with sufficiently low standard errors, that is, standard errors low enough to detect a 5% or higher annual increase or decline, ranged from five for Cyprus to 79 for France (Table 1). In total, we were able to calculate trends with sufficiently low standard errors for 90 of 99 species in our data set, for at least one of the regions (Supporting Information Data S1).
In 7 out of 10 regions, more species increased than decreased their occupied range (

| Multi-species Indicators
MSI of warm-dwelling species were increasing in all regions ( Figure 3).

Surprisingly, MSI of cold-dwelling species also increased in Sweden,
Britain, the Netherlands and North Rhine-Westphalia. In Flanders, F I G U R E 3 Multi-species Indicators (MSI) of warm-dwelling species (Species Temperature Index > 9.8°C) and cold-dwelling species (Species Temperature Index < 9.8°C) per geographical region (from north to south) and for Europe. The first year with data was set to 100.
Smoothed trend lines were plotted through the year effects to summarize overall change. Shaded areas represent confidence intervals. Please note that y-axes differ Wallonia and France, MSI of cold-dwelling species was stable, in Bavaria it declined, and in Andalusia, it was uncertain (Table 2). Cyprus has only one cold-dwelling species (Enallagma cyathigerum), which increased.
Comparing MSI trends of warm-dwelling and cold-dwelling species (Table 2) shows the former was significantly more positive in Britain, the Netherlands, Flanders, Wallonia, Bavaria and France.
At a European level, however, the MSI of warm-dwelling and colddwelling species both increased at a similar rate.

| Community Temperature Indices
CTI increased in all regions, except Cyprus (Table 3; Supporting Information Figure S1). The most significantly increasing CTI was

| D ISCUSS I ON
We indicate that climate change has changed dragonfly occurrence at the community level as well.

| Regional level
We hypothesized that (a) warm-dwelling species show more positive trends than cold-dwelling species and this was confirmed on average, up to 9.5 × 10 −3°C year −1 for the Netherlands (Table 3), this "warming" of regional communities evolves more slowly than the increase in temperature itself (1.1 × 10 −2°C year −1 , after correcting for the difference in latitudinal gradient between CTI and actual temperature; Devictor et al., 2012a), but the difference for the Netherlands is minimal. Thus, dragonflies in Europe are accumulating a substantial "climatic debt," that is, the difference between shifts in temperature and shifts in distribution (Devictor et al., 2012a;Menéndez et al., 2006), which varies between regions.
Ultimately, our hypothesis that (c) trends in regional CTI increase on a south-north gradient through Europe is rejected. Highest CTI increases were found for regions on a moderate latitude (the Netherlands, Flanders, Wallonia) and for Andalusia (although measured over a shorter time span), while lowest CTI increases were found for Britain, France and Bavaria. Regions differ in size and subsequently in latitudinal gradient. This may, in theory, affect regional occupancy trends (and thus regional CTI trends) to some extent, possibly limiting the validity of a comparison at the regional level. Calculating CTI across equally sized latitudinal bands would be a preferable approach, but requires a higher data density in some of our regions than is currently available. The MSI trends of cold-dwelling and warm-dwelling species (Table 2) do not show a structural difference between larger and smaller regions indicating that it is unlikely that CTI trend differences are confounded by differences in region size.

TA B L E 3 Slope of Community
Temperature index (CTI) per geographical region (from north to south) and for Europe TA B L E 2 Multi-species Indicator (MSI) trends of cold-dwelling species (Species Temperature Index < 9.8°C) and warm-dwelling species (Species Temperature Index > 9.8°C) per geographical region (from north to south) and for Europe Note. SE: standard error; p: probability value. Moderate increase = significant increase ≤ 5% (p < 0.05); stable = no significant change; moderate decline = significant decline ≤ 5% (p < 0.05); uncertain = no significant change and standard errors too large to detect a 5% trend if it had occurred.
shows that warm-dwelling species have slightly increased their share. To compare this outcome with the trends in CTI of European birds and butterflies (provided by Devictor et al., 2012a, as based on presence-absence data), we re-calculated the European CTI of dragonflies for the same period 1990-2008. Over these 18 years, CTI of European dragonflies increased with 2.4 × 10 −3°C year −1 , which is comparable with the increase in CTI of European butterflies (2.5 × 10 −3°C year −1 ) and considerably greater than the increase in CTI of European birds (1.9 × 10 −3°C year −1 ). This is in line with the well-known ability of dragonflies to quickly colonize new habitats (Corbet, 1999). Dragonflies should probably be considered as more dispersive than butterflies, which, for their part, may show a quicker community response at a local scale, due to their generally shorter life cycle. The net outcome of these opposing differences may have resulted in a similar CTI trend between dragonflies and butterflies.
The slower response of bird communities to climatic warming has been suggested by Devictor et al. (2012a) to be a consequence of their slower population turnover.
In conclusion, climate change has a considerable positive impact on the occurrence of dragonflies in several European regions.
However, at a continental scale, CTI's are changing only slowly so far, due to the relatively positive response of cold-dwelling species.

| Limitations of CTI
Several authors have highlighted the CTI as a useful tool for assessing the effect of climate change on the composition of communities (Devictor, Julliard, Couvet, & Jiguet, 2008;Lindström et al., 2013;Roth et al., 2014). However, our results show that a stable CTI does not necessarily mean that climate change is not affecting the occurrence of species. In Sweden, many dragonfly species have benefited from climate warming, including species of cool conditions. This has led to increasing MSI trends for both warm-dwelling and colddwelling species, while leaving CTI almost unaffected. We therefore recommend a reviewing of CTI in relation to MSI of warm-dwelling and cold-dwelling species, especially in high-latitude regions where temperatures may have limited species with low STI as well as high STI. In addition, we know that many dragonfly species have substantially expanded their range northwards (Boudot & Kalkman, 2015;Hickling, Roy, Hill, & Thomas, 2005;Ott, 2010). Dragonfly communities have changed as a result of these expansions, yet this is masked by an increase in other species resulting in a quite stable CTI. For example, it is likely that the reduction in organic pollution and nutrient input in the last quarter of the 20th century has compensated the effects of increasing temperature for species that are sensitive to low oxygen levels (Ketelaar, 2010;Termaat, Van Grunsven, Plate, & Van Strien, 2015). These limitations of CTI as an indicator of climate change are also relevant when calculations are based on local abundances instead of regional distributions, even though CTI trends based on abundances show a stronger response to climate change than when based on occupancy Virkkala & Lehikoinen, 2014).

| Threats of climate change
We were able to calculate European trends in occupancy for 87 species (88% of species occurring in our data set). Fifty-five of these species have increased from 1990 to 2015, while 32 have remained stable and none have declined. This is a remarkably positive outcome, given the fact that the conservation status of many freshwater organisms is known to have deteriorated globally (Collen et al., 2014;Dudgeon et al., 2006). Although we recognize that some species with a stable trend in occupancy ( These studies may give an estimation of distances covered by the majority of the studied population, but undoubtedly miss dispersal events by individuals over much longer distances (see also Suhling, Martens, & Suhling, 2017), leading to a severe underestimation of maximum dispersal abilities. These extreme dispersal events may be rare and seldom noticed, but they determine the pace at which species distributions may expand. Next to annual estimates of occupancy, occupancy models also provide annual estimates of persistence and colonization. These parameters may be more informative for future research on the effect of variation in species' dispersal abilities.

| Trends from distribution data
We based our trend calculations on readily available distribution data from ten European regions, using occupancy models to account for imperfect detection. These records allowed us to assess occupancy indices without the need for a standardized fieldwork programme.
As such, our method immediately informs about distribution trends and may serve as an "early warning system" for species with a deteriorating conservation status and, by proxy, the quality of freshwater habitats (Oertli, 2008). However, our study lacks data from eastern Europe and we have insufficient data from southern Europe adequately to represent that area. Unfortunately, 18 out of 19 dragonfly species considered to be threatened at European level are confined to southern or eastern Europe (Kalkman et al., 2018). Considering this, our European indices and trends may be biased to some degree at pan-European level. However, dragonfly data sets are rapidly growing in many countries, including several eastern and southern European countries (Boudot & Kalkman, 2015). Moreover, a network of European odonatologists has expanded over the past few years and the usefulness of a European dragonfly monitoring scheme is gaining attention. We are therefore confident that European indices will become more robust with future updates and will have a better geographical coverage.

| Future prospects
Overall, this study has shown that dragonflies present a suitable species group to gain better understanding of biodiversity changes and their causes, including climate change, and that suitable data needed for these analyses are becoming available. Dragonflies may therefore satisfy the need for a biodiversity indicator based on freshwater invertebrates (Feest, 2013). They are likely to represent other taxa which are primarily warm-adapted. Using opportunistic data analysed with occupancy models enables the assessment of species' distribution trends on both regional and European scale.
These trends inform about the state of freshwater habitats, which is urgently required (Darwall et al., 2018). Hence, we suggest adding dragonflies as an indicator group to the European biodiversity monitoring programme (European Environmental Agency, 2012), to invest in the extension of a European dragonfly recording network, and to encourage the centralization of European dragonfly distribution data.

ACK N OWLED G EM ENTS
We thank all dragonfly observers who contributed to this study by making their observations available to their regional data managing organizations. Data from Sweden were obtained from the