Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: Statistical analyses using field census data from 2009 to 2016

Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: Statistical analyses using field census data from 2009 to 2016

Journal Pre-proofs Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: statistical analyses usin...

644KB Sizes 0 Downloads 7 Views

Journal Pre-proofs Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: statistical analyses using field census data from 2009 to 2016 Kosuke Nakanishi, Tetsuyuki Uéda, Hiroyuki Yokomizo, Takehiko I. Hayashi PII: DOI: Reference:

S0048-9697(19)34490-0 https://doi.org/10.1016/j.scitotenv.2019.134499 STOTEN 134499

To appear in:

Science of the Total Environment

Received Date: Revised Date: Accepted Date:

21 July 2019 10 September 2019 15 September 2019

Please cite this article as: K. Nakanishi, T. Uéda, H. Yokomizo, T.I. Hayashi, Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: statistical analyses using field census data from 2009 to 2016, Science of the Total Environment (2019), doi: https://doi.org/10.1016/j.scitotenv.2019.134499

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Elsevier B.V. All rights reserved.

Title:

Effects of systemic insecticides on the population dynamics of the dragonfly Sympetrum frequens in Japan: statistical analyses using field census data from 2009 to 2016

Kosuke Nakanishi a*, Tetsuyuki Uéda b, Hiroyuki Yokomizo a, Takehiko I. Hayashi a

a

National Institute for Environmental Studies, Onogawa 16-2, Tsukuba, Ibaraki

305-8506, Japan b

Faculty of Bioresources and Environmental Sciences, Ishikawa Prefectural University,

Suematsu 1-308, Nonoichi, Ishikawa 921-8836, Japan

*Corresponding

author

[email protected] Tel.: +81-29-850-2144 Fax: +81-29-850-2920 ABSTRACT

Since the mid-1990s, populations of the common Japanese dragonfly Sympetrum frequens in rice fields have declined severely. Application of systemic insecticides—especially fipronil—to nursery boxes of rice seedlings is suspected to be the main cause of the decline. However, until now there have been insufficient population data to test the causality. We conducted a dragonfly survey from 2009 to 2016 in four prefectures of Japan and compiled the data to enable the comparison of

1

population growth rates along five main census routes over the years. We also estimated the use ratio of each insecticide applied to nursery boxes in rice fields (i.e., the area exposed to insecticide as a ratio of the total rice field area) by prefecture. We then statistically analyzed the effects of the insecticides on the dragonfly’s population growth rates, taking into account the potential confounding factors based on current knowledge. There was a significant negative association between the annual increase in use ratio of the sum of systemic insecticides (e.g., fipronil, imidacloprid, and chlorantraniliprole) and the annual population growth rate of S. frequens. This association suggests that systemic insecticide use affected the decrease in population density of the dragonfly, although some agronomic factors need to be further examined as potential confounders.

Keywords: Backdoor criterion; Confounding factor; Neonicotinoid; Odonata; Paddy field; Phenylpyrazole

2

1. Introduction Rice fields, including irrigation water systems, play important roles as alternatives to natural wetlands in sustaining biodiversity (Katayama et al., 2015; Lawler, 2001; Natuhara, 2013). There is, however, a trade-off between biodiversity conservation and rice productivity (e.g., Katayama et al., 2019). In general, insecticide use decreases biodiversity in rice ecosystems because of its adverse effect on a wide range of non-target organisms. Among many kinds of insecticides, the systemic insecticides neonicotinoids and fipronil (a phenylpyrazole) are the most widely used worldwide (Simon-Delso et al., 2015), although some rice pests have been resistant to these insecticides in recent years (reviewed by Furlan et al., 2018). Neonicotinoids alone accounted for more than 25% of the world insecticide market in 2014 (Bass et al., 2015). Because these insecticides are highly toxic to non-target organisms, they can have severe negative impacts on biodiversity and rice field ecosystems (reviewed by Pisa et al., 2015, 2017). Dragonflies (Odonata) at both the adult and the aquatic nymphal stage are common predators in rice fields. In Japan, rice fields occupy approximately 24,000 km2 (about half of the total cultivated area) (Ministry of Agriculture, 2018). These rice fields provide one of the main habitats (particularly breeding sites) for dragonflies. As many as 31 dragonfly species of a total of approximately 200 in Japan use rice fields (Uéda, 1998). In recent years, however, many dragonfly species are suspected to have been threatened with extinction because of degradation of agricultural habitats (Kadoya et al., 2009). Sympetrum frequens (Odonata: Libellulidae) is the species most commonly breeding in rice fields in Japan (Uéda, 1998). Population densities of this species decreased

3

drastically (e.g., to 1/100th of their levels in the early 1990s) in many regions of Japan in the 1990s (Futahashi, 2012; Uéda and Jinguji, 2013). Fipronil and imidacloprid (a neonicotinoid) are suspected to be the main causes of these population declines, because an association between increased use of these insecticides and population declines of S. frequens has been observed (Uéda and Jinguji, 2013) and these insecticides are lethally toxic to Sympetrum species (Jinguji et al. , 2009, 2013). In conventional Japanese rice cultivation, these insecticides are frequently used for nursery-box treatment of rice seedlings. Application of fipronil and imidacloprid to nursery boxes of rice seedlings is suspected to have the most severe impact on S. frequens because of the application timing and the dragonfly’s univoltine life cycle (Nakanishi et al., 2018); newly hatched nymphs are exposed to high concentrations of insecticides released from the rice seedlings just after transplantation (Sánchez-Bayo and Goka, 2006; Thuyet et al., 2013). The other types of insecticides applied to rice fields (i.e., aerial application) are not likely to affect mortality of the dragonfly, because when the insecticides are applied, S. frequens is in the adult stage and inhabiting mountains far away from rice fields. Nakanishi et al. (2018) estimated insecticide use for nursery box application during the period of sharp decline in S. frequens populations in the 1990s. They found a significant negative association between fipronil use and abundance of the dragonfly. However, this analysis used only one available data set from one site (Toyama Prefecture) for the decline period between 1993 and 2004. Because of the lack of replicated years and sites they were unable to statistically exclude the effects of possible confounding factors (e.g., year- and site-specific effects) and sampling variability. Therefore, the association they presented was not regarded as strong statistical evidence to support the hypothesis that fipronil and imidacloprid caused the population declines

4

of S. frequens in the 1990s. To estimate the causal effects of insecticides, we need to exclude confounding effects in statistical models. “Common cause” factors that affect both causal factors (i.e., insecticide application in this case) and the outcome (i.e., dragonfly’s population growth rates) can lead to confounding biases in causal estimations (Hernán and Robins, 2019). In estimations of causal effects, these common cause variables should be included (or, more unambiguously, the backdoor criterion should be satisfied; Pearl, 2009; Pearl and Mackenzie, 2018) in statistical models so as to control for confounding. “Common cause” variables can be formally identified by creating a causal diagram that represents the data-generation mechanism behind the observed data (Pearl, 2009). Here, we aimed to estimate the causal effects of insecticide use on nursery boxes on the population dynamics of S. frequens. We measured dragonfly abundance from 2009 to 2016 in five prefectures in north-central Japan. We then compiled the data to enable us to compare the dragonfly’s yearly population growth among five main census routes. We also estimated, by prefecture, the annual amounts of nine major insecticides, including fipronil and neonicotinoids, used for nursery box application. We then developed statistical models to estimate the causal effects of the insecticides, selecting covariates to control for confounding. This selection was based on a working causal diagram that we developed from currently available ecological and agricultural background knowledge. The relationships between the use of the insecticides and the dragonfly’s population growth rates were then estimated by using these statistical models.

5

2. Methods 2.1. Dragonfly species Sympetrum frequens is a common univoltine dragonfly distributed throughout Japan (Sugimura et al., 1999). It reproduces mainly in rice paddies between September and November, after the rice harvest (Uéda, 1988). During the breeding season, pairs of S. frequens fly to ovipositing sites from their roosts in tandem posture, mainly in the morning. Female S. frequens lay eggs in puddles in the rice paddies. The eggs overwinter and hatch the following spring in about April, soon after the start of irrigation and seedling transplantation. The nymphs emerge as adults in June and July and then migrate to mountainous areas far away from the rice fields. In late September, the adults return to the rice fields to start breeding (Inoue and Tani, 2010).

2.2. Dragonfly census We conducted a field survey to monitor the dragonfly abundance in five adjacent prefectures in north-central Japan, namely Niigata, Toyama, Ishikawa, Fukui, and Nagano. This survey included various census routes within and across years to reveal the status of population declines at various sites in these prefectures. We drove a car from 08:30 to 12:00 along country roads in rice-field areas, counting the numbers of tandem pairs of S. frequens passing anterior to the front window. We performed this route census on non-rainy days between mid-September and November (the breeding season of S. frequens in the study area; Uéda, 1988) from 2009 to 2016 (see Supplementary note S1 for details of the census method). Our original data from the above survey contained data from various census routes for different years. This made it difficult to directly compare data among different years,

6

because the base of comparison (i.e., the census routes) was not consistent across years. Therefore, by selecting data only from the main census routes that were consistent over years, we recompiled the data so that we could compare yearly population growth. All of the original census data will be reported separately (T. Uéda, in preparation). We used data compiled for five main census routes in four prefectures: route 1 for the cities of Joetsu and Myoko in Niigata Prefecture; route 2 for the cities of Tonami and Nanto in Toyama Prefecture; route 3 for the cities of Nonoichi and Hakusan in Ishikawa Prefecture; route 4 for the cities of Katsuyama and Ono in Fukui Prefecture; and route 5 for the cities of Fukui, Sabae, and Echizen in Fukui Prefecture (Fig. 1). We could not conduct the census on routes 2 and 5 in 2009 or on route 1 in 2016. Each census route was not always exactly the same because of changes in traffic and road conditions. The number of censuses for each route and each year ranged from 1 to 9 (Table 1), and the total length of each census route ranged from 34 to 118 km. We standardized the observed numbers of tandem pairs by the total route length of each census. We used these numbers (i.e., pairs/km) per census as a dragonfly population density index for the analyses.

132.3. Insecticide use We estimated the annual amounts of insecticides applied to nursery boxes of rice seedlings by prefecture as described by Nakanishi et al. (2018) and in accordance with the method of Yachi et al. (2016, 2017). Among various insecticides, we chose the nine dominant systemic insecticides used from 2009 to 2016 for the analyses: the neonicotinoids imidacloprid, dinotefuran, clothianidin, and thiamethoxam; the phenylpyrazole fipronil; the carbamates benfuracarb and carbosulfan; the thiocarbamate

7

cartap; and the anthranilic diamide chlorantraniliprole. Details of the estimating procedure are given in Supplementary note S2.

2.4. Models and analyses 2.4.1. Population model of a basic relationship between insecticides and dragonfly population growth We used the following population model of a basic relationship between the insecticide use ratio and the population growth rate of S. frequens:

𝑙=𝐿

λ𝑗,𝑖,𝑡 = 𝛽1∆𝐼𝑁𝑆𝐸1,𝑗,𝑡 +⋯ + 𝛽𝑘∆𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 + ∑𝑙 = 1𝛾𝑙∆𝑂𝐹𝑙,𝑗,𝑖,𝑡,

(Eq. 1)

where λ𝑗,𝑖,𝑡 is the annual population growth rate on census route i in prefecture j in year t. λ𝑗,𝑖,𝑡 is defined by λ𝑗,𝑖,𝑡 = ln (𝑁𝑗,𝑖,𝑡/𝑁𝑗,𝑖,𝑡 ― 1) (i.e., an equivalent expression in difference form is λ𝑗,𝑖,𝑡 = ln 𝑁𝑗,𝑖,𝑡 ― ln 𝑁𝑗,𝑖,𝑡 ― 1), where 𝑁𝑗,𝑖,𝑡 is the population density index (number of pairs/km) of S. frequens in year t on census route i in prefecture j. 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 is the use ratio of insecticide k in year t in prefecture j to which route i belongs. 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 was available only by prefecture. ∆𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 is the difference in 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 between year t and t – 1, defined by ∆𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 = 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 ― 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 ― 1. 𝛽𝑘 represents their coefficient. ∆𝑂𝐹𝑙,𝑗,𝑖,𝑡 represents the total differences in factors other than insecticides (l = 1 … L) between years t and t – 1 on census route i in prefecture j, and 𝛾𝑙 represents their coefficient. Factors that were considered in this term are described later. Note that Eq. 1 is a population model (i.e., not a statistical one; these are described later) and thus has no error term. The use of year-difference in 𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 and

8

𝑂𝐹𝑙,𝑗,𝑖,𝑡 in Eq. 1 comes from the modeling of ln 𝑁𝑗,𝑖,𝑡; that is, ln 𝑁𝑗,𝑖,𝑡 = 𝛽1𝐼𝑁𝑆𝐸1,𝑗,𝑡 𝑙=𝐿

+⋯ + 𝛽𝑘𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡 + ∑𝑙 = 1𝛾𝑙𝑂𝐹𝑙,𝑗,𝑖,𝑡. Biological interpretation of the use of year-differences is described in Supplementary note S3.

2.4.2. Selection of potential confounding factors for estimation of causal effects Confounding biases must be controlled for in statistical analyses in order to estimate the causal effects of explanatory variables (Hernán and Robins, 2019; Morgan and Winship, 2014). A confounding bias is caused by a “common cause” (or more unambiguously, an “unblocked backdoor path”; see Hernan and Robins, 2019, ch. 7) between a “cause” and a “response” variable (i.e., the insecticide use ratio and population growth rate, respectively). In general, common cause variables (hereafter, we call these variables “confounders”) must be included in statistical models to control for confounding. To specify confounders, we developed a working causal diagram (Fig. 2) that took into consideration the currently available ecological and agricultural background knowledge (described in Supplementary note S4) of insecticide use and the population dynamics of S. frequens. Note that our intention here was not to develop a perfect causal diagram, but to develop a working one to distinguish the factors that needed to be controlled for, and those did not need to be controlled for, in our causal analyses based on currently available knowledge. The working causal diagram specified at least three main potential confounders (Fig. 2; Supplementary note S4). First, “year” was a potential main confounder (Fig. 2). We therefore included year as an explanatory dummy variable to control for potential confounding bias in our statistical models. “Site” is also a potential main confounder (Fig. 2). In our analyses, however, we did 9

not include site factors (e.g., “prefecture” or “route”) themselves as explanatory 𝑙=𝐿

variables in the statistical models. This is because, in considering ∑𝑙 = 1𝛾𝑙∆𝑂𝐹𝑙,𝑗,𝑖,𝑡 in Eq. 1, the differences in any time-invariant factor between two consecutive years are canceled out (i.e., SITEj,i, t – SITEj,i, t–1 = 0 for site effects, where SITEj,i, t represents the site-specific factors in year t on route i in prefecture j). We assumed that the differences in use ratios (∆𝐼𝑁𝑆𝐸𝑘,𝑗,𝑡) and growth rates (λ𝑗,𝑖,𝑡) did not depend on site-specific factors, even though the absolute values of the use ratios and of the population densities may have depended on site-specific factors. The other potential main confounder was “agronomic practices” (e.g., midsummer drainage, crop rotation, and organic farming; Nakanishi et al., 2018) (Fig. 2). As far as we know, there were no apparent and systematic changes in agronomic practices in our survey area during our census period, although there were no available data on changes in agronomic practices. We therefore could not incorporate variables about local agronomic practices into the statistical models. We did not include factors related to biological interactions that can directly affect population growth rates of the dragonfly (e.g., variables related to prey and predator organisms). This was because these factors would not affect insecticide use and thus would not be confounders for the purposes of causal estimation. In controlling for confounding, these intermediate variables affected by the causal factor (i.e., insecticide use) itself should not be included in models, because their inclusion can cause bias in the estimation of causal effect (i.e., violation of condition (i) of the backdoor criterion; see Pearl, 2009, p. 79).

10

2.4.3. Statistical models On the basis of the above consideration of confounding factors and the population model (Eq. 1), we constructed a statistical model (model 1) for the multiple regression analysis:

λ𝑗,𝑖,𝑡 = 𝛼 + 𝛽𝑁𝑁𝐹𝑃∆𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡 + 𝛽𝐶𝐴𝑅𝑇∆𝐼𝑁𝑆𝐸𝐶𝐴𝑅𝑇,𝑗,𝑡 + 𝛽𝐶𝐻𝐿𝑂∆𝐼𝑁𝑆𝐸𝐶𝐻𝐿𝑂,𝑗,𝑡 + (Eq. 2) 𝛾𝑌𝐸𝐴𝑅𝑡𝑌𝐸𝐴𝑅𝑡 + 𝜀𝑗,𝑖,𝑡,

where 𝜀𝑗,𝑖,𝑡 represents an error term on route i in prefecture j in year t and α is an intercept. In this statistical model, we used three variables for insecticides, namely ∆ 𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡, ∆𝐼𝑁𝑆𝐸𝐶𝐴𝑅𝑇,𝑗,𝑡, and ∆𝐼𝑁𝑆𝐸𝐶𝐻𝐿𝑂,𝑗,𝑡. ∆𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡 represents the difference in the summations of the use ratios of neonicotinoids (i.e., imidacloprid, dinotefuran, clothianidin, and thiamethoxam) and fipronil between years t and t – 1 in prefecture j to which route i belongs. The reason for this aggregation of insecticides was that the use ratio of each of these insecticides was relatively small (Fig. 3), and it was difficult to detect the effects of each insecticide separately. Moreover, we aggregated neonicotinoids and fipronil together into ∆𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡 because these insecticides all affect nervous system transmission in insects and could be expected to have qualitatively similar adverse effects on S. frequens. ∆𝐼𝑁𝑆𝐸𝐶𝐴𝑅𝑇,𝑗,𝑡 is the difference in use ratios of cartap between years t and t – 1. We did not aggregate cartap into ∆ 𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡 because cartap likely does not have strong adverse effects on S. frequens in paddy fields (Jinguji and Uéda, 2015), although it does also affect nervous system transmission. ∆𝐼𝑁𝑆𝐸𝐶𝐻𝐿𝑂,𝑗,𝑡 is the difference in the use ratios of chlorantraniliprole 11

between years t and t – 1. Chlorantraniliprole may also not have strong adverse effects on dragonflies (Kasai et al., 2016). We did not aggregate chlorantraniliprole into ∆ 𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡 or ∆𝐼𝑁𝑆𝐸𝐶𝐴𝑅𝑇,𝑗,𝑡 because the mode of action of chlorantraniliprole differs from that of these insecticides (it acts on calcium channels in the muscle). We did not include benfuracarb and carbosulfan because their use ratios were small (less than 5%) in all the census years. Our preliminary analyses showed that inclusion of these minor insecticides did not qualitatively affect our analysis (Supplementary code S1). On the basis of our consideration of confounding factors (see section 2.4.2), we included “year” as an explanatory dummy variable (i.e., YEARt). This “year” variable can be interpreted as a composite variable (Supplementary note S4) reflecting 𝑙=𝐿

year-specific properties within ∑𝑙 = 1𝛾𝑙∆𝑂𝐹𝑙,𝑗,𝑖,𝑡 in Eq. 1. As an alternative way of aggregating these use ratios, we also used the following statistical model (model 2):

λ𝑗,𝑖,𝑡 = 𝛼 + 𝛽𝑇𝑂𝑇𝐴𝐿∆𝐼𝑁𝑆𝐸𝑇𝑂𝑇𝐴𝐿,𝑗,𝑡 + 𝛾𝑌𝐸𝐴𝑅𝑡𝑌𝐸𝐴𝑅𝑡 + 𝜀𝑗,𝑖,𝑡,

(Eq. 3)

where ∆𝐼𝑁𝑆𝐸𝑇𝑂𝑇𝐴𝐿,𝑗,𝑡 aggregates ∆𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃,𝑗,𝑡, ∆𝐼𝑁𝑆𝐸𝐶𝐴𝑅𝑇,𝑗,𝑡, and ∆𝐼𝑁𝑆𝐸𝐶𝐻𝐿𝑂,𝑗,𝑡 in Eq. 2; that is, ∆𝐼𝑁𝑆𝐸𝑇𝑂𝑇𝐴𝐿,𝑗,𝑡 represents the difference in the summations of the use ratios of all major insecticides (i.e., imidacloprid, dinotefuran, clothianidin, thiamethoxam, fipronil, cartap and chlorantraniliprole) between years t and t – 1 in prefecture j to which route i belongs. In these statistical models (Eq. 2 and Eq. 3), we assumed that 𝜀𝑗,𝑖,𝑡 in different years t were independent. This assumption means that the temporal autocorrelations over

12

years in these regression models were properly modeled by the included set of explanatory variables in the models. This assumption was statistically tested by the Durbin–Watson test for panel models (Croissant and Millo, 2015).

2.4.4. Statistical analysis Multiple linear regressions were conducted in applying the statistical models 1 (Eq. 2) and model 2 (Eq. 3). In these analyses, we used two metrics for population density (𝑁𝑗,𝑖,𝑡). First, we used maximum population density among the censuses in year t on route i in prefecture j for 𝑁𝑗,𝑖,𝑡. We used maximum population density because the results of the route-based census of S. frequens by car were generally very sensitive to the timing of the census. It was generally difficult to perfectly predict the time at which we would encounter a population of S. frequens on census routes. This sometimes led to low population densities in our survey data, not because the actual density was low but because the timing of the census was mismatched with the dragonfly movements to breeding sites. The use of maximum population density for the year has a clear merit in that it depends little on these low population density data due to mismatched timing. We therefore used maximum population density for 𝑁𝑗,𝑖,𝑡 as a first choice in our analysis. We used median population density for 𝑁𝑗,𝑖,𝑡 as a second choice. The use of median population density provides more robust metrics than the use of maximum population density against chance variabilities due to small sample size. We performed statistical analyses by using the statistical software R version 3.5.3 (R Core Team, 2018). We used the plm package (Croissant and Millo, 2015) for the Durbin–Watson test for panel models. Data and the R code are available in the supplementary materials (Supplementary data and code). 13

3. Results 3.1. Population dynamics of the dragonfly We plotted the population densities of S. frequens on the five routes from 2009 to 2016 (Fig. 4). The yearly fluctuation patterns of population densities varied among routes. Route 1 (in Niigata) had the largest median value in 2011 and 2012, whereas route 4 (in Fukui) had a largest one in 2014. There was a large difference in annual abundance between routes 4 and 5 in the same prefecture. The observed abundance varied widely among census days on each route and in each year, particularly on routes 1 and 4.

3.2. Estimated insecticide use ratios Fipronil and the neonicotinoids clothianidin and dinotefuran were dominant insecticides in the nursery-box treatment of rice seedlings in the early years of the study period (Fig. 3). Their use ratios then gradually decreased, except in Fukui, where fipronil remained dominant. In Toyama, the use ratios of imidacloprid and chlorantraniliprole increased dramatically in 2012 and were still high in 2016. The use ratio of chlorantraniliprole as an alternative to fipronil and neonicotinoids was dominant in Niigata and Ishikawa from about 2012 onward.

3.3 Relationship between insecticides and dragonfly When the maximum population density was used to calculate the population growth rate of S. frequens, the annual increase in total insecticide use ratio was significantly negatively associated with growth rate (Table 2: model 2: coefficient = −0.062, P = 0.045). When we used separate insecticide groups, there was a marginally significant

14

negative association between the sum of neonicotinoids and fipronil and population growth (Table 2: model 1: coefficient = −0.060, P = 0.065). When we used the median population density to calculate the population growth rate, we could not detect any significant association between the annual increase in use ratio of insecticides and population growth (Table 3). The results of the Durbin–Watson test for panel models showed that all our regression models had no temporal autocorrelation (Supplementary code S2).

4. Discussion 4.1. Effect of insecticides In the multiple regression analyses we found a significant negative association between the annual increase in use ratio of all systemic insecticides combined and the population growth of S. frequens. Our multiple regression models took into account a working causal diagram specifying the potential main confounders based on current knowledge. If all confounders are properly controlled for in these analyses (discussed later), the estimated association can be interpreted as a causal effect of the annual increase in insecticides use ratio on the population density of S. frequens. The negative effects of the increase in use of the insecticides on the dragonfly’s population density are consistent with existing knowledge. Nakanishi et al. (2018) reported that increased fipronil use had a negative association with the population growth of S. frequens from 1993 to 2004 in Toyama Prefecture. Many experimental studies have shown high lethal toxicity of fipronil and neonicotinoids to Sympetrum species (reviewed by Nakanishi et al., 2018). Only a few studies have shown an association between increased systemic

15

insecticides use and long-term population declines (e.g., over >1 year) of non-target organisms. Budge et al. (2015) reported a positive association between use of the neonicotinoid imidacloprid and honey bee colony loss rates. Woodcock et al. (2016) showed a positive association between neonicotinoid use and population extinction rates of wild bees. Forister et al. (2016) reported a negative association between neonicotinoid use and annual variation in the species richness of butterflies. Our study is among the first to demonstrate that such insecticides can negatively affect the long-term population dynamics of wild populations. Contrary to our original expectation, the negative effect of the annual increase in use ratio of all systemic insecticides combined (∆𝐼𝑁𝑆𝐸𝑇𝑂𝑇𝐴𝐿 in model 2 in Table 2) was almost the same as that of neonicotinoids and fipronil (∆𝐼𝑁𝑆𝐸𝑁𝑁𝐹𝑃 in model 1 in Table 2), even though the variable of increase in use ratio of all insecticides included insecticides with weak adverse effects on dragonflies (i.e., cartap and chlorantraniliprole). The first possible interpretation is that cartap and chlorantraniliprole had a certain level of negative impact on the dragonfly population in the field. The negative coefficients of these insecticides (model 1 in Table 2) may support this interpretation, although the results of mesocosm experiments suggest that these adverse effects are limited (Jinguji and Uéda, 2015; Kasai et al., 2016). To examine this disparity between such experimental results and our statistical analyses of the field data, we need more information about effects on the dragonfly in real rice fields under the application of these insecticides.

4.2. Potential confounding effects of other agronomic practices In conventional management of rice cultivation in Japan, most farmers use systemic 16

insecticides on their nursery boxes. In contrast, in environmentally friendly rice cultivation (e.g., organic farming), systemic insecticides tend not to be applied to nursery boxes. Therefore, the use or non-use of systemic insecticides can be closely associated with the type of cultivation regime. In addition, many agronomic practices other than insecticide use (e.g., the use of midsummer drainage or the method of tillage) differ between conventional and environmentally friendly regimes. The variable of increased annual use ratio of all systemic insecticides therefore can also work as a surrogate variable for unmeasured variables associated with these agronomic practices. A substantial decrease in the emergence percentage of S. frequens in a conventional rice cultivation regime compared with that in an environmentally friendly regime has been reported from a field survey, although individual agronomic factors were not separated (Baba et al., 2019). Because we could not include agronomic factors related to conventional regimes (e.g., the use of midsummer drainage or the method of tillage) in the models, the negative effects of these factors might have partly been included in—that is, might have biased—the effects of annual increase in total insecticide use ratio. This means that the effects of insecticides can reflect factors other than the toxic effects of the insecticides themselves. This could be another possible interpretation of the result that the inclusion of insecticides with potentially low adverse effects did not decrease the negative effect of the sum total of insecticides on the dragonfly (Table 2). No statistical data about other agronomic practices in the study area during the census period were available to us. Conducting field surveys in the future to identify the effects of agronomic factors related to the different types of cultivation regimes will help us to evaluate the effects of these unmeasured potential confounders.

17

4.3. Limitations of field census data The results of our regression analyses differed between the two cases of maximum (Table 2) and median (Table 3) population density. One possible reason for this difference was that median population densities were likely to have been underestimated by low population density data owing to mismatches between the timing of high population density and that of the census. These underestimations may have obscured the effects of the insecticides. On the other hand, maximum population densities are likely to be heavily influenced by the number of censuses performed. Our data had inconsistency in the number of censuses on each route and in each year (Table 1). This inconsistency could have increased the variance of our data, especially when we used maximum densities, although the data set showed a relatively good fit in our regressions (i.e., the residuals showed no outlier and were roughly flat over explanatory variables; see residual plots in Supplementary note S5). This greater variance might also have affected our estimates of the effects of the insecticides, although they were not expected to cause systematic bias to the estimates because they would have been independent of insecticide use.

4.4. Future work One of the limitations of our analyses is that data on the amounts of insecticides shipped were available only by prefecture, whereas the dragonfly density data were recorded on routes in local areas of each prefecture. We set two census routes (routes 4 and 5) in Fukui Prefecture, and the population density of the dragonfly differed greatly between these two routes in some years. This local variation in dragonfly population density

18

could have resulted from local-level differences in the use ratios of insecticides. Interviews with members of local agricultural cooperatives suggested that cartap was used frequently in the rice fields around route 4, whereas fipronil was used frequently around route 5. Quantitative data on the insecticide use around each route were not available. As future research, we are planning to interview local farmers to obtain information on insecticide use at a more local scale. We also plan to conduct field surveys of dragonfly population density within the corresponding areas to evaluate the causal effect of insecticides at a finer scale. In addition, we will conduct a survey (e.g., a comparison of the vital rates of dragonflies under different management systems) in rice paddies to evaluate the impacts of agronomic factors other than systemic insecticide application to nursery boxes. The results of this study, combined with the results of these future field studies, will provide a scientific basis for effectively restoring and conserving dragonfly populations in paddy fields.

Acknowledgments We thank Dr. Noriko Kusamitsu (The Japan Association of Rural Solutions for Environmental Conservation and Resource Recycling) for her help with the field survey. We are grateful to Dr. Takashi Nagai (Institute for Agro-Environmental Sciences, NARO), Dr. Ryo Futahashi and Dr. Masashi Kamo (National Institute of Advanced Industrial Science and Technology), and Dr. Koichi Goka (National Institute for Environmental Studies) for their valuable technical advice. We also thank Dr. Keiichi Fukaya (National Institute for Environmental Studies) for his comments on the statistical analyses and Dr. Yoshitaka Imaizumi, Ayumi Matsubara, and Fumie Tanaka

19

(National Institute for Environmental Studies) for providing the data on insecticides, as well as technical support. This study was supported by the Environment Research and Technology Development Fund (4-1701) of the Environmental Restoration and Conservation Agency of Japan.

References Baba, Y.G., Kusumoto, Y., Tanaka, K., 2019. Positive effect of environmentally friendly farming on paddy field odonate assemblages at a small landscape scale. J. Insect Conserv. 1–8. https://doi.org/10.1007/s10841-019-00132-2 Bass, C., Denholm, I., Williamson, M.S., Nauen, R., 2015. The global status of insect resistance to neonicotinoid insecticides. Pestic. Biochem. Phys. 121, 78–87. https://doi.org/10.1016/j.pestbp.2015.04.004 Budge, G.E., Garthwaite, D., Crowe, A., Boatman, N.D., Delaplane, K.S., Brown, M.A., Thygesen, H.H., Pietravalle, S., 2015. Evidence for pollinator cost and farming benefits of neonicotinoid seed coatings on oilseed rape. Sci. Rep. 5, 12574. https://doi.org/10.1038/srep12574 Croissant, Y., Millo, G., 2015. Panel Data Econometrics in R : The plm Package. J. Stat. Softw. 27, 1–43. https://doi.org/10.18637/jss.v027.i02 Forister, M.L., Cousens, B., Harrison, J.G., Anderson, K., Thorne, J.H., Waetjen, D., Nice, C.C., De Parsia, M., Hladik, M.L., Meese, R., van Vliet, H., Shapiro, A.M., 2016. Increasing neonicotinoid use and the declining butterfly fauna of lowland California. Biol. Lett. 12, 20160475. https://doi.org/10.1098/rsbl.2016.0475 Furlan, L., Pozzebon, A., Duso, C., Simon-Delso, N., Sánchez-Bayo, F., Marchand, P.A., Codato, F., Bijleveld van Lexmond, M., Bonmatin, J.-M., 2018. An update of

20

the Worldwide Integrated Assessment (WIA) on systemic insecticides. Part 3: alternatives to systemic insecticides. Environ. Sci. Pollut. Res. 1–23. https://doi.org/10.1007/s11356-017-1052-5 Futahashi, R., 2012. Recent decline of red dragonflies in Toyama Prefecture. Nat. Insects 47, 10–15. Hernán, M.A., Robins, J.M., 2019. Causal Inference. Chapman & Hall/CRC, Boca Raton. Inoue, K., Tani, K., 2010. All about red dragonflies. Tombow Publishing, Osaka. Jinguji, H., Thuyet, D.Q., Uéda, T., Watanabe, H., 2013. Effect of imidacloprid and fipronil pesticide application on Sympetrum infuscatum (Libellulidae: Odonata) larvae and adults. Paddy Water Environ. 11, 277–284. https://doi.org/10.1007/s10333-012-0317-3 Jinguji, H., Uéda, T., 2015. Can the use of more selective insecticides promote the conservation of Sympetrum frequens in Japanese rice paddy fields (Odonata: Libellulidae)? Odonatologica 44, 63–80. Jinguji, H., Uéda, T., Goka, K., Hidaka, K., Matsura, T., 2009. Effects of imidacroprid and fipronil insecticide application on the larvae and adults of Sympetrum frequens (Libellulidae:Odonata). Trans. Japanese Soc. Irrig. Drain. Rural Eng. 77, 35–41. Kadoya, T., Suda, S., Washitani, I., 2009. Dragonfly crisis in Japan: A likely consequence of recent agricultural habitat degradation. Biol. Conserv. 142, 1899– 1905. https://doi.org/10.1016/J.BIOCON.2009.02.033 Kasai, A., Hayashi, T.I., Ohnishi, H., Suzuki, K., Hayasaka, D., Goka, K., 2016. Fipronil application on rice paddy fields reduces densities of common skimmer and scarlet skimmer. Sci. Rep. 6, 1–10. https://doi.org/10.1038/srep23055

21

Katayama, N., Baba, Y.G., Kusumoto, Y., Tanaka, K., 2015. A review of post-war changes in rice farming and biodiversity in Japan. Agric. Syst. 132, 73–84. https://doi.org/10.1016/j.agsy.2014.09.001 Katayama, N., Osada, Y., Mashiko, M., Baba, Y.G., Tanaka, K., Kusumoto, Y., Okubo, S., Ikeda, H., Natuhara, Y., 2019. Organic farming and associated management practices benefit multiple wildlife taxa: A large‐scale field study in rice paddy landscapes. J. Appl. Ecol. 1365-2664.13446. https://doi.org/10.1111/1365-2664.13446 Lawler, S.P., 2001. Rice fields as temporary wetlands: a review. Isr. J. Zool. 47, 513– 528. https://doi.org/10.1560/X7K3-9JG8-MH2J-XGX1 Ministry of Agriculture, Forestry and Fisheries (Japan), 2018. Statistics of Crops [WWW Document]. URL http://www.maff.go.jp/j/tokei/kouhyou/sakumotu/index.html (accessed 12.26.18). Morgan, S.L., Winship, C., 2014. Counterfactuals and Causal Inference : Methods And Principles For Social Research, 2nd ed. New York. Nakanishi, K., Yokomizo, H., Hayashi, T.I., 2018. Were the sharp declines of dragonfly populations in the 1990s in Japan caused by fipronil and imidacloprid? An analysis of Hill’s causality for the case of Sympetrum frequens. Environ. Sci. Pollut. Res. 25, 35352–35364. https://doi.org/10.1007/s11356-018-3440-x Natuhara, Y., 2013. Ecosystem services by paddy fields as substitutes of natural wetlands in Japan. Ecol. Eng. 56, 97–106. https://doi.org/10.1016/j.ecoleng.2012.04.026 Pearl, J., 2009. Causality: Models, Reasoning, and Inference, 2nd ed. Cambridge University Press, New York.

22

Pearl, J., Mackenzie, D., 2018. The Book of Why: The New Science of Cause and Effect. Penguin Books, London. Pisa, L., Goulson, D., Yang, E.-C., Gibbons, D., Sánchez-Bayo, F., Mitchell, E., Aebi, A., Sluijs, J. van der, MacQuarrie, C.J.K., Giorio, C., Yim Long, E., McField, M., Van Lexmond, M.B., Bonmatin, J.-M., 2017. An update of the Worldwide Integrated Assessment (WIA) on systemic insecticides. Part 2: impacts on organisms and ecosystems. Environ. Sci. Pollut. Res. https://doi.org/DOI 10.1007/s11356-017-0341-3 Pisa, L.W., Amaral-Rogers, V., Belzunces, L.P., Bonmatin, J.M., Downs, C.A., Goulson, D., Kreutzweiser, D.P., Krupke, C., Liess, M., Mcfield, M., Morrissey, C.A., Noome, D.A., Settele, J., Simon-Delso, N., Stark, J.D., Van Der Sluijs, J.P., Van Dyck, H., Wiemers, M., 2015. Effects of neonicotinoids and fipronil on non-target invertebrates. Environ. Sci. Pollut. Res. 22, 68–102. https://doi.org/10.1007/s11356-014-3471-x R Core Team, 2018. R: A language and environment for statistical computing [WWW Document]. R Found. Stat. Comput. Vienna, Austria. URL https://www.r-project.org/ (accessed 9.18.18). Sánchez-Bayo, F., Goka, K., 2006. Ecological effects of the insecticide imidacloprid and a pollutant from antidandruff shampoo in experimental rice fields. Environ. Toxicol. Chem. 25, 1677–1687. https://doi.org/10.1897/05-404R.1 Simon-Delso, N., Amaral-Rogers, V., Belzunces, L.P., Bonmatin, J.M., Chagnon, M., Downs, C., Furlan, L., Gibbons, D.W., Giorio, C., Girolami, V., Goulson, D., Kreutzweiser, D.P., Krupke, C.H., Liess, M., Long, E., Mcfield, M., Mineau, P., Mitchell, E.A., Morrissey, C.A., Noome, D.A., Pisa, L., Settele, J., Stark, J.D.,

23

Tapparo, A., Van Dyck, H., Van Praagh, J., Van Der Sluijs, J.P., Whitehorn, P.R., Wiemers, M., 2015. Systemic insecticides (Neonicotinoids and fipronil): Trends, uses, mode of action and metabolites. Environ. Sci. Pollut. Res. 22, 5–34. https://doi.org/10.1007/s11356-014-3470-y Sugimura, M., Ishida, S., Kojima, K., Ishida, K., Aoki, T., 1999. Dragonflies of the Japanese Archipelago in Color. Hokkaido University Press, Sapporo. Thuyet, D.Q., Watanabe, H., Motobayashi, T., Ok, J., 2013. Behavior of nursery-box-applied fipronil and its sulfone metabolite in rice paddy fields. Agric. Ecosyst. Environ. 179, 69–77. https://doi.org/10.1016/j.agee.2013.07.012 Uéda, T., 1988. Diversity in life history of the dragonfly Sympetrum frequens (Odonata: Insecta). Bull. Ishikawa Agric. Coll. 18, 98–110. Uéda, T., 1998. Dragonfly communities in paddy fields, in: Ezaki, Y., Tanaka, T. (Eds.), Conservation of Wetland Environments: A View from Biological Communities. Asakura Book Co., Tokyo, pp. 93–110. Uéda, T., Jinguji, H., 2013. The ecological impact of the insecticides fipronil and imidacloprid on Sympetrum frequens in Japan. Tombo 55, 1–12. Van Der Sluijs, J.P., Amaral-Rogers, V., Belzunces, L.P., Bijleveld Van Lexmond, M.F., Bonmatin, J.M., Chagnon, M., Downs, C.A., Furlan, L., Gibbons, D.W., Giorio, C., Girolami, V., Goulson, D., Kreutzweiser, D.P., Krupke, C., Liess, M., Long, E., Mcfield, M., Mineau, P., Mitchell, E.A., Morrissey, C.A., Noome, D.A., Pisa, L., Settele, J., Simon-Delso, N., Stark, J.D., Tapparo, A., Van Dyck, H., Van Praagh, J., Whitehorn, P.R., Wiemers, M., 2015. Conclusions of the worldwide integrated assessment on the risks of neonicotinoids and fipronil to biodiversity and ecosystem functioning. Environ. Sci. Pollut. Res. 22, 148–154.

24

https://doi.org/10.1007/s11356-014-3229-5 Woodcock, B.A., Isaac, N.J.B., Bullock, J.M., Roy, D.B., Garthwaite, D.G., Crowe, A., Pywell, R.F., 2016. Impacts of neonicotinoid use on long-term population changes in wild bees in England. Nat. Commun. 7, 12459. https://doi.org/10.1038/ncomms12459 Yachi, S., Nagai, T., Inao, K., 2016. Development of a simple method for estimating the usage of paddy insecticides with each application method. Japanese J. Pestic. Sci. 41, 1–10. https://doi.org/10.1584/jpestics.W15-31 Yachi, S., Nagai, T., Inao, K., 2017. Analysis of region-specific predicted environmental concentration of paddy pesticides at 350 river flow monitoring sites. Japanese J. Pestic. Sci. 42, 1–9. https://doi.org/10.1584/jpestics.W17-05

25

Figure captions

Fig. 1. Study area and examples of each daily census route in the four prefectures

Fig. 2. Working causal diagram used to distinguish factors to be controlled in the causal analyses between insecticide use and population decline of Sympetrum frequens

Fig. 3. Estimated use ratios of insecticides in nursery-box treatment of rice seedlings in the four prefectures

Fig. 4. Annual median numbers of tandem pairs of Sympetrum frequens along five routes in the four prefectures. Lower and upper bars indicate minimum and maximum values, respectively, in each year.

Highlights • Dragonfly abundance in Japan’s agricultural landscapes has decreased. • We monitored the population density of the dragonfly Sympetrum frequens. • We estimated systemic insecticide use and analyzed its effect on the dragonfly. • The results suggest a negative effect of systemic insecticide use on population size. • Some agronomic factors need to be further examined as potential confounders. Table 1. Numbers of censuses along each route in four prefectures

26

Route 1

Route 2

Route 3

Route 4

Route 5

(Niigata)

(Toyama)

(Ishikawa) (Fukui)

(Fukui)

2009

3

0

1

1

0

2010

2

2

3

3

3

2011

2

4

5

4

3

2012

2

3

4

5

5

2013

1

2

5

3

4

2014

1

3

4

4

3

2015

1

5

6

3

1

2016

0

3

9

4

1

Total

12

22

37

27

20

Year

Table 2. Summary of results of multiple regression analyses using the population growth rate of Sympetrum frequens calculated with maximum population density as a response variable Model 1

Model 2

Explanatory variable

Estimated

Estimated SE

P value

coefficient

SE

P value

coefficient

∆INSENNFP

–0.060

0.031

0.065







∆INSECART

–0.144

0.157

0.369







∆INSECHLO

–0.073

0.062

0.253







∆INSETOTAL







-0.062

0.029

0.045*

27

YEAR2011

1.259

0.717

0.093

1.074

0.599

0.086

YEAR2012

0.732

0.656

0.277

0.593

0.588

0.323

YEAR2013

-0.316

0.702

0.657

-0.524

0.589

0.383

YEAR2014

0.669

0.665

0.326

0.747

0.591

0.218

YEAR2015

–0.073

0.653

0.912

-0.187

0.596

0.756

YEAR2016

–0.293

0.671

0.667

-0.411

0.615

0.510

NNFP: sum of neonicotinoids and fipronil; CART: cartap; CHLO: chlorantraniliprole; TOTAL: sum of all insecticides *: P < 0.05 Durbin–Watson test for panel models: Model 1, DW = 2.548, P = 0.739; Model 2, DW = 2.637, P = 0.812

Table 3. Summary of results of multiple regression analyses using the population growth rate of Sympetrum frequens calculated with median population density as a response variable Model 1

Model 2

Explanatory variable

Estimated

Estimated SE

P value

coefficient

SE

P value

coefficient

∆INSENNFP

–0.006

0.036

0.862







∆INSECART

–0.130

0.182

0.482







∆INSECHLO

0.010

0.072

0.888







∆INSETOTAL







-0.012

0.034

0.731

28

YEAR2011

1.842

0.831

0.037*

1.793

0.697

0.017*

YEAR2012

0.981

0.761

0.211

0.906

0.685

0.198

YEAR2013

0.456

0.814

0.581

0.254

0.686

0.714

YEAR2014

0.582

0.771

0.458

0.794

0.687

0.259

YEAR2015

0.816

0.757

0.293

0.652

0.693

0.356

YEAR2016

0.575

0.778

0.468

0.515

0.716

0.479

NNFP: sum of neonicotinoids and fipronil; CART: cartap; CHLO: chlorantraniliprole; TOTAL: sum of all insecticides *: P < 0.05 Durbin–Watson test for panel models: Model 1, DW = 2.210, P = 0.349; Model 2, DW = 2.213, P = 0.328

29

N

0

25

Niigata Route 1

36°50′ N

Toyama Ishikawa Route 2

Fukui

Route 3

36°00′ N

Route 5

136°00′ E

Route 4

137°00′ E

138°00′ E

50 km

Site Meteorological factors

Year

Agronomic practices

Insecticides

Dragonfly

a. Niigata

b. Toyama 80

Use ratio (%)

Use ratio (%)

80 60 40 20 0

60 40 20 0

2009 2010 2011 2012 2013 2014 2015 2016

2009 2010 2011 2012 2013 2014 2015 2016

Year

Year

c. Ishikawa

d. Fukui 80

Use ratio (%)

Use ratio (%)

80 60 40 20 0

60 40 20 0

2009 2010 2011 2012 2013 2014 2015 2016

2009 2010 2011 2012 2013 2014 2015 2016

Year

Year

Fipronil

Imidacloprid

Dinotefuran

Clothianidin

Chlorantraniliprole

Cartap

Benfuracarb

Carbosulfan

Thiamethoxam

Route # (Prefecture) 1 (Niigata) 2 (Toyama) 3 (Ishikawa) 4 (Fukui) 5 (Fukui)

Pairs/km

60

40

20

0 2009

2010

2011

2012

2013

Year

2014

2015

2016