Different roles of concurring climate and regional land-use changes in past 40 years’ insect trends

All statistical analyses were performed through R version 4.1.0⁵⁰. Besides the explicitly mentioned packages, the R packages cowplot⁵¹, data.table⁵², dplyr⁵³, ggplot2⁵⁴, itsadug⁵⁵, purrr⁵⁶, raster⁵⁷, sf⁵⁸, sfheaders⁵⁹, tibble⁶⁰ and tidyr⁶¹ were key for data handling, data analysis, and plotting. Posterior distributions were summarised through means and credible intervals (CIs). CIs are the highest density intervals, calculated through the package bayestestR⁶². To summarise multiple posterior distributions, 5000 Monte Carlo simulations were used.

Study region

The study included data from the whole of Switzerland. As an observation unit for records, we chose 1 × 1 km squares (henceforth squares). Switzerland covers 41,285 km², spanning a large gradient in elevation, climate and land use. It can be divided into five coarse biogeographic regions based on floristic and faunistic distributions and on institutional borders of municipalities⁶³ (Fig. 1b). The Jura is a mountainous but agricultural landscape in the northwest (~4200 km², 300–1600 m asl; annual mean temperature: ~9.4 °C, annual precipitation: ~1100 mm (gridded climate data here and in the following from MeteoSwiss (https://www.meteoswiss.admin.ch), average 1980–2020, at sites ~500 m asl.)). The Jura is separated from the Alps by the Plateau, which is the lowland region spanning from the southwest to the northeast (~11,300 km², 250–1400 m asl, mostly below 1000 m asl; ~9.5 °C, ~1100 mm). It is the most densely populated region with most intensive agricultural use. For the Alps, three regions can be distinguished. The Northern Alps (~10,700 km², 350–4000 m asl; ~9.2 °C, ~1400 mm), which entail the area from the lower Prealps, which are rather densely populated and largely used agriculturally, up to the northern alpine mountain range. The Central Alps (~11,300 km², 450–4600 m asl; ~9.5 °C, ~800 mm) comprise of the highest mountain ranges in Switzerland and the inner alpine valleys characterised by more continental climate (i.e., lower precipitation). Intensive agricultural land use is concentrated in the lower elevations and agriculture in higher elevations is mostly restricted to grassland areas used for summer grazing. The Southern Alps (~3800 km², 200–3800 m asl; ~10.4 °C, 1700 mm) range from the southern alpine mountain range down to the lowest elevations of Switzerland and are clearly distinguished from the other regions climatically, as they are influenced by Mediterranean climate, resulting in, e.g., milder winters. Besides differences between biogeographic regions, climate, land use and changes therein vary greatly between different elevations (Supplementary Fig. S9). To account for these differences, we split the regions in two elevation classes at the level of squares. All squares with a mean elevation of less than 1000 m asl were assigned to the low elevation, whereas squares above 1000 m asl were assigned to the high elevation (no squares in the Plateau fell in the high elevation). This resulted in nine bioclimatic zones (Fig. 1b), for which separate species trends were estimated in the subsequent analyses. The threshold of 1000 m asl enabled a meaningful distinction based on the studied drivers (climate and land-use change) and was also determined by the availability of records data (high coverage in all nine bioclimatic zones).

Species detection data

We extracted records of butterflies (refers here to Papilionoidea as well as Zygaenidae moths), grasshoppers (refers here to all Orthoptera) and dragonflies (refers here to all Odonata) from the database curated by info fauna (The Swiss Faunistic Records Centre; metadata available from the GBIF database at https://doi.org/10.15468/atyl1j, https://doi.org/10.15468/bcthst, https://doi.org/10.15468/fcxtjg). This database unites faunistic records made in Switzerland from various sources including both records by private persons and from projects such as research projects, Red-List inventories or checks of revitalisation measures. Only records with a sufficient precision, both temporally (day of recording) and spatially (place of recording known to the precision of 1 km² or less), were used for analyses. Besides temporal and spatial information, information on the observer and the project (if any) was obtained for each record. All records made by a person/project on a day in a square were attributed to one visit, which was later used as replication unit to model the observation process (see below).

We included records from the focal time range 1980–2020. Additionally, we included records from 1970–1979 for butterflies in occupancy-detection models to increase the robustness of mean occupancy estimates. We excluded the mean occupancy estimates for these additional years from further analyses to cover the same period for all groups. Prior to analyses, following the approach in ref. 26, we excluded observations of non-adult stages and observations from squares that only were visited in 1 year of the studied period, because these would not contain any information on change between years⁶⁴. This resulted in 18,018 squares (15,248 for butterflies, 9870 for grasshoppers, 5188 for dragonflies) and 1,448,134 records (879,207 butterflies, 272,863 grasshoppers, 296,064 dragonflies) that we included in the analyses (Supplementary Fig. S2). The three datasets for the different groups were treated separately for occupancy-detection modelling, following the same procedures for all three groups. To determine detections and non-detections for each species and visit, which could then be used for occupancy-detection modelling, we only included visits that (a) did not originate from a project, which had a restricted taxonomic focus not including the focal species, (b) were not below the 5% quantile or above the 95% quantile of the day of the year at which the focal species has been recorded²⁶ and (c) were from a bioclimatic zone, from which the focal species was recorded at least once.

Occupancy-detection models

We used occupancy-detection models^65,66 to estimate annual mean occupancy of squares for the whole of Switzerland and for the nine bioclimatic zones for each species (i.e., mean number of squares occupied by a species), mostly following the approach in ref. 26. We fitted a separate model for each species, based on different datasets for the three groups. We included only species that were recorded in any square in at least 25% of all analysed years. Occupancy-detection models are hierarchical models in which two interconnected processes are modelled jointly, one of which describes occurrence probability (ecological process; used to infer mean occupancy), whereas the other describes detection probability (observation process)⁶⁵. The two processes are modelled through logistic regression models. The occupancy model estimates occurrence probability for all square and year combinations, whereas the observation model estimates the probability that a species has been detected by an observer during a visit. More formally, each square i in the year t has the latent occupancy status z_i,t, which may be either 1 (present) or 0 (absent). z_i,t depends on the occurrence probability ψ_i,t as follows

The occupancy status is linked to the detection/non-detection data y_i,t,j at square i in year t at visit j as

where p_i,t,j is the detection probability.

The regression model for occurrence probability (occupancy model) looked as follows

with μ_o being the global intercept, elevation_i being the scaled elevation above sea level and α_o1,i, α_o2,i and γ_r(i),t being the random effects for fine biogeographic region (12 levels, Supplementary Fig. S10; these were again defined based on floristic and faunistic distributions and followed institutional borders⁶³), square and year. The random effects for fine biogeographic region and square were modelled as follows:

and

The random effect of the year was implemented with separate random walks per zone following ref. 67, which allowed the effect to vary between the nine bioclimatic zones, while accounting for dependencies among consecutive years. Conceptually, in random walks, the effect of 1 year is dependent on the previous year’s effect, resulting in trajectories with less sudden changes between consecutive years. This was implemented as follows:

with

The regression model for detection probability (observation model) looked as follows

where μ_d is the global intercept, yday_j is the scaled day of the year of visit j, shortlist_j and longlist_j are dummies of a three-level factor denoting the number of species recorded during the visit (1; 2–3; >3), and expert_j, project_j, targeted_project_j and redlist_j are dummies of a five-level factor denoting the source of the data. The source might either be a common naturalist observation (reference level), an observation by an expert naturalist, an observation made during a not further specified project, an observation made in a project targeted at the focal species or an observation made during a Red-List inventory. An expert naturalist was defined as an observer that contributed a significant number of records, which was defined as the upper 2.5% quantile of all observers arranged by their total number of records, and that made at least one visit with an exceptionally long species list, which was defined as a visit in the upper 2.5% quantile of all visits arranged by the number of records. The proportions of records originating from these different sources changed across years, but change was not unidirectional and differed among the investigated groups (Supplementary Fig. S11), such that accounting for data source in the model should suffice to yield reliable estimates of occupancy trends. α_d1,t is a random effect for year, which was modelled as

The occupancy and observation models were fitted jointly in Stan through the interface CmdStanR⁶⁸. Four Markov chain Monte Carlo chains with 2000 iterations each, including 1000 warm-up iterations, were used. Priors of the model parameters are specified in Supplementary Table S5. For the prior distribution of global intercepts, a standard deviation of 1.5 was chosen to not overweight the extreme values on the probability scale. To ensure that chains mixed well, Rhat statistics for annual mean occupancy estimates were calculated through the package rstan⁶⁹. For Switzerland-wide annual estimates (n = 18,140), 98.0% of values met the standard threshold of 1.1 (99.9% of values <1.55 and all values <1.56). For the annual estimates per bioclimatic zone (n = 116,844), 98.3% of values met the 1.1 threshold (99.9% of values <1.55 and all values <1.62). To check the validity of the model results, we compared species trends estimated for butterflies to trends estimated from standardised samplings for the same region in a national monitoring programme and found them to be clearly correlated (Supplementary Fig. S3).

From the models of all three insect groups, we extracted the posterior distribution of the predicted annual mean occupancy per bioclimatic zone and for the whole of Switzerland for all 390 species. We used linear models to quantify species trends of mean occupancy against years for each draw of the posterior distribution. On the one hand, we determined global species trends in mean occupancy across the whole of Switzerland for the time range 1980–2020. On the other hand, we determined short-term species trends in mean occupancy per bioclimatic zone for all consecutive 5-year intervals (8 intervals starting with 1980–1985). This was necessary to being able to analyse species trends against a set of uncorrelated climate and land-use change variables (see below). The length of the intervals was chosen to be ecologically meaningful while representing the variability in species trends adequately. Longer intervals might miss relations because short-term variability in drivers and trends is flattened out and analyses would lack replication of climate variables. Still, the main findings could be confirmed in additional analyses with 10-year intervals (Supplementary Table S6 and Supplementary Fig. S12).

Climate and land-use change

We selected climatic variables from the commonly used set of 19 bioclimatic variables⁷⁰ to represent climate change. Selection was restricted to variables that are potentially most meaningful to explain insect trends⁷¹ and not strongly correlated, resulting in three variables: annual mean temperature (BIO1) representing changes in absolute temperature; temperature seasonality (BIO4, defined as the standard deviation of monthly temperature means), representing changes in annual temperature cycles; and precipitation of the warmest quarter (BIO18) representing changes in precipitation during summer months (thus termed summer precipitation throughout). We determined yearly values for these variables for the whole of Switzerland from high-resolution mean monthly temperature and total monthly precipitation values reported by MeteoSwiss (https://www.meteoswiss.admin.ch) at a 1.25-degree minute grid (~2.3 × 1.6 km). Then, we calculated yearly means per bioclimatic zone while only considering grid cells intersecting with 1 × 1 km squares included in any of the occupancy-detection models. To infer climate change, we used linear models to determine trends in the three variables across consecutive 5-year intervals (same as for short-term species trends). To account for lags in effects of climate change⁷² and to prevent large biases in trend estimates due to single extreme years, we also included the 5 years preceding each 5-year interval in the linear models for climatic variables (spanning 10 years in total, first interval 1975–1985).

We derived land-use variables at regional scale from the national agricultural statistics and censuses that were recorded by the Federal Statistical Office for the time range 1955–2020, from which we used data on total agricultural area, grassland area, livestock numbers and area of different crops (arable and permanent). These data were available at the level of 2172 municipalities in yearly resolution (1996–2020) or 5- to 10-year resolution for earlier years (Supplementary Table S7). The cover of other land-use types, such as forests, might be important for the studied species, but data were not available at the necessary temporal resolution. As additional analyses of available data of other land-use types indicated that changes of the most relevant land-use types (forests, built-up areas) were strongly related to changes in agricultural area (Supplementary Fig. S8), inclusion of agricultural area in our analyses sufficiently covered the main land-use changes in the last 40 years in the study area. We summarised livestock numbers into number of livestock units (LSU) based on commonly used weighting factors⁷³. To aggregate data at the level of the nine bioclimatic zones, we distributed municipality-level data to agricultural land within municipalities based on the spatial distribution of agricultural areas within the municipalities⁷⁴. Note that data from the national censuses are attributed to the municipality of the farm of the landowner, resulting in some degree of misclassification, which is, however, expected to level out at the aggregation level used for later analyses. We then aggregated the resulting spatially explicit data to the bioclimatic zones, considering only squares that were included in any of the occupancy-detection models to not overrepresent land use in largely unvisited regions (mainly high alpine regions). Because all included variables have little year-to-year variation, but rather change gradually, we used generalised additive models to fill data gaps and to predict yearly values for the time range of interest (1980–2020) for all variables and bioclimatic zones (function stan_gamm4 in the package rstanarm⁷⁵, four Markov chain Monte Carlo chains with 2000 iterations each, including 1000 warm-up iterations). Based on these data, we derived three variables of land-use change at the level of bioclimatic zones. First, we determined the proportion of agricultural area as the proportion of total agricultural land in relation to the total study area. Second, we quantified grassland-use intensity as the mean number of LSU per grassland area. As such, it represents a general index of grazing, mowing and fertilisation pressure on the available grasslands. To check the validity of this measure, we compared trajectories of regional grassland-use intensity to grassland-use intensity inferred from satellite imagery, which showed very similar regional trends (Supplementary Fig. S13). Third, we determined crop-use intensity by attributing mean insecticide application rates to the most common crop types (including orchards and vineyards) (following ref. 76). Using these values, we calculated insecticide application rates for the total cropped area in a bioclimatic zone and aggregated it to an overall application rate per zone. To represent a mean insecticide pressure in the landscapes, we set the aggregated rates in proportion to the total study area per zone (note that here general additive models were only used on the aggregated application rate and not on the single crop types). Finally, we used changes across the investigated 5-year intervals as variables of land-use change.

Species traits

Various traits have been linked to species’ susceptibility to global change^{23,31,49,77,78,79,80,81}. Here, we included two traits in our analyses of species trends that were expected to be strongly sensitive to the investigated drivers (climate and land-use changes). First, the temperature niche of each species may determine its response to climate change. We estimated the temperature niche following the Species Temperature Index approach⁸² based on distributional records from the GBIF database (GBIF.org, https://doi.org/10.15468/dl.t6ha3h for butterflies, https://doi.org/10.15468/dl.reemkv for grasshoppers, https://doi.org/10.15468/dl.czbrmq for dragonflies). To reduce sampling biases, we only considered records from Europe and aggregated them at the Common European Chorological Grid Reference System (CGRS) grid. We extracted mean temperature (1970–2000) from WorldClim 2⁸³ at a 2.5 min spatial resolution and aggregated it at the CGRS grid. For each species, we determined the temperature niche as the mean temperature of the grid cells where it was recorded at least once. Second, we quantified habitat specialisation for all species, as specialist species, e.g., in terms of habitat or feeding preference, have repeatedly been shown to respond particularly strongly to land-use change and intensification^77,79,80. We estimated habitat specialisation based on a database of ecological preferences⁸⁴, which compiles information on habitat preference specific for Switzerland for all studied species based on available literature and expert knowledge. We extended the database to include Zygaenidae based on relevant literature^85,86 and expert knowledge. For each species, the preferred habitats are given from a list of typical habitats, which are defined separately for each of the studied insect groups based on their habitat ranges (Supplementary Table S8). For grasshoppers, adult habitat preference is specified, whereas larval habitat preference is specified for butterflies and dragonflies. We defined a continuous habitat specialisation index I_spec,i for species i as

where n_i is the number of preferred habitats of species i and N_group(i) is the total number of habitats listed for the insect group. To obtain a global measure for habitat specialisation for each single species, we standardised this index to values between 0 (least specialised) and 1 (most specialised) within insect groups.

Analysis of species trends

We analysed trends in species mean occupancy at 5-year intervals in a regression model with changes in climate and land use (at the same time intervals) as well as species traits as explanatory variables (Fig. 3a). Besides the main effects of the three climate change variables (annual mean temperature, temperature seasonality, summer precipitation), the three land-use change variables (total agricultural area, grassland-use intensity, crop-use intensity), the species traits (temperature niche, habitat specialisation) and the elevation (low or high), we included a set of interactions, which were expected to affect species trends (Fig. 3a). First, because absolute values of climate and land-use variables differ most considerably between the two elevations (Supplementary Fig. S9), we included interactive effects between the elevation and climate as well as land-use change variables. Because crop-use intensity in high elevations was minimal due to generally low amounts of cropped area, we did not include the interactive effect between elevation and crop-use intensity. Second, following the same line of reasoning as for global change variables, we included interactive effects of species traits and elevation. Third, we included all interactive effects of climate and land-use change variables to test for synergistic or antagonistic effects between the two global change drivers. Fourth, because the temperature niche of a species is expected to affect its reaction to climate change, we included the interactive effects between temperature niche and the three climate change variables. Finally, because specialised species have been shown to be particularly susceptible to land-use change and intensification, we included the interactive effects between habitat specialisation and the three land-use change variables. In addition, we included separate intercepts for the three groups to account for potential differences in mean trends and a factorial covariate for the 5-year time interval to account for variance in trends not covered by the climate and land-use change variables. Also, we included the bioclimatic zone and species identity as random effects. For species identity, we not only included random intercepts, but also random slopes in respect to the global change drivers and their interactions, because different species were expected to react differently to these drivers. To meet the normality assumption for the residual distribution, we transformed species trends by taking the square root of their absolute values, while keeping their original sign. We scaled all continuous variables to standard deviation 1 and centred to mean 0 at the level of their recording prior to analyses. We fitted the model in Stan through the R interface rstan⁶⁹ (four Markov chain Monte Carlo chains with 2000 iterations each, including 1000 warm-up iterations; priors in Supplementary Table S5). To assess the model fit, we did posterior predictive model checking of residual distribution and checked for temporal autocorrelation in the residuals (Supplementary Fig. S14).

Not all study species are equally associated to agriculturally influenced habitats and might thus not be equally sensitive to agricultural land-use change. While crop-use intensity is expected to affect all habitats in a region due to the spread of insecticides to neighbouring terrestrial habitats⁸⁷ or water bodies⁸⁸, the effect of the other land-use change variables (total agricultural area, grassland-use intensity) on habitats that are not agriculturally influenced is less clear. As a result, we used three different versions of the above model to account for a potential bias. In the first version of the model, we made sure that all parameters that included the respective land-use change variables (total agricultural area, grassland-use intensity) were estimated based on observations from a subset of species, namely species that are at least partly associated to agriculturally influenced habitats. All other parameter estimates, however, were based on observations from the full set of species. The subset of species was defined based on the database of ecological preferences that was also used to quantify habitat specialisation⁸⁴. From the list of habitats, we defined agriculturally influenced habitats (Supplementary Table S8) and included only species that occur in at least one of them in the species subset. It contained 183 butterfly species (out of 215), 93 grasshopper species (out of 103) and no dragonfly species. To avoid underestimating effects of land-use changes due to the inclusion of species not directly associated to agriculturally influenced habitats when comparing climate and land-use effects (see scenario predictions below), the second version of the model did only include the species associated to agriculturally influenced habitats for all parameter estimates, i.e., was a full model run on a subset of the data. However, changes in an agricultural area or grassland-use intensity might also affect the species not directly associated to agriculturally influenced habitats. For example, a change in agricultural area indirectly also means a change in other habitats such as forests (Supplementary Fig. S8) and increasing grassland-use intensity might result in higher nutrient loads of adjacent water bodies⁸⁹. Thus, the third version of the model was a full model with all species, where the full set of species was included for all parameter estimates. The findings based on the different model versions were largely consistent. We report parameter estimates from the first version in the main manuscript and results from the other versions in Supplementary Tables S3 and S4.

Additional sensitivity analyses were done based on the first model version. Because the data-based species selection also included some species for which mean occupancy estimates might be biased due to different reasons (migratory or (re)introduced species and species with uncertain taxonomic status or very difficult identification; Supplementary Table S1), we did sensitivity analyses, in which we excluded these species from the trend analyses. At the same time, mean occupancy estimates might be less reliable for species with only few records. Thus, we also excluded the species with the lowest record numbers (lowest 20% of species per group) for these sensitivity analyses. We found model results to be robust to these exclusions (Supplementary Table S9). Furthermore, mean occupancy estimates for some species in some bioclimatic zones might be unreliable because of low record numbers in these zones. In additional sensitivity analyses, we thus excluded species-zone combinations with very few records (<41, i.e., on average less than 1 record per year and zone) and again found model results to be robust (Supplementary Table S10).

To understand how the relations of short-term species trends to global change drivers explained the observed changes in species’ mean occupancy across the 40 studied years (long-term trends), we predicted 40-year species trends for different climate and land-use change scenarios from the parameter estimates of the regression model. To prevent an underestimation of the land-use change effects, we used the second model version (i.e., only including species associated to agriculturally influenced habitats) for these model predictions. First, the basic scenario was composed as such that, starting from the original dataset, all climate and land-use change variables were held constant at their absolute 0 (uncentred variables), representing no change. All other variables (elevation, time interval, bioclimatic zone, temperature niche, habitat specialisation) were kept unchanged. Second, we defined additional scenarios on top of the basic scenario such that always one climate or land-use change variable was following its measured trajectory. Third, to account for interactive effects between climate and land-use change, we defined scenarios in which always two variables (one climate change, one land-use change) followed their measured trajectories. Finally, we defined a scenario in which all variables followed their measured trajectories. From all scenarios, we made predictions at the replication unit of the original dataset (one prediction per species and bioclimatic zone and time interval). We then back-transformed these short-term trend predictions (reversed root transformation) to represent mean occupancy changes across a 5-year interval, summed them per species and bioclimatic zone to represent 40-year trends and finally averaged long-term trends per species while accounting for the number of squares per bioclimatic zone to yield an estimate of the Switzerland-wide mean occupancy change across 40 years for each species. We applied the same procedure to the observed short-term trend estimates. To evaluate how well the different scenarios reflected the observed trends, we determined R² values (squared Pearson correlations) for the comparison of predicted and observed long-term trends as a measure of explained variance for all scenarios.

Source: Ecology - nature.com