Climate warming promotes pesticide resistance through expanding overwintering range of a global pest

Insect preparation

We collected 200–300 larvae and pupae of the diamondback moth from cabbage and cauliflower fields in Wuhan, Beijing, and Shenyang in late September from 2008 to 2012. We mixed all collections into one stock colony because of no geographic differentiation in this species from these sites⁵⁹. We reared all individuals on cabbage leaves spread evenly across five screen cages (35 × 35 × 15 cm) in growth chambers at constant temperature (25 ± 1 °C) with 15-h light:9-h dark photoperiod, and relative humidity set at 60 ± 10%. We moved any new pupae to new screen cages for adult emergence. Emerged adults were fed with 10% honey solution in cotton balls. To collect eggs, we dipped four small pieces of laboratory film (7 × 5 cm) in fresh cabbage juice for 3–5 s and hung the treated film pieces on the top of each screen cage. To further enlarge the population for our experiments, we reared these insects in the artificial diet in plastic boxes at 25 ± 1 °C. We transferred 200 eggs to the surface of 120 g artificial diet (Southland Products Incorporated, USA) in each plastic box (10 × 10 × 9 cm). The hatched larvae dropped to the surface of the artificial diet and fed on it. Once individuals developed into 3rd or 4th instar larvae, pupae or adults, they were exposed to 10 °C for 24 h (to simulate gradually reduced temperatures in late autumn and allow a thermal acclimation) just before they were placed to the low-temperature regimes for overwintering tests. Overall, we obtained >7000 larvae and >8000 pupae for the laboratory experiment, and >8000 larvae, >8000 pupae, >4000 adults for the field experiment. We have compared the life history traits of insects reared on the artificial diet with natural host plants (cabbage leaves), they performed similarly (Peng and Li, unpublished data).

Laboratory experiment of winter survival

Site selection

To identify what factor determines winter survival under different winter thermal conditions, we conducted a laboratory experiment that simulated temperature regimes of 10 selected sites across a latitudinal gradient in China (Fig. 1, Supplementary Table 1) at which this species is known to occur and damage cruciferous crops during the growing season.

Temperature treatment

To simulate the winter temperatures in the 10 geographically distinct sites (Fig. 1a, Supplementary Table 1), we collected daily mean temperatures during winter (November to next April of 1966–2010) at each site from China Meteorological Data Service Centre (http://data.cma.cn/en). Then, we fitted a polynomial model to the temporal changes of winter daily mean temperatures for each site (Fig. 1b). To simplify the logistics of temperature control procedures, we set all temperature regimes in combinations of linear decline, horizontal maintenance and linear increase to mimic the polynomial changes of winter temperatures in the 10 sites, and adjusted temperature every 10 days as needed (Fig. 1c). We controlled the winter temperature changes of the 10 sites with climate chambers (RXZ-280B, Jiangnan Ltd., Ningbo, China) and refrigerators (Royalstar BCD-246GER) according to curves in Fig. 1c.

Experimental protocols

We conducted a winter survival experiment with 10 low-temperature regimes (Fig. 1c). We exposed 6050 larvae and 7150 pupae to 10 low-temperature regimes according to the experiment design (Fig. 1c). Then we sampled 55 larvae and 65 pupae every 10 days from each temperature regime resulting in 11 sampling points. Sampled larvae were placed at 25 °C for 1.5 days to observe the survival based on if their body kept fresh green⁶⁰ and appendage moved after touching with a brush³⁴. The pupae were placed at 25 °C, RH 70–80% and photoperiod of 16 L:8D for emergence to determine the survival (emergence rate). These samples were not returned to the temperature treatments. Thus, no individual was measured more than once and each sample interval represents an independent observation.

Field survival experiments across 12 geographic sites

To verify the cold survivals from the laboratory simulation and identify the best predictor under natural conditions, we conducted field experiments to explore the winter survival for multiple years at various geographic sites in China (Fig. 1a, Supplementary Table 1). The diamondback moth overwinters either in remaining cabbage plants or in fallen leaves (post-harvest conditions) in regions without standing cabbage crops in the winter. We tested the winter survival of larvae, pupae and adults in the caged cabbage plants or post-harvest conditions in fallen leaves on the soil surface at each site for 3–4 months. We transferred 30 larvae, 30 pupae or 30 adults from our stock rearing to a cabbage plant in the field. Then each plant was covered with a screen cage to avoid disturbance and contain focal individuals (see Supplementary Fig. 2). We set 6–8 cages for larvae, pupae and adults, respectively, in a field in November or early December. After an exposure of 1, 2, 3 and 4 months, we collected 2 cages of larvae, 2 cages of pupae and 2 cages of adults from the field at each sampling point and kept individuals in the laboratory (25 ± 1 °C, RH 65–75%, L:D = 16:8 h) for two days. We checked the survival status of the larvae based on the change in body coloration (i.e. if the larval body kept fresh green colour)⁶⁰, pupal survival based on whether adults could emerge from the pupae, and adult survival based on if their appendage moved after touching with a brush.

To simulate the field microenvironment of post-harvest conditions in winter, we filled half of a glass jar (diameter = 5.5 cm, height = 14 cm) with moist soil. Then, we transferred 30 larvae or 30 pupae to the soil surface, covered the insects with leaves, and then covered the glass jar with a nylon net (see Supplementary Fig. 2). We buried 6–8 jars for larvae and pupae, respectively, and kept the top of the jar at ground surface level at each site in November and early December. Because almost all adults died in few days within the jar, we did not test the adult survival in post-harvest conditions. After an exposure of 1, 2, 3 and 4 months, we took 2 jars of larvae and 2 jars of pupae per sampling period from the field and placed them in the laboratory with 25 ± 1 °C, RH 65–75%, L:D = 16:8 h for 2 days. The survival status of the larvae and pupae was checked with the same procedures as the overwintering tests on caged cabbage plants. Note that as in the standing plant experiment, no individual was tested more than once assuming that each observation is independent at the replicate level.

Modelling and predicting winter survival

Model development

Our goal was to identify key metrics that best predict the winter survival of the diamondback moth across a climatic gradient. To achieve this goal, we took several steps. First, we fit a set of predictive models to the laboratory experiments to identify which metric and model best describes survival under controlled conditions. We focused on three alternative predictors: the lowest daily mean temperature (MinDT_mean), mean temperature (DT_mean) combined with exposure days, and low-temperature degree-days (LTDD). The MinDT_mean model assumes that survival can simply be predicted as a function of the lowest temperature an individual experienced during its exposure time. The DT_mean model assumes that survival depends on both the average temperature individuals experience below the cold threshold for survival (11 °C)³² and exposure duration (note that exposure time varied systematically in 10-day increments). Finally, the LTDD model predicts survival depending on coldness below the cold threshold. We calculated LTDD by summing up negative deviations of daily mean temperatures from the cold threshold (11.0 °C) during each exposure period for each simulated geographic site (Fig. 1c). To detect potential relationships, we fit each model using three different functions, i.e. linear, exponential and sigmoid models to describe the survival probability (Supplementary Table 2, Fig. 1). We estimated parameters of models in SigmaStat 3.5 and compared model fit using R² and AIC values (see detailed models in Supplementary Table 2).

Field validation of survival models

To validate winter survival models derived from the laboratory (see models in Supplementary Table 2) for complex and variable field conditions, we compared model predictions to observed survivals in field experiments across 12 different geographic sites over 5 years (Fig. 1a, Supplementary Table 1). To make the connection, we first collected daily mean temperatures recorded at the nearest weather stations to our field sites from China Meteorological Data Service Centre. We then calculated MinDT_mean, DT_mean and exposure days, and LTDD for each site for each treated period and input these values into these laboratory models to predict winter survival. Note, that because the coefficients were calculated from the laboratory experiment, predictions are completely independent of survival observed under field conditions. During the model validation, we excluded the field data of south China, e.g. Guangzhou, Changsha and Wuhan where the warmer temperatures allowed moths to continue their regular life cycle during the whole winter, resulting in unrealistic winter survival. We also excluded replicates in which glass jars were filled with water and destroyed the tested insects. We used linear regression to compare predicted survival with field observations. The validity of each model was evaluated based on the variance explained, slopes of linear regressions and prediction bias (i.e. deviation from unity slope). Finally, we selected the exponential model driven by LTDD as the model to predict the global distribution of winter survival due to its lowest AIC value (Supplementary Table 2) and the least bias (Supplementary Table 3, Fig. 2) among all models.

Global prediction of overwintering range shift

To extrapolate our winter survival predictions to a global scale under present and future climate conditions, we downloaded global historical daily mean temperature data for 50 years (1967–2016) from Berkeley Earth (1° × 1° grid, http://berkeleyearth.org/data/). We added 1, 2, 3, 4, 5 and 6 °C to mean temperatures of 2012–2016, respectively, to represent the different future warming scenarios³⁷. Then, we calculated the annual LTDD in the northern hemisphere with Eq. (1) and in the southern hemisphere with Eq. (2). For x_i,j < x₀,

where x_i,j is the daily mean temperature at each grid, i is the year, j is the Julian day, and x₀ = 11.0 °C is the cold threshold for the diamondback moth survival. If x_i,j > x₀, we excluded the x_i,j for the calculation LTDD. For Eq. (1), we started the calculation of LTDD from July 1st (Julian date 182), ended on June 30th of next year (Julian date 181) to cover the whole low-temperature season in the northern hemisphere cross the calendar year. We used LTDD for every year during past conditions to our validated survival model (LTDD-dependent exponential model) and further calculated the expected corresponding yearly winter survival and 5-year mean survival. Since the diamondback moth only feeds on Brassicaceae plants⁶¹, we incorporated host availability to refine the pest distributions. We retrieved Brassicaceae occurrence data during 1967–2016 (3,720,971 records) from the Global Biodiversity Information Facility (GBIF) database (www.gbif.org), and excluded unknown and duplicate records; 919,808 records were retained to model the global distribution of host plants. We used a dataset of eight selected bioclimatic variables as described in a previous Brassicaceae biogeographic study⁶², including isothermality (bio3), temperature seasonality (bio4), min temperature of coldest month (bio6), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), precipitation seasonality (bio15), precipitation of warmest quarter (bio18), precipitation of coldest quarter (bio19) from Worldclim dataset⁶³ (http://worldclim.org). We ran the species distribution model using the Maxent algorithm in R package dismo⁶⁴. Model outputs were presented in grid ranks of host plant presence probability from 0 (unsuitable) to 1 (most suitable). Based on the known distribution of Brassicaceae, we only included grid cells with Brassicaceae presence probability ≥0.3 for our final survival and distribution analysis to ensure the presence of the host plant and mapped them with Arcmap 10.2 (Environmental Systems Research Institute) (see Fig. 3a, e). To show spatial-temporal changes in the geographic distribution of winter survival, we quantified the historical change (expansions or contractions) in the overwintering range based on the total numbers of grids for each year between 1967 and 2016 relative to the baseline area in 1967 (see Fig. 3f) and further calculated average changes of every 5 years (see Fig. 3b). We selected sites in the overwintering marginal belt (with winter survival between 1 and 5%) in the baseline year (1967), calculated the annual LTDD of these sites from 1967 to 2016, and built the linear trend of annual LTDD for years 1967–2016 (see Fig. 3g). We predicted distribution changes for future scenarios (added 1, 2, 3, 4, 5 and 6 °C to the current mean temperatures of 2012–2016) relative to the baseline area of 1967–1971 (see Fig. 3c, d).

Meta-analysis linking pesticide resistance to overwintering type

Data preparation: literature search and selection criteria

We performed a comprehensive literature survey to collect data on pesticide resistance of the diamondback moth worldwide. We searched for publications in databases of ISI Web of Science, Scopus and China National Knowledge Infrastructure (CNKI) using keywords “pesticide resistance” in combination with “diamondback moth” or “Plutella xylostella” and expanded references in the selected papers. We reviewed titles, abstracts and in many cases the full articles for relevance and agreement with our inclusion criteria. Studies were included if they (1) monitored the pesticide resistance of field populations, (2) used the leaf dip bioassay method to test pesticide resistance which is the most commonly used method recommended by Insecticide Resistance Action Committee (IRAC, http://www.irac-online.org); (3) provided resistance ratio of field populations. Resistance ratio (abbreviated as RR) is the magnitude of pesticide resistance and is commonly calculated by dividing the median lethal concentration (LC50) of a tested field population by LC50 of the susceptible population (without exposure to pesticide). The LC50 is commonly estimated from a concentration-mortality curve of a given pesticide. The preliminary literature search resulted in 2151 studies out of which 62 matched these criteria. A PRISMA diagram describing details of our literature search is available in Supplementary Fig. 4.

Data preparation: data extraction

We extracted data from each selected publication, the names of pesticides, sampling locations and years of field populations, number of tested individuals in a bioassay, resistance ratio of field populations (RR), LC50 of field populations (LC50_field) and susceptible populations (LC50_susceptible), and 95% confidence intervals (CIs) of LC50_field and LC50_susceptible. Some studies generated results from multiple types of pesticides with the same field population, each of which was considered as a different entry. Finally, we gathered 1806 entries for pesticide resistance of field populations of the diamondback moth.

Data preparation: calculation of the weighted effect size

We conducted a meta-analysis to test if pesticide resistance levels vary across different types of overwintering sites. To account for differences in sample sizes and variances in resistance ratios across studies, we calculated the corrected (weighted) resistance ratio for each study following the method in Hedges et al.⁶⁵. We calculated the logarithm of resistance ratio (logRR) to present the effect size for each entry and further calculated the weighted effect size (wlogRR) by

where w is the weighting factor of each entry, with w = 1/sqrt(V_logRR)⁶⁶. To consider the contribution from both field and susceptible population, the pooled variance V_logRR was calculated as follows⁶⁵:

where LC50_field and LC50_susceptible, SE_field and SE_susceptible, n_field and n_susceptible, are LC50, the standard error of LC50 and sample size for field population and susceptible population, respectively. SE_field and SE_susceptible can be calculated from their own confidence intervals (95% CI)⁶⁷:

where CI_{upper limit} is the upper limit and CI_{lower limit} is the lower limit of the 95% CI for LC50.

We used the prognostic method⁶⁸ to estimate V_logRR for entries that miss either 95% CI or LC50 based on the average V_logRR of the other complete entries.

Data preparation: potential moderator variables

Several factors could influence pesticide resistance besides overwintering temperatures. The effective temperature degree-days (ETDD) may change the annual number of generations, the intensity of pesticide application, and thus the selection stress⁴⁷, e.g. between 7.4 and 33 °C for the diamondback moth³². In addition, the variety of pesticides used in a study may also affect the resistance levels through their mode of actions (the lethal mechanism) and cross-resistance^69,70. To account for these potentially confounding factors, we collected the mode of action for each variety of pesticides from IRAC, and calculated LTDD, ETDD and overwintering type for each of the 1806 original records. We collected data for daily mean temperatures for each site from Berkeley Earth. For each location, we calculated the mean annual LTDD, ETDD and winter survival average across the 5 years before the sample. We split the sampling sites into three types based on predicted winter survival of the diamondback moth: (1) the permanent (overwintering) sites: locations with the mean winter survivals ≥5%, (2) marginal sites: locations with the mean winter survivals 1–5%, (3) transient (non-overwintering) sites: locations with the mean winter survivals <1%. We mapped all three category sites in a world map with Arcmap 10.2 (see Fig. 4d).

Data analysis: impact of overwintering types

To isolate the effects of overwintering types (permanent, marginal, transient) on resistance from ETDD and pesticide variety, we fit a linear mixed model using lmer function from R package lme4. We used this model with the weighted effect size (wlogRR) as the response variable, overwintering type and pesticide variety as fixed factors, ETDD as a covariate, and combination of sampling location and year as a random factor. We used a Wald Chi-square test to evaluate the significance of each fixed factor, covariate and their interactions using Anova function from R package car. R scripts are provided as Supplementary Software.

To show the resistance level in clearer biological meanings, the resistance ratios (RR) were classified into five levels⁷¹, (1) susceptible: RR ≤ 1, (2) tolerance/low resistance: 1 < RR < 10, (3) moderate resistance:10 ≤ RR < 100, (4) high resistance: 100 ≤ RR < 1000, (5) extremely high resistance: RR ≥ 1000. We calculated the incidence frequencies of these five resistance levels, respectively, and showed the frequencies in a pie chart for the three overwintering types (see Fig. 4a). Furthermore, we calculated the weighted mean resistance ratio⁶⁶ for the permanent, marginal and transient overwintering sites, respectively.

Data analysis: distribution of high-level resistance

Pesticide resistance varies markedly from region to region due to different histories of insecticide application on the farm level, which could obscure the differences in statistics between different regions. However, practitioners of pest control are especially concerned about high levels of pesticide resistance in their regions. Thus, we built quantile regression models driven by LTDD to predict the global distribution of resistance levels (log₁₀-transformed resistance ratio) in R package quantreg (see Supplementary Fig. 3), especially we focused on the significant 0.85 quantile model (see Fig. 4b) for the distribution of the top 15% resistance levels in China and the world in 2012–2016. Areas with no host plants have no resistance problems, so we excluded the grids where Brassicaceae were absent. For each grid cell of the top 15% resistance level, we kept the resistance level as the value if the Brassicaceae presence probability was ≥0.3, otherwise setting it 0. Then we excluded grids with 0 value and mapped with Arcmap 10.2 (see Fig. 4c, d). R scripts are provided as Supplementary Software.

Publication bias

We checked for publication bias by drawing a funnel plot of the effect size (logRR) vs. study size (sample size in bioassay) and tested this potential bias statistically using Kendall’s rank correlation⁷². We found no indication for a publication bias (see details in Supplementary Methods for publication bias assessment).

Statistical analyses were done in R v4.0.5 and RStudio 1.1.463.

Ethics statement

All applicable international, national, and/or institutional guidelines for the care and use of animals were followed. All procedures performed in studies involving animals were in accordance with the ethical standards or practice of the Institute of plant protection, Chinese Academy of Agriculture Sciences, Beijing, China.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com