More than 75 percent decline over 27 years in total flying insect biomass in pro...
source link: https://journals.plos.org/plosone/article?id=10.1371%2Fjournal.pone.0185809
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
Materials and methods
Biomass data.
Biomass data were collected and archived using a standardized protocol across 63 unique locations between 1989 and 2016 (resulting in 96 unique location-year combinations) by the Entomological Society Krefeld. The standardized protocol of collection has been originally designed with the idea of integrating quantitative aspects of insects in the status assessment of the protected areas, and to construct a long-term archive in order to preserve (identified and not-identified) specimens of local diversity for future studies. In the present study, we consider the total biomass of flying insects to assess the state of local entomofauna as a group.
All trap locations were situated in protected areas, but with varying protection status: 37 locations are within Natura2000 sites, seven locations within designated Nature reserves, nine locations within Protected Landscape Areas (with funded conservation measures), six locations within Water Protection Zones, and four locations of protected habitat managed by Regional Associations. For all location permits have been obtained by the relevant authorities, as listed in the S1 Appendix. In our data, traps located in nutrient-poor heathlands, sandy grasslands, and dune habitats provide lower quantities of biomass as compared to nutrient nutrient-rich grasslands, margins and wastelands. As we were interested in whether the declines interact with local productivity, traps locations were pooled into 3 distinct habitat clusters, namely: nutrient-poor heathlands, sandy grassland, and dunes (habitat cluster 1, n = 19 locations, Fig 1A), nutrient-rich grasslands, margins and wasteland (habitat cluster 2, n = 41 locations, Fig 1B) and a third habitat cluster that included pioneer and shrub communities (n = 3 locations).
Most locations (59%, n = 37) were sampled in only one year, 20 locations in two years, five locations in three years, and one in four years, yielding in total 96 unique location-year combinations of measurements of seasonal total flying insect biomass. Our data do not represent longitudinal records at single sites, suitable to derive location specific trends (e.g. [28]). Prolonged trapping across years is in the present context (protected areas) deemed undesirable, as the sampling process itself can negatively impact local insect stocks. However, the data do permit an analysis at a higher spatial level, i.e. by treating seasonal insect biomass profiles as random samples of the state of entomofauna in protected areas in western Germany.
Malaise traps were deployed through the spring, summer and early autumn. They operated continuously (day and night), and catches were emptied at regular intervals, on average every 11.2 days (sd = 6.3). We collected in total 1503 trap samples, with an average of 16 (4–35) successive catches per location-year combination (Table 1). Between 1989 and 2016, a total of 53.54kg of invertebrates have been collected and stored, over a total trap exposure period of 16908 days, within an average of 176 exposure days per location-year combination. Malaise traps are known to collect a much wider diversity of insect species (e.g. [29–31]) as compared to suction traps (e.g. [28]) and are therefore considered superior as a method of collecting flying insects. On the basis of partial assessments, we can assume that the total number of insects included in 53.54 kg biomass represents millions of individuals.
For each year, the number of locations sampled, the number of location re-sampled, total number of samples, as well as mean and standard deviation of exposure time at the trap locations (in days) are presented.
The sampling was standardized in terms of trap construction, size and design (identical parts), colors, type of netting and ground sealing, trap orientation in the field as well as slope at the trap location. Hence none of the traps differed in any of these field aspects. Our trap model was similar to the bi-colored malaise trap model by Henry Townes [32, 33]. The traps, collecting design, and accompanying methods of biomass measurement as designed and applied by the Entomological Society Krefeld are described elsewhere [34–36] and in S2 Appendix.
Trap catches were stored in 80% ethanol solution, prior to weighing, and total insect biomass of each catch (bottle) was obtained based on a standardized measurement protocol by first subtracting fluid content. In order to optimally preserve samples for future species determination, the insects were weighed in an alcohol-wet state. First, the alcohol concentration in the vessels was stabilized to 80%, while this concentration was controlled with an areometer over a period of at least two days. In order to obtain biomass per sample with sufficient accuracy and comparability, the measuring process was fixed using a standardized protocol [34]. For this purpose the insects of a sample were poured onto a stainless steel sieve (10cm diameter) of 0.8 mm mesh width. This sieve is placed slightly obliquely (30 degrees) over a glass vessel. The skew position accelerates the first runoff of alcohol and thus the whole measuring procedure. The drop sequence is observed with a stopwatch. When the time between two drops has reached 10 seconds for the first time, the weighing process is performed with a laboratory scale. For the determination of the biomass, precision scales and analytical scales from Mettler company were used with an accuracy of at least 0.1g and controlled with calibrated test weights at the beginning of a new weighing series. In a series of 84 weightings of four different samples repeating this measurement procedure, an average deviation from the mean value of the measurement results of 0.4 percent was observed (unpublished results).
Weather data.
Climate change is a well-known factor responsible for insect declines [15, 18, 21, 37]. To test if weather variation could explain the observed decline, we included mean daily temperature, precipitation and wind speed in our analysis, integrating data from 169 weather stations [38] located within 100km to the trap locations. We examined temporal trends in each weather variable over the course of the study period to assess changes in climatic conditions, as a plausible explanation for insect decline. Estimates of each weather variable at the trap locations were obtained by interpolation of each variable from the 169 climate stations.
We initially considered mean daily air temperature, precipitation, cloud cover, relative air moisture content, wind speed, and sunshine duration. However, only temperature, precipitation and wind speed were retained for analysis, as the other variables were significantly correlated with the selected variables [R(temperature, cover) = −43.2%, R(temperature, sunshine) = 53.4%, R(precipitation, moisture) = −47.3%] and because we wanted to keep the number of covariates as low as possible. Additionally, we calculated the number of frost days and the sum of precipitation in the months November- February preceding a sampling season. We used spatio-temporal geostatistical models [39, 40] to predict daily values for each weather variable to each trap location. Amongst other methods, the geostatistical approach is considered a superior interpolation method in order to derive weather variables to trap locations [41]. Uncertainty in interpolated variables such as wind speed is usually associated with altitude differences. However, as our trap locations are all situated in lowland areas with little altitude variation, we do not expect a large error in our interpolations at trap locations.
We decomposed the daily values of each weather variable into a long-term average trend (between years), a mean seasonal trend, and a yearly seasonal anomaly component (S2 Fig), modeled using regression splines [42] while controlling for altitude of weather stations. The remaining residual daily values of each station were further modeled using a spatio-temporal covariance structure. For example, temperature T, on given day t, of a given year k at a given trap location s is modeled as: (1) where fk(k) is the long-term trend over the years (a thin plate regression spline), ft(t) the mean seasonal trend within years (a penalized cyclic cubic regression spline), r(k, t) the mean residual seasonal component, which measures annual anomaly in mean daily values across selected stations, and a is the linear coefficient for the altitude h effect. The spatio-temporal covariance structure Cs, t, fitted independently to the residuals of each weather variable model, allowed us to deal with lack of independence between daily weather data within and between stations, as well as to interpolate to trap locations using kriging. Altitude of trap locations was extracted from a digital elevation models at 90m resolution [43].
Land use data.
Land use variables (and changes therein) were derived from aerial photographs [44] taken within two distinct time periods (between 1989–1994, and between 2012–2015), and allowed us to characterize land use composition at surroundings of the traps, as well as changes over time. We distinguished cover of forests, agricultural areas, natural grassland, and surface water. For each trap location, aerial photographs were manually processed, polygons extracted and categorized, and their surface area calculated with a radius of 200 meter. Preliminary analysis of the relationship between log biomass and landuse variables, on a subset of the trap locations, indicated that land use elements at 200m radius better predicted insect biomass than elements at 500 and 1000m radius, similar to findings elsewhere for wild bees [45]. Land use variables were measured at a coarse temporal resolution, but fortunately cover the temporal span of insect sampling. To link the cover of a given land use variable to the insect biomass samples in a particular year, we interpolated coverage between the two time points to the year of insect sampling using generalized linear models with a binomial error distribution, a logit link, and an estimated dispersion parameter. Mean distributions of land use at each of the two time points are depicted in S3A & S3B Fig.
Habitat data.
Plant inventories were conducted in the immediate surroundings (within 50m) of the trap, in the same season of insect sampling. These data permitted the assessment of plant species richness (numbers of herbs, shrubs and trees) and environmental conditions based on average Ellenberg values [46–48], as well as changes therein over time. Each Ellenberg indicator (we considered nitrogen, pH, light, temperature and moisture) was averaged over all species for each location-year combination. We examined annual trends in each of the above-mentioned variables in order to uncover potential structural changes in habitat characteristics over time. Species richness was analyzed using mixed-effects generalized linear models [49] with a random intercept for trap location and assuming a Poisson distribution for species richness, and a normal distribution for mean Ellenberg indicator values. Although a Poisson distribution fitted tree and shrub species adequately, (residual deviance/ degree of freedom = 0.94 and 1.04 respectively), severe overdispersion was found for herb species richness (residual deviance/ degree of freedom = 2.16). Trend coefficients of richness over time between a Poisson mixed effects model and a negative binomial model were comparable but differed in magnitude (Poisson GLMM: −0.034 (se = 0.003), vs NB GLMM −0.027 (se = 0.006)). Although the fit is not perfect in the case of herb richness, we believe our trend adequately describes direction of change over time. Mean changes in plant species richness are depicted in S3C Fig.
Insect biomass model
The temporal resolution of the trap samples (accumulated over several days) is not directly compatible with the temporal distribution of the weather data (daily values). Additionally, variable exposure intervals between trap samples is expected to induce variation in trapped biomass between samples, and hence induce heteroscedasticity. Furthermore, biomass data can numerically only be positive on the real line, and we require a model to reflect this property of the data. Because of the unequal exposure intervals however, log-transforming the response would be inappropriate, because we require the sum of daily values after exponentiation, rather that the exponent of the sum of log-daily biomass values. In order to indirectly relate biomass to daily weather variables, to account for the variation in time exposure intervals over which biomass was accumulated in the samples, and to respect the non-negative nature of our data, we modeled the biomass of each catch as the sum of the expected (but unobserved) latent daily biomass. The mass m of each sample j, at site s in year k, is assumed to be distributed normally about the sum of the latent expected daily mass (zt, s, k), with variance : (2) subject to where τ1 and τ2 mark the exposure interval (in days) of biomass collection of each sample j. The latent daily biomass itself is represented by a log normal distribution, in which coefficients for covariates, random effects and residual variance are all represented on the log scale. In turn, daily biomass is modeled as (3) (4) where c is a global intercept, X a design matrix of dimensions n×p (number of samples × number of covariates; see Model analysis below), βx the corresponding vector of coefficients that measure the weather, habitat and land use effects, and log(λ) a mean annual population growth rate parameter. The random term (us) denotes the location-specific random effect assumed to be distributed normally about zero . The exponentiation of the right hand side of Eq (3) ensures expected values to be positive.
The expected residual variance of each sample , is expressed as the sum of variances of daily biomass values ().
(5)The variances of daily biomass should respect the non-negative nature of the data as well. Additionally, we are interested in being able to compare the residual variance with the random effects variance, and this requires them to be on the same scale. Therefore, we expressed the variance of the daily biomass as a function of the variance of the logarithm of the daily biomass. Using the method of moments: (6) where v represents the residual variance of daily log-biomass.
Analysis
We developed a series of models each consisting of a set of explanatory variables that measure aspects of climate, land use and local habitat characteristics. Significant explanatory variables in these models were combined into a final model, which was then reduced to exclude insignificant effects. An overview of which covariates were included in each model is given in Table 2.
The year covariate yields the annual trend coefficient.
Weather effects explored were daily temperature, precipitation and wind speed, as well as the number of frost days and sum of precipitation in the preceding winter. Habitat effects explored tree and herb species richness, as well as average Ellenberg values for nitrogen, pH, light, temperature and moisture, per location-year combination. Land use effects explored the fractions of agricultural area, forest, grass, and surface water in a radius of 200m around the plot location.
Parameter values are obtained by the use of Markov chain Monte Carlo (MCMC) methods by the aid of JAGS (Just Another Gibbs Sampler [50]) invoked through R [51] and the R2Jags package [52]. JAGS model scripts are given in S1 Code, while data are given in S1 and S2 Dataset. For each model, we ran 3 parallel chains each consisting of 24000 iterations (first 4000 discarded), and kept every 10th value as a way to reduce within chain autocorrelation. We used vague priors for all parameters, with uniform distributions for the residual and random effect variance components, and flat normal distributions (with very high variance) for all other parameters. Covariates in X were standardized prior to model fitting, with the exception of year (values 1–26), and land use variables (proportions within 0–1 range).
For all models, we computed the Deviance Information Criterion [53] (DIC) as well as the squared correlation coefficient (R2) between observed and mean posterior estimates of biomass on the log scale. Results are given in Table 3. Parameter convergence was assessed by the potential scale reduction factor [54] (commonly ), that measures the ratio of posterior distributions between independent MCM chains (in all models, all parameters attained values below 1.02). For all models, we confirmed that the posterior distribution of the trend coefficient did not confound any other variable by plotting the relevant posterior samples and computing pairwise correlation coefficients.
For each model, the number of parameters, the Deviance Information Criterion, the effective number of parameters (pD), calculated R2 and difference in DIC units between each model and the model with lowest ΔDIC. See Table 2 for covariates included in each model.
Our basic model included habitat cluster (3 levels), a quadratic effect for day number, an annual trend coefficient measuring the rate of biomass change, and the interactions between the annual trend coefficient and the day number variables. Next we developed 3 models each consisting of either weather variables (S1 Table), land use variables (S2 Table), or habitat variables. Because interactions between the annual rate of change and land use variables seemed plausible, a fourth model was developed to include these interactions (S3 Table). Finally, all significant variables were combined into our final model (Table 4), which included effects of an annual trend coefficient, season (linear and quadratic effect of day number), weather (temperature, precipitation, number of frost days), land use (cover of grassland and water, as well as interaction between grassland cover and trend), and habitat (number of herb and tree species as well as Ellenberg temperature).
For each included variable, the corresponding coefficient mean, standard deviation and 95% credible intervals are given. P-values were calculated empirically based on posterior distributions of coefficients.
Our estimate of decline is based on our basic model, from which we can derive seasonal estimates of daily biomass for any given year. The basic model includes only a temporal (annual and seasonal effects, as well as interactions) and a basic habitat cluster distinction (additive effects only) as well as a random trap location effect. We here report the annual trend coefficient, as well as a weighted estimate of decline that accounts for the within season differences in biomass decline. The weighted insect biomass decline was estimated by projecting the seasonal biomass (1-April to 30-October) for years 1989 and 2016 using coefficients our basic model, and then dividing the summed (over the season) biomass of 2016 by the summed biomass over 1989.
Using our final model, we assessed the relative contribution (i.e. net effect) of the explanatory variables to the observed decline, both combined and independently. To this aim we projected the seasonal daily biomass for the years 1989 and 2016 twice: first we kept covariates at their mean values during the early stages of the study period, and second we allowed covariate values to change according to the observed mean changes (see S2 and S3 Figs). Difference in the total biomass decline between these two projections are interpreted as the relative contribution of the explanatory variables to the decline. The marginal (i.e. independent) effects of each covariate were calculated by projecting biomass increase/decline as result of the observed temporal developments in each variable separately, and expressing it as percentual change.
Our data provide repetitions across years for only a subset of locations (n = 26 out of 63). As such, spatial variation in insect biomass may confound the estimated trend. To verify that this is not the case, we fitted our basic model (but excluding the day number and year interaction to avoid overparameterization) to the subset of our data that includes only locations that were sampled in more than one year. Seasonal profiles of daily biomass values are depicted in S4 Fig. Finally, we reran our basic model for the two (of the three) habitat clusters (for which sufficient data existed; see Biomass Data) separately in order to compare the rate of decline between them (S5 Fig).
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK