Winners and losers over 35 years of dragonfly and damselfly distributional change in Germany

Recent studies suggest insect declines in parts of Europe; however, the generality of these trends across different taxa and regions remains unclear. Standardized data are not available to assess large‐scale, long‐term changes for most insect groups but opportunistic citizen science data are widespread for some. Here, we took advantage of citizen science data to investigate distributional changes of Odonata.


| INTRODUC TI ON
Recent studies suggest long-term declines of insect populations in different parts of Europe (Conrad et al., 2006;Hallmann et al., 2017;Homburg et al., 2019;Valtonen et al., 2017). Such declines are a major conservation concern, especially because they could have a broad range of cascading effects for other species (Cardoso et al., 2020;Hallmann et al., 2014). However, many studies on insect change are based on local or small-scale datasets and their representativeness of large-scale patterns is unclear. Assessing change over large spatial scales is difficult for most insect taxa due to a lack of standardized monitoring. Nonetheless, effective conservation policies strongly rely on knowledge of the large-scale trends of species. To facilitate conservation decision-making, there is an urgent need to make use of all available data to estimate large-scale and long-term changes in insect populations and communities.
While large-scale standardized insect monitoring is rare, at least beyond butterflies (van Swaay et al., 2008), large amounts of opportunistic data, without a common sampling protocol, are collected by citizen science (CS; Chandler et al., 2017). Citizen science data are associated with numerous statistical challenges but they have the advantage of a large spatial coverage and a reasonable time span.
Moreover, as some citizen scientists are active year-round, CS data also tend to capture a larger proportion of the biological community, including rare species, than standardized data (Bradter et al., 2018).
As CS data have become more accessible, there has been simultaneous development of methods to analyse such data (Isaac et al., 2014;Rapacciuolo et al., 2017;van Strien et al., ,2010van Strien et al., , , 2013. Occupancydetection models have emerged as one approach that is robust to different ways citizen scientists collect data, by explicitly modelling heterogeneity in survey effort and species detectability (Isaac et al., 2014).
Odonata are a good case study for the application of occupancydetection models to study long-term change because they are subject to extensive CS recording. Recent studies of dragonflies in different European countries suggest many species have increased in occurrence (Powney et al., 2015;Termaat et al., 2019). In general, there is evidence of recent population increases of many freshwater organisms in parts of Europe, thought to be due to better waste water treatment enabling recovery from previous impacts of water pollution (Van Klink et al., 2020;van Strien et al., 2016), especially following the implementation of the European Water Framework Directive (WFD) in the 2000s (Giger, 2009;Hering et al., 2010).
Climate changes signals have also been apparent by the northward expansion of southern species (Hickling et al., 2006;Ott, 2010).
Many studies on biodiversity change have focused only on the simple long-term mean trend of a species but this approach can mask a diversity of more complex temporal responses (Baranov et al., 2020;Outhwaite et al., 2020). Moreover, simple models can produce trend estimates driven by fluctuations in particular years quadrants. We also compiled data on species attributes that were hypothesized to affect species' sensitivity to different drivers and related them to the changes in species' distributions. We further developed a novel approach to cluster groups of species with similar patterns of distributional change to represent multispecies indicators.
Results: More species increased (45%) than decreased (29%) or remained stable (26%) in their distribution (i.e. number of occupied quadrants). Species showing increases were generally warm-adapted species and/or running water species, while species showing decreases were cold-adapted species using standing water habitats such as bogs. Time series clustering defined five main patterns of change-each associated with a specific combination of species attributes, and confirming the key roles of species' temperature and habitat preferences. Overall, our analysis predicted that mean quadrant-level species richness has increased over most of the time period.
Main conclusions: Trends in Odonata provide mixed news-improved water quality, coupled with positive impacts of climate change, could explain the positive trends of many species. At the same time, declining species point to conservation challenges associated with habitat loss and degradation. Our study demonstrates the great value of citizen science and the work of natural history societies for assessing large-scale distributional change.

K E Y W O R D S
biodiversity monitoring, citizen science, long-term change, occupancy-detection models, range-shifting, trait-based (Seibold et al., 2019). Recent analyses of invertebrate changes in the UK found both time periods of increases and time periods of decreases, which varied among taxa (Macgregor et al., 2019;Outhwaite et al., 2020). For instance, UK freshwater organisms decreased between 1970 and the mid-1990s but increased from then until 2010 (Outhwaite et al., 2020).
We used opportunistic data, including CS data, to study Odonata distributional change between 1980 and 2016 in Germany, which has the highest Odonata species richness in Europe (Brockhaus et al., 2015;Kalkman et al., 2018). We studied change in terms of the long-term trends but also the specific temporal patterns of change over time. To explore possible drivers of change, we identified species attributes (i.e. species-specific characteristics) associated with distributional changes.
This analysis was in part exploratory, to identify predictive attributes, and in part hypothesis-driven, for some key attributes assumed to link to specific drivers. We hypothesized that sensitivity to climate change may be affected by species' temperature preferences, while sensitivity to environmental management would be affected by species' habitat preference. Specifically, we predicted increases of warm-adapted and early spring species and decreases of cold-adapted species, as well as increases of running water species, because these are thought to most benefit from climate change and improved water and environmental management (Termaat et al., 2015). We hypothesized to find the strongest recovery of running water species especially during the 2000s when the WFD activities began, but recovery may also be seen earlier due to improved wastewater treatment in the 1990s. By contrast, increases associated with climate change were expected to be already visible from the 1980s. Finally, we tested whether communities had become more dominated by widespread and generalist species in line with biotic homogenization (Powney et al., 2015).

| Data compilation
In collaboration with the Society of German speaking Odonatologists (GdO) and various conservation agencies representing different federal states, we compiled 1,198,708 occurrence records on Odonata across Germany. The aggregated dataset comprised heterogeneous data, collected by both official and voluntary nature conservation organizations, without a common sampling protocol; however, experienced naturalists collected most of the data (Brockhaus et al., 2015;Mauersberger et al., 2013;Petzold & Fritzlar, 2014;Trockur, 2013).
Available data usually included information on observer or project, date of observation, life stage of species and geographic coordinates or ordinance survey quadrant (Meßtischblatt Quadrant, MTBQ) of c.
5 × 5 km (Goertzen & Suhling, 2019). We supplemented our dataset with observations from iNaturalist-removing duplicates and filtering to "research grade" (i.e. have been verified as reliable) to ensure the highest data quality, but these comprise only a small part of the total database (4.3% of the records, exclusion of these records does not affect our results).

| Data processing
We excluded: (a) data before 1980 because there were relatively few; (b) survey quadrants if they had not been visited in at least two separate years, so that they contain some information on change between years (Outhwaite et al., 2018), (c) observations of freshwater larvae because they must be sampled differently to adults and exuviae that are found on land (93% of the observations were of adults) and (d) species seen in less than 25% of years (4 species: Coenagrion hylas, Gomphus simillimus, Lestes macrostigma, Onychogomphus uncatus) due to insufficient data to estimate a trend. We set the bar low for inclusion so that we could examine the trends of rare species, but later we examined the uncertainty of the model estimates. The least observed included species was Oxygastra curtisii, with 43 records over 11 years while the most observed species was Ischnura elegans with 88,486 records over all 37 years. However, most species were seen in almost all years (median = 37, lower quartile = 36).
We included seasonal migratory species, such as Anax ephippiger and Sympetrum fonscolombii, even though they do not overwinter in Germany. Our trend estimation (described below) focuses on the trend among years during the main flight period of each species in Germany. Overall, 81 species of Odonata were included in the last published atlas for Germany (Boudot & Kalkman, 2015;Brockhaus et al., 2015)-our analysis included 77 species after applying the above exclusion criteria (Table S1). Following these filtering steps, our dataset was reduced to 1,073,129 records ( Figure 1).

| Species attribute data
We selected species attributes for which there are available data and are known to show large differences among species. In particular, we were interested in attributes that might affect species' vulnerability to climate change and habitat/land-use change.
Distribution: We estimated species' European geographic range size as the number of occupied grids (50 × 50 km) in a published atlas (Boudot & Kalkman, 2015).
Species temperature preference: Species' temperature preferences were calculated by overlaying each species' European distribution with an average temperature map from E-OBS v. 19e (Cornes et al., 2018) following other studies (Jiguet et al., 2007). For each species, we calculated the mean of the mean daily temperatures of occupied grid cells.
While we call this variable "temperature preference," its calculation did not aim to estimate species' optimal temperatures but rather to place species on a gradient from those preferring cooler temperature to those preferring warmer temperatures.
Life history: Data on voltinism, that is number of generations per year, were compiled from Corbet et al. (2006), complemented by expert knowledge within the co-author team. We applied a weighted mean of fuzzed-coded species affinities (values assigned to multiple categories reflecting the relative commonness of that category for the species, summing to 10 across all categories) to voltinism categories: multivoltine (coded as 5), bivoltinie (4), univoltine (3), semivoltine (2) and partivoltine (1). This weighted mean ranged between 1 (for a fully partivoltine species) and 5 (for a fully multivoltine species).
Phenology: Mean start dates of the main flight period were taken from Boudot and Kalkman (2015), which provided start date at a resolution of monthly tertiles. Species' phenologies vary geographically but the data presented were usually for Bavaria, southern Germany.
However, like for temperature preference, the aim was to create a variable that placed species on a gradient from those appearing early in the year to those appearing later.
Habitat: Main habitat preferences were classified according to descriptions in Dijkstra (2006) and Boudot and Kalkman (2015).
Each species was coded to whether they use the following habitats: streams, rivers, ponds, lakes, ditches, canals, fenland, bogs, forest and quarries/pits. Morphology: Hind wing length (median of lower and upper values) was taken from Dijkstra (2006).

Threat-level:
We compiled data on the 2015 red list classification for each species in Germany . We aligned the German categories with the international IUCN categories following Jansen et al. (2020).
Species attribute data are provided in Dryad datafile 1.

| Annual occupancy estimates
The finest unit of our analysis was a visit, referring to species observations collected on a given date in a given survey quadrant by a given observer/project (van Strien et al., 2010). Therefore, we organized the data into a list of species records seen on a given visit. Following others (Kéry et al., 2010), species' absences (non-detections) were inferred from observations of other species on a given visit. We used occupancy-detection models that account for imperfect detection to analyse the presence/absence (detection/non-detection) of species on a visit, which have been used in previous studies using opportunistic CS data (Outhwaite et al., 2020;van Strien et al., 2013) and tested in simulation studies (Isaac et al., 2014). By analysing the data at the level of a visit, the increase in the total annual number of records ( Figure 1) does not bias the trend estimation (Isaac et al., 2014).
Non-detections occur when a species inhabits a site but is not detected by the observer during a visit. Detection probability here also includes recording probability (i.e. the citizen scientist does not necessarily record all species that they detect). Detection probability is estimated by making an assumption about closure: a period during which species' occupancy does not change. We assume closure during the flight period of each species within the same year. The number of times a species was/was not seen during this period of closure informs on detectability. Hence, repeated visits to the same survey quadrants within the same year are required to estimate detection probabilities. In our dataset, on average, 51% of survey quadrants were visited at least twice in a given year (range over years between 49% and 58%). Moreover, revisited survey quadrants were on average visited 3 times within a year (range of mean over years between 2.3 and 3.9 times). We estimated the flight period of each species as between the lower 5% and upper 5% of days of year when each species was seen across all records. We only fit the model for each species to the subset of records during the flight period to meet the closure assumption. Models were fit to each species separately.
Letting z i,t refer to the true occupancy status for a species in a survey quadrat i in a year t, we modelled occurrence probability ( ) F I G U R E 1 (a) Spatial distribution of Odonata occurrence records for each c. 5 × 5 km survey quadrant in Germany. Shading refers to the number of occurrence records between 1980 and 2016 (subset to quadrats with records in at least 2 years). (b) Temporal distribution of the number of annual occurrence records (top) and annual number of visited survey quadrants (bottom). The red dashed vertical line denotes the start of our study period as a function of site and year variation. Year variation was modelled by including year as a fixed effect factor. Site variation was modelled as a series of random terms: ecoregion variation (Bundesamt für Naturschutz, 2008) at two spatial scales (level 1:7 coarser-scale regions and level 2: finer-scale 487 regions) and survey quadrant variation. While we were only interested in the mean yearly variation in occupancy, we included these spatial terms to account for the known large amount of spatial variation in species' occurrences across Germany. Hence, our occurrence model was: Detection probability (p) was modelled for each visit j to a given quadrant in a given year and assumed to depend on year and day of year (yday, accounting for species' phenology). Following Outhwaite et al. (2019), survey effort was modelled as a function of list length, that is number of species reported on a visit (a categorical variable with three levels: a single list (1 species, 31% visits), a short list (2-3 species, 24% visits) or a longer list (4 or more species, 45% visits, set as the reference level).
The observed detection data for a given species, y (0 or 1, for non-detection or detection), on each visit are then assumed to be drawn from a Bernoulli distribution conditional on the presence of the species in that quadrant and year: Using the estimated occurrences (z i,t ), we calculated a number of derived parameters: (a) the proportion of survey quadrants that were occupied by a species in each year (hereafter "occupancy proportion"); (b) the slope of a regression line through the annual occupancy proportions for each species (hereafter "trend"); and (c) ratio of occupancy proportion between the first and last year (Available in Dryad datafiles 1 and 2). These statistics were calculated during model fitting, and hence, the uncertainty in occupancy estimates was retained in each parameter.
The models were fit by Bayesian inference using JAGS, a program for fitting hierarchical models using Markov chain Monte Carlo simulation. We used vague priors for most parameters but a random walk prior, to share information across years, for the year fixed effect on occupancy (Outhwaite et al., 2018). We used 3 chains with 30,000 iterations, discarding the first 15,000 as burn-in. Model convergence was assessed by the Rhat statistic and traceplots. A small fraction of the annual predictions did not meet the standard threshold of <1.1 for Rhat (1.6% of annual occupancy estimates) but these Rhats were still reasonably close (all <1.26) and their exclusion had little effect on our overall results (see "Long-term trends"). Moreover, we carried our posterior predictive checks by calculating a Bayesian p-value. This was calculated as a Pearson chi-square statistic between the observed total number of species' detections in each year and the predicted number generated by simulating values from the fitted model. Bayesian p-values close to 0 or 1 would indicate poor model fit (Kéry & Royle, 2015). Our values were on average 0.46 across species (interquartile range = 0.43, 0.50). Predicted time series plots were also examined by Odonata experts in Germany to check for plausibility. The model code is provided in the Supporting information.
We ran a series of sensitivity analyses to check whether our results were robust to some modelling decisions. First, we ran dynamic-occupancy models that explicitly model changes in occupancy between years as either caused by persistence probability when the site was occupied in the previous year or by colonization probability when the site was not occupied in the previous year (Kéry et al., 2013). For most species, these models produced similar patterns to the simpler occupancy model above, but for the rarer species, the parameters of the dynamic-occupancy model showed lower convergence and large uncertainty, especially for the earlier years when there were less data. Hence, we used our simpler occupancy model that produced satisfactory results for all species. Second, we explored a more complex detection model in which the effect of day of year (linear and quadratic effects) was allowed to vary among years, to allow for possible changes in phenology, for example due to climate change. The species trend estimates of this model were highly correlated with the estimates from the simpler model (r = .99); hence, we proceeded with the simpler model. Last, we re-ran our original analysis with data limited to the set of quadrants that were surveyed at least once in each decade of our study period (e.g. at least once in the 80s, 90s, 00s and 10s), hence controlling for any effects of variation in site selection over time. This analysis used data from 13% of the original survey quadrants and again the estimated trends were highly correlated with those estimated using the full dataset (r = .97) but with larger uncertainties.

| Long-term trends
We analysed the interspecific variation in distribution trends using a linear regression model, with species attributes as explanatory variables and species' trend estimates as the response. Correlations among species' attributes were examined to check for any possible collinearity issues (all |r| < .6). We first conducted an exploratory analysis with each of the habitat variables in simple regression models. Using these models, we identified habitat variables that explained variation in trends (i.e. if p < .05). We then combined the selected habitat variables along with the other variables (temperature preference, voltinism, flight start date and wing length) in a multiple regression model of long-term trends. We performed step-wise deletion, removing insignificant attributes, to identify the best model.
We also used a linear model to compare the trends of species with different red list status. We re-ran our trait analyses excluding the species with a low number of detection records (Anax ephippiger, Boyeria irene and Oxygastra curtisii) and for which the Rhat values indicated the weakest convergence of annual occupancy and trend estimates, but the results did not change.
As species do not necessarily provide independent data points due to shared evolutionary histories, we checked whether the results of the regression were consistent after accounting for species relatedness. We used the taxonomic classification to build a simple phylogenetic tree with equal branch lengths for each taxonomic rank (Paradis & Schliep, 2019). We then included the tree as a correlation structure (corPagel-based on Brownian motion) in a generalized least-squares model (Paradis & Schliep, 2019). As this had little effect on the effect sizes of the attributes, and the estimated phylogenetic signal was close to zero (likelihood ratio test between models with and without the correlation structure, p =.36), we present the simpler model without this correlation structure.  Figure S4).

| Assemblage-level properties
To examine the implications of the species-level changes for the total Odonata species pool (referred to here as assemblage-level), we aggregated the predicted occupancy proportions of all species. We calculated for each year: mean species richness and diversity (Shannon index) using all species' occupancy proportions and community-weighted means for range size (mean range size of each species weighted by its occupancy proportion). We retained the uncertainty of the species' estimates by repeating the calculations for 1,000 random draws from the posterior distributions of the species' occupancies and calculating the upper and lower 2.5% quantiles.
We used R version 3.6.3 for all analysis. Statistical significance was assessed when 95% confidence intervals did not overlap zero.

| Long-term trends
Of the 77 species, we found that 35 were significantly increasing in (vulnerable, endangered, critically endangered) had more negative population trends compared with species of least concern or near threatened (p <.05, Figure 2c). Mean species detection probability was 0.25 (interquartile range across species = 0.20, 0.31) ( Figure S2).
In a multiple regression model, the most important attributes explaining variation in species' long-term trends were temperature preference, flight start date, wing length and river use (Figure 3).
Temperature preference was positively associated with species' trends: warm-adapted species increased while cool-adapted species decreased (Figure 4). Species that appear as flying adults earlier in the year and those with longer wings (in absolute terms) also had more positive trends. In the analysis of habitat associations, species occupying river habitats tended to increase, while species associated with bog habitats tended to decrease (Figure 4). While bog use was negatively associated with species' trends in the simple regression model, it was no longer significant in the multiple regression model ( Figure 3). This was probably because bog species were associated with colder temperature preferences (r = −.54), leading to a more uncertain effect of bog use after accounting for temperature preference. Voltinism had little importance (Figures 3 and 4). Together, temperature preference, wing length, flight start date and river use explained 27% of the variation in trends among species. . Sensitivity analysis defining one fewer cluster showed that clusters 4 and 5 might be grouped together as a declining species group ( Figure S4B). While defining one additional cluster suggested a group of seven species that primarily increased in the 1980s ( Figure S4C).

| Assemblage-level consequences
Predicted mean species richness per quadrat generally increased over the time period; however, some periods of decline were apparent dur- F I G U R E 2 Estimated nationwide trends in Odonata species' distributions. (a) long-term trends (y-axis, interpreted as mean annual change in % occupied sites, 0 = no change) for each species (x-axis), ordered by their magnitude of trend; (b) ratio in the estimated number of occupied quadrants between the last (2016) and first (1980) year (y-axis, 1 = no change, 2 = doubling) for each species (x-axes), ordered by their magnitude of change, and (c) boxplots (outliers are shown as solid points) of the association between long-term trend and red list classification (number of species in each group is shown in brackets). See Figure S1 for time series of individual species and Dryad datafile 1 for data on the trends of each species F I G U R E 3 Effect of species' attributes on their long-term distribution trends tested in a multiple regression model. Continuous variables (all except river and bog use) were scaled to units of 2 standard deviations to facilitate comparison with the binary habitat variables. The dashed red line is the line of no effect. Shown are the mean effect ±95% confidence interval. Table S2 provides the plotted values F I G U R E 4 Relationships between each attribute and species' distribution trends, each point is a species. The blue line is a fitted simple regression line. Boxplots are shown for river and bog use. The dashed horizontal line is the line of stable trends. Statistical significance at the 5% level is inferred when the 95% confidence intervals for the effects do not overlap zero-these are shown in Figure 3. The effects were significant for temperature preference, flight state date, wing length and river use (marginally) F I G U R E 5 Time series clusters and associated attributes. (a) Each cluster reflects a common pattern of change in occupancy over time for a group of species. The index represents the annual mean occupancy estimate relative to 1980. The number of species in each cluster was as follows: 35, 10, 9, 12 and 11. (b) The plots below shown boxplots or barcharts for attribute values within each cluster. See Table S3 for list of species in each cluster and Figure S3 for the species time series grouped by cluster. (a) was repeated by removing species with the most extreme occupancy changes in each cluster to check they were not driving the mean patterns-but similar patterns were found ( Figure S4)

| D ISCUSS I ON
Freshwater habitats have faced multiple anthropogenic threats, including eutrophication, acidification, climate change and canalization (Vörösmarty et al., 2010). Globally, freshwater vertebrate species are reported to be declining (He et al., 2019). Hence, our findings of many stable or increasing Odonata species since 1980 in Germany might seem surprising. However, our results are consistent with other studies on Odonata species that show increasing trends in Europe (Powney et al., 2015;Termaat et al., 2019;van Strien et al., 2016) and more generally positive trends in biomass of freshwater insects found in some datasets (Van Klink et al., 2020).
Hence, our results suggest a more complex range of species' distribution changes than the current simplistic narrative of declining insect populations. Additionally, our findings highlight the value of opportunistic data, especially from CS, in combination with statistical tools, for assessing large-scale biodiversity change.
Climate change probably plays a key role in the success of many Odonata species in Europe. As highly mobile organisms, many Odonata species may have adaptive capacity to respond to climate change, demonstrated by range-shifts reported in other countries (Flenner & Sahlen, 2008;Hickling et al., 2005) and may be responding stronger to climate change than many other terrestrial species (Hassall, 2015). We find that formerly rare warm-adapted species such as Crocothemis erythraea and Erythromma viridulum have undergone large range expansions across Germany. Increases in the occurrence of species typically appearing earlier in the year might also be linked to warming temperature. Earlier and longer reproduction seasons may increase the potential for more than one generation within a year (Braune et al., 2008) or allow early breeders to monopolize resources ahead of later-breeding species (Suhling & Suhling, 2013).
Additionally, we found that wing length was an important predictor.
One of the biggest winners has been the Emperor, Anax imperator, a strong flier with long wings (Rüppell, 1989).
Our findings may also reveal the impacts of improved environmental management, especially for rivers. Many running water species were increasing 1980 onwards even though the EU Water Framework Directive (WFD), which aimed to improve water quality, was not adopted until 2000 (Hering et al., 2010). This is probably because there was a range of other conservation and environmental management projects to improve water quality prior to 2000 in Germany (Detering, 2000;Giger, 2009). These projects included expansion of water purification plants; improved watercourse management (less removal of vegetation and disturbance of sediments); Fauna-Flora-Habitat Directive activities that targeted specific species and conservation measures to improve degraded wetlands.
Positive trends of dragonflies in the Netherlands were also thought to partly reflect habitat improvements (Termaat et al., 2015). River restoration projects in Europe have also allowed some recovery of other taxa, such as fish, though not necessarily to former historic states (Thomas et al., 2015). The success of riverine species may also reflect some synergism in the impacts of climate change and environmental management, because improvements in water quality may have facilitated climate change-driven range expansion by increasing the establishment success of immigrants (Braune et al., 2008).
Despite improvements in some freshwater habitats, we also identified a significant number of declining, "loser" species. Decreases of cold-adapted species could represent range contraction due to physiological stress under warming climates; however, more likely, Germany. The strong associations between species attributes and distribution trends support the use of trait-based vulnerability assessments for conservation decision-making (Conti et al., 2014).
However, increasing occurrences for some species represent expansion into new regions; while for other species, increases rather reflect recovery to formerly occupied parts of the range. The former group may be regarded as "neonatives" (Essl et al., 2019) and contribute to the development of novel species interactions and assemblages (Carrasco et al., 2018), with currently unclear repercussions for established/native species (Flenner & Sahlen, 2008;Suhling & Suhling, 2013). Decreasing distributions of other species have led to a decline in mean species richness during the 2010s. Standstills in the recovery of Odonata have been reported in the Netherlands (van Grunsven et al., 2020) and in the UK (Outhwaite et al., 2020).
Ongoing monitoring and synthesis of habitat assessments are needed to assess the likely cause.

| Data limitations
As our analysis is based on opportunistic data, concerns should always remain about the robustness of the trends. During our study period, there was a large increase in the reporting of species' observations. Simulation tests of occupancy models have shown them to be robust to different scenarios of survey effort change over time because analysis is performed on observations at the visit level (i.e. detections from the same date and site; Isaac et al., 2014); however, we cannot rule out that different changes in data collection across time and space affected some of our results. Also, our analysis only focused on changes in distributions and not changes in abundance.
For some species, increases or decreases in abundance may not yet have translated into changes in occupancy. We also only examined changes from 1980 onwards-most likely more historical data would highlight the earlier negative impacts of past water pollution and enable assessment of how much species have been able to recover to their former historical ranges (Goertzen & Suhling, 2019;Outhwaite et al., 2020).

| CON CLUS IONS
Using an extensive citizen science dataset, our analysis revealed a complex picture of positive and negative, linear and nonlinear distribution changes of Odonata species in Germany over the past 35 years. Climate change, habitat change and environmental management have probably all played a role. Cold-adapted habitat specialists of standing water habitats are likely to be most vulnerable to further environmental change, while increases of species associated with river habitats signal the conservation success that can be achieved by better environmental management. Overall, our study demonstrates the value of the intensive recording efforts of citizen scientists in the past, often promoted and coordinated by natural history societies, and highlights the need to support these efforts in the future, especially given signs of ongoing Odonata declines over the last decade.

DATA AVA I L A B I LT Y S TAT E M E N T
Species' annual occupancy estimates, trait data and distribution trend estimates are available in Dryad datafiles (https://doi. org/10.5061/dryad.73n5t b2wp).

ACK N OWLED G EM ENTS
This analysis was made possible through the efforts of many volunteer recorders over the past decades to whom we are very grateful. We thank especially the GdO-Gesellschaft

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13274.