1. Introduction

The development of pesticide resistance in invertebrate pests represents rapid contemporary evolution in response to human-induced stressors (ffrench-Constant et al., 2004, McKenzie and Batterham, 1994, Tabashnik et al., 2013, Gould et al., 2018). Since the first instance of pesticide resistance was documented in 1914, cases have exponentially increased to over 17870, involving more than 625 pest species and 360 compounds, and posing threats to agricultural production, human health, and the environment (Sparks et al., 2021, Melander, 1914, Gould et al., 2018, Fritz, 2022). This has led to attempts to stop the evolution of resistance, or at least slow its appearance within populations (Gould et al., 2018). A key task in resistance management is to understand the evolution of resistance in field populations, where a given mutation may have originated once and spread via gene flow or evolved through multiple independent substitutions, and where resistance-associated mutations can pre-date (standing genetic variation) or arise subsequent to (de novo mutation) the introduction of pesticide selection pressures (Ffrench-Constant, 2007, Daborn et al., 2002b, Kersten et al., 2023).

A challenge in understanding the evolution of resistance is that the genetic basis of pesticide resistance observed in laboratory-generated resistant strains often does not reflect that in field populations (Chen et al., 2023, Walsh et al., 2022, Ffrench-Constant, 2013, Legan et al., 2024). Field populations are typically much larger and more genetically heterogeneous than laboratory populations. Because of the greater potential availability of rare mutations with large effect sizes in field populations, monogenic resistance to pesticides involving these rare mutations may occur relatively more frequently than in laboratory-generated strains (McKenzie and Batterham, 1994), though small effect loci can also remain important in field resistance (Schmidt et al., 2024, Lucas et al., 2023). In contrast, artificial selection in the laboratory usually takes place in populations with much smaller effective population sizes, as well as a constant intensity and pattern of selection expected to favor the accumulation of multiple changes of small effect (Ffrench-Constant, 2013). Where alleles with large effects are selected in the laboratory, these may not necessarily reflect the full set of resistant mutations in the field (Legan et al., 2024). Such mismatches between laboratory-selected and field-evolved resistance may lead to an inaccurate understanding of resistance evolution (Ffrench-Constant, 2013, Bras et al., 2022, Legan et al., 2024).

A further challenge in understanding resistance evolution in the field is that the genetic basis of resistance is often examined retrospectively well after the emergence of resistance in the field. Retrospective studies can fail to capture processes that have occurred in the intervening period (Andreev et al., 1999, Raymond et al., 1991, Ffrench-Constant et al., 1993b, Adesanya et al., 2021), such as the introduction of new alleles through gene flow or de novo mutations, and fitness costs as well as changing selection regimes through pesticide rotations that can lead to resistance alleles being transient (Schmidt et al., 2010, Ffrench-Constant, 2013, Kreiner et al., 2022, Kreiner et al., 2019, Groeters and Tabashnik, 2000). Preserved specimens before and following the introduction of pesticides can allow for resistance processes to be investigated retrospectively, but samples are often unavailable (Hartley et al., 2006, D’Costa et al., 2011). Ongoing genetic monitoring of populations can help in overcoming some of these issues, particularly if population processes affecting a species are understood and if the genetic basis and new resistant alleles are identified as they evolve in the field (Van Leeuwen et al., 2012, Pélissié et al., 2022, Pezzini et al., 2024, Taylor et al., 2021).

Here, we take advantage of the recent evolution of resistance in the two-spotted spider mite (TSSM), Tetranychus urticae (Acari: Tetranychidae), to a novel acaricide, cyetpyrafen, to provide a detailed understanding of contemporary resistance evolution. TSSM is a global super pest that has evolved resistance to almost all chemical modes of action and is the most resistant invertebrate pest worldwide (De Rouck et al., 2023, Van Leeuwen et al., 2010). Succinate dehydrogenase inhibitor (SDHi) acaricides are the most recently developed mode of action (IRAC group 25) for mite control, and include four chemicals: cyetpyrafen, cyenopyrafen, cyflumetofen, and pyflubumide (Hayashi et al., 2013, Nakahira, 2011). Laboratory-selected resistance has already been identified to these new acaricides, and is conferred by target-site mutations of the SDH coding sdhB and sdhD genes, with mutations identified at three positions (Njiru et al., 2022, Sugimoto et al., 2020). In China, three acaricidal SDH inhibitors have been registered over the last decade: cyflumetofen (2013), cyetpyrafen (2017) and cyenopyrafen (2019). Since its release in 2017, cyetpyrafen has become the most popular pesticide for controlling spider mites, though low-level resistance was detected in the field only one year after its release (Chen et al., 2019).

We have now tracked resistance through repeatedly monitoring field-evolved resistance of TSSM to cyetpyrafen at a time when resistance has developed rapidly to high levels across geographically distant populations. We investigated resistance using a combination of population genomics, laboratory selection, bioassays, and studies on historically preserved specimens before the introduction of cyetpyrafen. We find that multiple and recurrently evolved mutations as large-effect alleles drive the rapid development of high-level resistance to cyenopyrafen in field populations across a widespread geographical range. We also demonstrate that these mutations likely evolved de novo after the introduction of cyetpyrafen, though they may also have originated from extremely rare mutations segregating in populations before cyetpyrafen selection. Two of the mutations from the field were also identified from a laboratory selection experiment. Our results reveal diverse resistance mechanisms encoded by major effect alleles underlying pesticide resistance in field populations. We emphasize the benefits of integrating genomics with other approaches in longitudinal studies to capture the dynamics of rapid evolutionary changes.

2. Results

2.1 Strong resistance evolves rapidly in field and lab-selected populations

To monitor the evolution and spread of cyetpyrafen resistance in the field, we conducted bioassays on 46 mite population samples collected from various regions in China between 2020 and 2024 (Figure 1, Tables S1-S2). The susceptibility of these population samples to cyetpyrafen was compared with previously reported data for populations collected between 2017 and 2018 during the early release stage of cyetpyrafen in China (Chen et al., 2019, Gong et al., 2017). In 2017-2018, all seven tested populations remained susceptible, with LC₅₀ values ranging from 1.4 to 21.6 mg/L (Figure 1a). In 2020, one population from Shouguang City, Shandong Province (SDSG2), showed resistance with an LC₅₀ > 32,000 mg/L, while the other 13 populations remained susceptible (Figure 1b). By 2021, high levels of resistance to cyetpyrafen quickly became widespread in field populations: eight of 10 tested populations had an LC₅₀ higher than 1606.2 mg/L (Figure 1c). In 2024, 11 of the 13 tested populations evolved resistance with an LC₅₀ higher than 1498.6 mg/L (Figure 1d). Overall, the susceptibility of all populations was bimodal, with populations being highly susceptible (LC₅₀ ranges from 0.61 to 30.1 mg/L, x̄ = 6.72 mg/L) or highly resistant (LC₅₀ > 1498.6 mg/L, x̄ > 5892.5 mg/L) and an 877-fold difference between mean LC₅₀s of the susceptible and resistant groups. While only a single sample was resistant in 2020, by 2021 resistance was observed across the entire range of China (Figure 1c). This pattern indicates that field populations of TSSMs have rapidly developed high levels of resistance to cyetpyrafen across a wide geographical area.

Figure 1

We conducted artificial selection to further investigate resistance evolution, using two laboratory strains derived from the same field population and reared in parallel. The field population for selection was collected from Hangzhou, Zhejiang province (ZJXS1), with an LC₅₀ of 21.6 mg/L. One laboratory strain, ‘LabR’, was exposed repeatedly to cyetpyrafen, while the other, ‘LabS’, was not. LabR underwent a gradual increase in cyetpyrafen resistance for the first 15 selection rounds (31 generations), but a much sharper increase over the remaining 18 selection rounds (35 generations). By the end of the experiment, LabR possessed an LC₅₀ of 54,336 mg/L, a 2,521-fold increase in resistance relative to the ZJXS1 field population. In contrast, the LC₅₀ of the LabS strain decreased to 1.21 mg/L over the same period (Figure 1e, Table S3), due to unknown reasons, perhaps linked to adaptation to the laboratory environment. However, this is a very minor change in resistance compared to the large increase seen in the selected strain.

Additional bioassays were used to test whether the selected LabR strain showed cross-resistance to other pesticides that target proteins in the mitochondrial electron transport chain. For the other two mitochondrial complex II inhibitors (IRAC group 25), the LC₅₀ of cyenopyrafen increased 21,254-fold (3.76 to 79918 mg/L), and the LC₅₀ of cyflumetofen increased 69,861-fold (6.47 to 452,003 mg/L). For the mitochondrial complex III inhibitors, the LC₅₀ of bifenazate (group 20D) increased 2-fold (9.87 to 18 mg/L), and the LC₅₀ of pyridaben (group 21A) increased 3-fold (3,763.4 to 10,352 mg/L). The results indicate that cyetpyrafen selection led to high cross-resistance to other mitochondrial complex II inhibitors but weak cross-resistance to the mitochondrial complex III inhibitors (Table S3).

2.2 Selective sweep signals are clear but heterogeneous

To identify genes potentially associated with cyetpyrafen resistance, we conducted pooled whole-genome resequencing on 22 mite population samples, including historic samples of eight geographical populations collected in 2017 before the first field releases of cyetpyrafen, 12 samples (10 resistant and two susceptible) collected afterwards, and two laboratory populations (Tables S1, S4). To allow for variation in the genetic basis of resistance and genetic backgrounds among populations, we chose the Population Branch Excess (PBE) analysis (Yassin et al., 2016) method to scan selection signals across the genome; we used two susceptible populations (NMHS2 and LabS) as nonfocal populations to correct for the varied genetic backgrounds, and one resistant population as the focal population to identify population-specific signals of pesticide selection. No clear and coincident selection signal was identified among historical population samples (Figure 2a). The laboratory-selected resistant LabR population had clear selection signals at two genes, sdhB and sdhD, that encode two subunits of succinate dehydrogenase, the target of SDHi acaricides (Figure 2b). None of the 10 field-collected resistant populations had selection signals at both of these genes, but four (BJDX6, ZJHZ1, ZJHZ2, and LNDD) had selection signals at sdhB, and one (BJCP4) had a selection signal at sdhD. Five of the resistant populations (QHHD, GXNN, SDSG3, SDQZ and YNYX) showed no selection signal at either gene; this same pattern was observed in the one population sample (HNCS2) that remained susceptible after cyetpyrafen was introduced. Additional PBE peaks were observed at other genomic regions in historical, resistant, and susceptible populations; these could reflect demographic processes or other selective forces not considered in this study.

Figure 2

Selection for cyetpyrafen resistance could impact patterns of genetic diversity in the regions of interest or more widely across the genome if selection decreases overall population size. We investigated these impacts focusing on genome-wide patterns of nucleotide diversity (π) and on Tajima’s D. Five resistant populations with selection signals identified by PBE analysis (Figure 2b) showed corresponding decreases in π and Tajima’s D near the two SDH genes (Figure S1). These patterns were not observed in the overall genetic diversity across the genome (Figure S2, Table S4), suggesting that resistant populations had not experienced strong bottlenecks. Five resistant populations without selection signals did not show obvious reduced genetic diversity or Tajima’s D values around the SDH genes (Figure S1) suggesting that these had not experienced strong selection or selective sweeps.

2.3 Multiple mutations have arisen in the field

To identify specific mutations conferring resistance, we examined the two SDH genes, sdhB and sdhD, where selective sweep signals were detected in some populations. We used whole-genome resequencing data from pooled individuals; we also conducted Sanger sequencing of partial sdhB (216 bp of 1001 bp) and sdhD (155 bp of 1,008 bp) genes to validate the pool-seq data and examined mutation frequencies in additional populations; this involved 1268 individuals from 52 populations, including samples collected before and after the release of cyetpyrafen (Table S1). Allele frequencies identified by pool-seq and Sanger sequencing were correlated (Figure S4). A total of 15 amino acid substitutions were identified, including six mutations on five amino acid residues of sdhB: H146Q (CAT to CAA), S212I (AGC to ATC), H258Y (CAT to TAT), I260T (ATC to ACC), I260V (ATC to GTC), and A285S (GCA to TCA), and nine mutations on three amino acid residues of sdhD: D116G (GAC to GGC), D116E (GAC/T to GAA), D116N (GAC to AAC), R119C (CGT to TGT), R119L (CGT to CTT), R119G (CGT to GGT), R119P (CGT to CCT), R119H (CGT to CAT), and P120L (CCC to CTC) (Figure 3b, Table S2). The residue R119 on sdhD was a notable extreme case, involving five amino acid changes. The second codon position of R119 (G) was substituted by the other three nucleotides, while the first codon position (C) was substituted by two other nucleotides (T and G), indicating that the first two codon positions of R119 are almost saturated in terms of possible nucleotide substitutions.

Figure 3

Multiple mutations were commonly found in a single population, while none of the mutations were fixed in any field population (Figure 2b, Table S2). Populations sampled in 2020 had at most two mutations, while those sampled in 2021 and 2024 had at most five. Populations sampled in 2020, 2021 and 2024 had on average 1.17, 2.46 and 2.55 mutations, respectively. We also examined the other two SDH subunit genes, sdhA and sdhC, using pool-seq data, but no amino acid replacements were identified.

We checked the evolutionary conservation of these mutations across animals by aligning homologous sequences from ∼700 animal species, encompassing nematodes, insects, fishes, and mammals (Figs. 3a and 3b). The protein sequences of sdhB and sdhD were highly conserved among various animal species, particularly around the mutant residues H146, S212, I260, D116, and R119 observed in this study. In a structural model of human SDH, eight amino acid residues were predicted at the ubiquinone binding site (Q-site), corresponding to 211, 215 and 260 on mite subunit sdhB, 58, 63, 67 and 71 on sdhC, and 117 on sdhD (Njiru et al., 2022, Du et al., 2023). All detected mutations in the mite are located around or in the Q-site, although differences exist between human and mite sequences (Figure 3c).

Although cyetpyrafen only started to be used in 2017, it is possible that these potential resistance mutations were present as standing genetic variation before this time. To investigate this, we obtained Sanger sequencing and pool-seq data from 2317 female mites collected before the cyetpyrafen release and 4631 after its release (Figure 3d, Table S2). Of the 15 putative resistance mutations, only one (A285S) was observed in the pre-release specimens, which was present in population SDRZ (unclear resistance status), but was not detected in specimens collected after the cyetpyrafen release. Seven mutations were first detected in 2020 (I260T, D116G, D116N, R119C, R119L, R119P, and R119H), three mutations were first detected in 2021 (H146Q, S212I, and I260V), and the remaining four mutations were first detected in 2024 (H258Y, D116E, R119G, and P120L) (Figure 3d, Table S2). Among the target-site mutations, I260V and R119C were the most common, found in 65.9% and 46.3% (respectively) of the populations collected after the release of cyetpyrafen, with their frequency increasing annually (Figure 3d). In the laboratory-selected population LabR, two mutations (I260V and R119C) were fixed, but neither was detected in either the field population from which it was sourced (ZJXS1) or the unselected susceptible laboratory population derived from ZJXS1 (LabS). These results indicate that the putative resistance mutations were not likely to be present as standing genetic variation unless they were segregating at very low frequencies, and they may have evolved recently through de novo substitutions although this cannot be demonstrated definitively.

2.4 Recurrent evolution of identical mutations

When a specific resistance mutation is found in multiple populations, this can be the result of migration between populations or to multiple independent origins of the same mutation. We evaluated these evolutionary hypotheses using two methods: TreeMix (Pickrell and Pritchard, 2012) and haplotype network analysis.

We used TreeMix to investigate whether populations with identical alleles were more genetically similar or had recent histories of migration. TreeMix builds maximum likelihood trees from population allele frequencies and then adds optional migration edges between populations. We set the SCCD population from southwest China as the root population as it had the highest genetic diversity, and we built trees that followed a strict drift model (i.e. no migration) and trees with up to 10 migration edges. The strict drift model explained most of the variance in relatedness among populations (mean over 1000 replicates of 0.916). The best number of historical migration events was one. However, this migration edge was placed between the laboratory-selected population and a field clade different to the one used for the laboratory population (Figure 4a), an impossible evolutionary scenario, suggesting that migration among populations is low. The consensus tree over the 1000 replicates assigned most of the populations into northern and southern groups by geography, but not grouped by resistant status (Figure 4a). We also observed several instances where geographically distant populations grouped together. Some southern populations were placed with the northern group (YNYX), and vice versa (BJDX6); these populations may have been recently established via long-distance dispersal, possibly by the transport of seedlings.

Figure 4

The above-mentioned mutations in sdhB and sdhD genes were detected in six populations from the southern group and seven populations from the northern group. The mutation I260V was detected in every field population sampled after the release of cyetpyrafen (Figure 4a, Table S2), using pool-seq data. I260T was only detected in the northern group and R119L was mostly detected in the southern group. The mutation H146Q was found in four adjacent populations (ZJHZ1, ZJHZ2, GXNN and BJCP4). This geographical heterogeneity in the distribution of resistance mutations suggests multiple mechanisms of resistance evolution underly these patterns, involving spread via migration between populations (e.g. H146Q) as well as possible independent origins of the same mutation (e.g. I260V).

We explored these ideas further using haplotype network analysis. When resistance spreads via migration, these mutations should have identical haplotypes, unless recombination has taken place which will introduce a second haplotype that segregates locally alongside the initial haplotype. We constructed the haplotype network for the partial sdhB (Figure 4b) and sdhD (Figure 4c) genes determined from Sanger sequencing. Since haplotypes could not be determined when two loci were heterozygous, we detected haplotypes from sequencing data with at most one heterozygous locus. The genes showed high genetic diversity, with 11 (sdhB) and 24 (sdhD) haplotypes detected. We assumed that common haplotypes located at the centre of the network are the ancestral haplotypes (H1 and H4 in sdhB; H1, H2, H3, and H11 in sdhD). We assumed all other haplotypes were derived. Each derived haplotype was only one nucleotide different from the nearest ancestral haplotype. Five of the nine derived sdhB haplotypes carried a putative resistance mutation, as did 13 of the 22 derived sdhD haplotypes.

Of the 10 putative resistance mutations detected through Sanger sequencing, all but four (H258Y, D116N, R119G, and P120L) were found as multiple haplotypes (Figs. 4b and 4c). To test whether these haplotypes were produced by multiple origins or by recombination, we inferred the possible means by which recombination could produce each haplotype, and investigated the geographical distributions of the requisite ancestral and derived haplotypes. For mutations in sdhB, we inferred that recombination would have to involve the H4 ancestral haplotype as donor, while the H11 ancestral haplotype was required for sdhD (Figure S3). Thus, recombination could have produced the sdhB haplotypes H5 (H4 x H6) and H8 (H4 x H11), and the sdhD haplotypes H17 (H11 x H10), H16 (H11 x H13) and H21 (H11 x H20). However, only three populations (GXNN, YNKM3, and YNYX) had both of the two recombining haplotypes required and one derived haplotype (Table S5). Additionally, the proximity of the mutations was such that multiple recombination events would be required within the same 5 – 60 bp regions, which seems unlikely even with a high recombination rate of (say) 10 cM/Mb where two mutations might only be expected to occur once in 3.6 x 10⁹ individuals across 60 bp. These results do not indicate an obvious role for recombination in generating the derived haplotypes.

All haplotypes carrying putative resistance mutations in the sdhB gene arose in 2020-2021, along with all haplotypes of the sdhD gene apart from H6, H7, H16 and H23 (Figure 4d, Table S5). By 2024, the sdhB H6 haplotype became the most prevalent haplotype, while the non-resistant sdhB H3 haplotype was almost lost despite it being prevalent between 2013 and 2021 (Figure 4d, Table S5). Some individual populations had multiple haplotypes of the same putative resistance mutation (Figure 4d), and different haplotypes of the same mutation were often found in adjacent populations (Figure 4a, Figure 4d, Table S5). In the ZJXS1 population used to produce the LabR and LabS strains, four haplotypes of the sdhB gene (H1, H3, H4 and H7) and six haplotypes of the sdhD gene (H1-H4, H11 and H24) were found, none of which had putative resistance mutations. After 66 generations of cyetpyrafen selection, two resistant haplotypes were fixed in the LabR population: H6 of sdhB and H10 of sdhD, which were likely derived from H1 of sdhB and H2 of sdhD, respectively (Table S5). These results illustrate the rapid and geographically widespread emergence of resistant alleles in these populations, linked to several mutations found in a range of different haplotypes.

2.5 Single mutations have a large effect on resistance

We next characterized the function of each mutation by itself and in combination with other mutations. We first tested the population-level association between allele frequencies of the mutations and levels of resistance, using the five mutations present in at least ten populations (I260T, I260V, D116G, R119C, and R119L) and resistance indicated by survival at 1000 mg/L cyetpyrafen. When each mutation was analyzed separately, a weak association was found between survival and the frequency of the mutant allele (Table S6, linear regression R² = 0.091-0.303). When we considered the frequency of individuals having any of the five mutations, we found a strong positive association between survival and the frequency of predominant allele (R² = 0.707), the frequency of individuals carrying at least one resistant allele (heterozygous and homozygous mutant genotypes) (R² = 0.747), and the frequency of individuals with homozygous mutant genotypes (R² = 0.606) (Figure 5a-5c). These results showed that single mutations explained a low proportion of resistance variation across populations. But when considering all five mutations together, a substantial proportion of resistance variation was explained.

Figure 5

Near iso-female lines were created and tested for susceptibility to cyetpyrafen (Figure 5d). These carried I260V, D116G, R119E, R119L, and one selected line carried both I260V and R119C. When treated with 1000 mg/L cyetpyrafen, mortality of the susceptible line was 100%, but was lower than 25.2% in all other lines, indicating that each mutation and a combination of two mutations conferred high resistance to cyetpyrafen. These results demonstrate that cyetpyrafen resistance is often conferred by different large-effect mutations in the sdhB and sdhD genes.

3. Discussion

Pesticide resistance can evolve rapidly, in contrast to traditional views of adaptation which often emphasize a gradual process (Fisher, 1930). Due to the time taken for resistance to be recognized, evolutionary studies of pesticide resistance are often post hoc. Here, we conducted field monitoring of resistance to a novel acaricide, cyetpyrafen, following its 2017 release, with the opportunity to track resistance as it develops. We demonstrated a longitudinal picture of resistance development in the field involving multiple and recurrently evolved large-effect mutations.

Currently available pesticides mainly act by direct binding to receptors and enzymes, resulting in the inhibition of physiological functions associated with the nervous system, energy production system, growth and development, and midgut absorption (Sparks et al., 2021). Mutations that alter the binding kinetics of pesticides on these proteins and thus potentially confer pesticide resistance have been widely identified in agricultural pests and disease vectors (De Rouck et al., 2023, Acford-Palmer et al., 2023). Resistant mutations of some target genes are highly conserved both among populations and among species, perhaps due to the functional constraints of these targets (Ffrench-Constant et al., 1993a, Hawkins et al., 2018). However, multiple mutations on single pesticide targets have been identified in the same pest species and among different species. These cases often involve the voltage gated sodium channel (VGSC) gene, where at least 5 and 10 mutations have been cumulatively reported across studies on TSSM and the ectoparasite Dermanyssus gallinae respectively (De Rouck et al., 2023).

Here we found 15 mutations on eight amino acid residues in TSSM that all arose within a few years following the introduction of cyetpyrafen. Notably, five of the substitutions involve R119 of the sdhD gene. As far as we know, this is the highest number of both mutations and altered residues ever found in a single pesticide target from the same pest species, all emerging in a short period following initial detection of resistance. These mutations were located at the Q-site of SDH targeted by acaricidal inhibitors of SDH. Among these mutations, I260T was previously found in a cyflumetofen-resistant population (Sugimoto et al., 2020), while I260V (with S56L in subunit C) was found in a pyflubumide and cyenopyrafen-resistant population (Maeoka and Osakabe, 2021). R119L has previously been reported in cyflumetofen resistant and susceptible strains of TSSM but had not yet been validated as having a role in resistance (Sugimoto et al., 2020). These results provide information on the nature of novel mutations in a pesticide target and the high level of resistance to cyetpyrafen that they can confer. It would be interesting to further validate the impacts of the mutations across chemicals to establish patterns of cross resistance given that the same mutation may confer a different level of resistance to different chemicals (Maeoka and Osakabe, 2021).

Our temporal collections of TSSM preceding cyetpyrafen release in the field indicated that the resistant mutations were not at an appreciable frequency prior to pesticide selection. Sequencing of 13 populations (involving a total of 2666 TSSMs) from 2013 to 2018 did not recover any of the identified mutations. This points to a possible role of recurrent de novo mutation, although resistant alleles could have been present at very low frequencies before pesticide selection which can be difficult to detect (Hawkins et al., 2018). In laboratory populations of TSSM, continuous selection pressure resulted in the appearance of the de novo mutation H258Y, which confers resistance to pyflubumide and cyenopyrafen (Njiru et al., 2023). Our selection experiments also led to new mutations (I260V, R119C) arising over 66 generations (2 years) of cyetpyrafen selection in the LabR population. Although it is hard to completely exclude the possibility of contamination during laboratory selection, because TSSM is tiny and can easily be dispersed passively, the haplotype analysis showed that the two mutations in the lab-selected population likely arose de novo. In the same way, multiple de novo mutations could have contributed to the rapid evolution of resistance in different populations and regions as documented in our field surveys.

Although few studies have traced resistance evolution in the field, they suggest that both de novo mutations and standing genetic variation can contribute to resistance, with their relative importance depending on the specific pesticides and pests involved (Hawkins et al., 2018). Resistant alleles can originate from a single event and spread to distant populations by migration (Raymond et al., 1991, Daborn et al., 2002a), but resistant mutations can also have multiple origins (Guan et al., 2021, Ilias et al., 2014, Troczka et al., 2012, Kreiner et al., 2022, Schmidt et al., 2024). Direct evidence for multiple origins comes from the presence of the same mutation on multiple haplotypes, assuming that recombination can be excluded as an explanation because insufficient time has elapsed after a resistant mutation is introduced into a new background (Semenov et al., 2019). In our study, the requisite recombinant haplotypes were rarely found in the same population, and we therefore suspect that recombination is unlikely to be the source of patterns of haplotype variation over short DNA sequences. We suspect that populations of mites like TSSM have a strong evolutionary potential to evolve in situ rather than relying on gene flow, and that there is a repeatability of evolutionary changes when an organism is exposed to the same selection pressure, a long-standing issue in evolutionary biology (Konecna et al., 2021). The recurrent appearance of the same mutation in different genetic backgrounds suggests that the likelihood of resistance evolution is high in TSSM and not dependent on gene flow, and this is further increased by the fact that resistance can develop through several mutational options.

Because field populations of TSSM can use multiple mutations to reach the same phenotypic endpoint (resistance), it is unlikely that all resistant mutations will be identified from laboratory selection studies particularly if these involve a single laboratory population. Even in the case where multiple populations of TSSM have been selected for resistance (to cyenopyrafen and the related chemicals cyflumetofen and pyflubumide), QTL mapping still only identified a few of the relevant mutations providing resistance to each chemical (Sugimoto et al., 2020). Additionally, mutations identified in laboratory selection may not be involved in resistance in the field because of fitness costs. For instance, H258Y conferred resistance to cyenopyrafen and pyflubumide in laboratory-selected populations, but likely reinforced cyflumetofen binding and toxicity (Njiru et al., 2022) to produce a fitness cost (Njiru et al., 2023). Relying on a limited set of resistance alleles identified in the laboratory for monitoring may bias inferences of resistance evolution in the field when many alleles can produce resistance. Genetic and evolutionary studies on field populations are then essential to develop a full picture of resistance evolution (Konecna et al., 2021, Lucas et al., 2023, Schmidt et al., 2024, Kersten et al., 2023, Pezzini et al., 2024).

Theoretical and empirical studies have showed that field-evolved resistance often involves large-effect alleles with a simple genetic basis, while laboratory selection may result in polygenic genetic architectures, due to differences in the occurrence of rare alleles (more likely in large field populations) and different selection pressures (more variable with intense periods in the field on occasions) (Ffrench-Constant, 2013, Roush and McKenzie, 1987, Ffrench-Constant et al., 1990, McKenzie and Batterham, 1994, Crow, 1954). This could potentially be mitigated by establishing laboratory selection lines using large starting populations from the field with pre-existing resistance alleles (Gould et al. 1997). Until recently, few studies have directly compared the genetic basis of resistance in laboratory populations with field-evolved resistance across a large geographical scale (McKenzie et al., 1992, Ffrench-Constant et al., 1990, Crow, 1954, Legan et al., 2024). Here we identified large-effect alleles in both field- and laboratory-selected resistant populations. High-level resistance was generated after lab selection with cyetpyrafen for 33 of 65 generations, almost at the same rate as the field populations where high-level resistance was observed three years after the commercial release of cyetpyrafen. In field populations, it seems that multiple allelic options can favored by selection due to parallel adaptation, and that the same allele may also end up residing in different haplotypes in different areas. Through migration these alleles/backgrounds may eventually end up in the same population to produce a polyallelic basis of adaptation where particular alleles may then not reach high frequencies unless one haplotype is favored (Ralph and Coop, 2010).

Pesticide resistance provides many fascinating examples of rapid and recent evolution under high selection pressure, which can be leveraged to address fundamental issues about the evolutionary basis of adaptations to new environments (Hawkins et al., 2019), and also examine the reliability of commonly used research methods. The development of high-throughput sequencing and population genomics methodologies provides an opportunity to detect selective signatures across the genome (Van Etten et al., 2020, Cohen et al., 2022, Pélissié et al., 2022). Although pesticide resistance is a rapid evolutionary process under strong selective pressures, identifying candidate genes directly from genomic signals remains difficult when adaptation occurs via soft sweeps (Pennings and Hermisson, 2006, Schrider et al., 2015). Soft sweeps may be common in adaptive evolution if the effective mutation rate (2N_eμ) >> 0.01, e.g., through large population sizes (N_e), and/or high allelic mutation rates (Hermisson and Pennings, 2005, Pennings and Hermisson, 2006). The multiple haplotypes and distinct evolution history of SDH genes among resistant field populations suggest a prevalence of soft selective sweeps through recurrent mutation as implicated in pesticide resistance in other pest systems, such as Drosophila (Garud and Petrov, 2016), and the oomycete Plasmopara viticola (Delmas et al., 2017). This is congruent with the ‘messier’ signals found in field populations relative to the artificially selected laboratory resistant population. Sweeps around the SDH genes are incomplete, or other contributing loci are undergoing selective sweeps, which may further reduce signal to detect adaptation (Pennings and Hermisson, 2006, Schrider et al., 2015). Fewer mutation types and haplotypes along with higher mutation frequency are more likely to be detected. Clearer signals of selective sweeps could emerge from haplotype-based linkage-disequilibrium methods using individual-level genomic data (Hermisson and Pennings, 2017).

In conclusion, we have captured a detailed picture of contemporary evolutionary processes under pesticide selection in field populations of an invertebrate pest. We demonstrate that the rapid evolution of pesticide resistance is driven by an unprecedented number of population-specific mutations on target genes that arose recurrently. These findings reveal complex origins and diverse genetic bases of pesticide resistance in the fields, challenging the molecular monitoring of resistance and the implementation of effective resistance management strategies. Finally, we call for future studies on field-evolved resistance to provide theoretical and practical insights into effective pest control, leveraging the genomics era.

4. Materials and Methods

4.1 Samples

We used 72 field-collected and laboratory-selected populations of TSSM for bioassays, whole-genome resequencing, Sanger sequencing of target genes, and establishment of mutant strains (Table S1). Samples of TSSM were collected across China from 2013 to 2017 and from 2018 to 2024 before and after the release of cyetpyrafen. More than 500 adults were collected from at least three points separated by 100 meters for each population. Samples were kept in 100% alcohol or transported live to the laboratory for chemical bioassays. Because TSSM has a haplodiploid sex-determination system, only female adults were used for bioassays and genotyping. A field population, ZJXS1, was collected in 2018 and selected for 60 generations by spraying cyetpyrafen. The TSSM strain were reared on French bean, Pemphigus vulgaris L. under 25±1 °C, 60±5% relative humidity and photoperiod of L16h: D8h. We continuously selected this population in the laboratory for another 40 generations by killing 60% of individuals every two generations, to obtain a laboratory-selected resistant strain.

4.2 Bioassay

Bioassays were conducted to test the susceptibility of TSSM to acaricides using a Potter Spray Tower (Burkard Scientific, London, UK). We used transparent plastic cups (6 cm in diameter and 4 cm in height) to keep the female adult mites. A layer of 0.2% agar was placed on the bottom of the cups to avoid leaves drying out. Fresh leaves of P. vulgaris were cut into a circle with a diameter of 6 cm to fit the cup and put onto the agar. The borders of the leaves were sealed using 0.2% agar. Thirty female adults were carefully moved into a plastic cup by suction trap. Two ml acaricide solution was applied through a Potter spray tower (Burkard Manufacturing Co. Ltd., UK) at 68.9 kPa and a sedimentation time of 30 s. Each population was exposed to between five and eight pre-estimated concentrations of acaricides, diluted with water containing 0.1% Triton X-100 (Beijing Solar BioScience and Technology Limited Company, China). The 0.1% Triton X-100 solution was used as a control treatment to correct the potential influence of adjuvants in the acaricides. The treated mites were kept under 25 °C, 60-70% relative humidity, and a photoperiod of 16 h:8 h (L:D). Mortality was assessed under a Stemi 305 stereomicroscope (Zeiss, Germany) 48 hours after spraying. The testing TSSMs were considered dead if no movement of appendages was observed when they were prodded with a fine brush at 24 hours and 48 hours after spraying. The 50% lethal concentration (LC₅₀) and its 95% confidence intervals were calculated using probit regressions implemented in the statistical software DPS v12.01 (Tang and Zhang, 2013).

4.3 Pooled whole-genome resequencing and SNP calling

Between 50 and 400 individuals per population were used for pooled genome resequencing to obtain population genomic data for diversity estimates and demographic analyses. DNA libraries with insert sizes of 400 bp were constructed and sequenced on the Illumina NovaSeq 6000 platform to obtain 150 bp paired-end reads. Raw reads were filtered using fastp v0.23.0 (Chen et al., 2018). The filtered reads were mapped to our recently assembled chromosome-level genome (Cao et al., 2024) using bwa-mem2 2.2.1 with default parameters (Vasimuddin et al., 2019). SAMtools v1.14 (Li et al., 2009) was used to sort, remove duplication reads, remove reads with mapping quality < 20, and convert bam files to mpileup files. Indels in mpileup files were filtered out (with a 5-bp margin on each side) to avoid spurious SNP calling using scripts from PoPoolation (Kofler et al., 2011a). Then mpileup files were converted using scripts of PoPoolation2 (Kofler et al., 2011b).

4.4 Population genetic diversity and demography analyses

Genetic diversity in each population was quantified by estimating nucleotide diversity (π) and Tajima’s D across 5 kb non-overlapping windows, using PoPoolation (Kofler et al., 2011a). To reduce the impact of read depth on estimates of nucleotide diversity, each alignment was converted into a population-specific pileup file using SAMtools, and reads were subsampled without replacement up to a depth of 50 (Hoban et al., 2016). The mean F_ST across all windows for each pairwise population comparison was calculated using PoPoolation2 (Kofler et al., 2011a) across 5 kb non-overlapping windows. To infer population relationships and historical migration, a rooted maximum likelihood tree was constructed with a strict drift model (no migration) using Treemix v1.13 (Pickrell and Pritchard, 2012). The SCCD population was set as root and we tested migration events from 0 to 10 with at least 100 runs for each case. The best number of migration events was determined by the second-order rate of change in likelihood (Δm) across incremental migration events (m) using the Evanno method (Evanno et al., 2005) embedded in R package ‘optM’ (Fitak, 2021). For the best m, we performed 1000 runs to generate 1000 maximum likelihood trees and obtained a consensus tree using the program SumTrees version 4.10 (Sukumaran and Holder, 2010).

4.5 Genome scan of outlier SNPs

We investigated the genomic signals of selection around the four SDH genes (± 1 Mb) and across the whole genome. The main statistic used to identify outlier regions was the population branch excess (PBE) statistic (Yassin et al., 2016). The PBE statistic is based on comparisons of F_ST values among three populations, comprising a focal population and two nonfocal populations (Pool et al., 2017). In the analysis, each resistant population was set as the focal population, and two susceptible populations (LabS and NMHH2) were set as the nonfocal populations. PBE values were calculated for 5 kb non-overlapping sliding windows. Windows with top 5% PBE value were considered outliers. Patterns of nucleotide diversity (π) estimates and Tajima’s D values around the SDH genes were also examined. We expected resistant populations to have lower nucleotide diversity and Tajima’s D values near the SDH genes due to directional selection favoring resistant haplotypes. The coding effects of SNPs within the SDH genes were obtained using SnpEff v4.3 (Cingolani et al., 2012).

4.6 SDH homology modeling

The three-dimensional (3D) structures of the wild and mutational SDH from TSSM were modeled with SWISS-MODEL (Waterhouse et al., 2018) using 8GS8.1 (Du et al., 2023) as a template. To examine the conservation of SdhB and SdhD genes, we collected the protein sequences of two genes of different animals in the protein database of NCBI (https://www.ncbi.nlm.nih.gov/protein). Then, MUSCLE in MEGA v11 (Tamura et al., 2021) were used to perform sequence alignment and trim protein sequences around mutation points. After removing identical sequences of the same species, the multiple sequence alignment sequences were used to determine conservation and visualized with WebLogo 3 (Crooks et al., 2004).

4.7 PCR amplification and sequencing of target genes

Sanger sequencing was used to characterize the spatiotemporal variation of these target-site mutations in populations before and after the release of cyetpyrafen (Table S1). PCR primers were designed for sequences of the sdhB gene (sdhb_F: GAGTCTAGCCGCAGTCTGAT, sdhb_R: ACTTGGGACCTGCTGTACTC) and the sdhD gene (sdhd_F: GCGATGGATAAGGCATAGACG, sdhd_R: TCTGGTCAGTGTGGTTCCTC) to capture the candidate mutation sites. Polymerase chain reaction (PCR) was conducted in the Mastercycler pro system (Eppendorf, Germany) under an annealing temperature of 54 °C for 1 min using Promega Go Taq™ Master Mixes (Madison, Wisconsin, USA). Amplified products were purified and sequenced directly from one strand using an ABI 3730xl DNA Analyzer by Sangon Biotech Co. Ltd. (Shanghai, China).

4.8 Haplotype analysis

Haplotypes of each gene were identified from Sanger sequencing. We constructed haplotype networks for sdhB and sdhD from individuals that contained homozygous resistant genotypes. Haplotype networks were constructed using the R package pegas (Paradis, 2010), and were visualized with the R package ggplot2.

4.9 Association analysis between mutation frequency and resistance level

Regressions were run to assess the strength of association between target-site mutation frequencies and mean population-level resistance phenotypes. Since LC₅₀ of some populations could not be calculated or were inaccurate due to high resistance levels, the survival rate after spraying with 1000 mg/L cyetpyrafen, which is about 6.7 to 10-fold of field-recommended dose, was used to quantify the resistance level of TSSM to cyetpyrafen. Because multiple target-site mutations existed, we evaluated the associations using different metrics of target-site mutation frequencies: (i) the frequency of the predominant mutation in each population, (ii) the frequency of individuals carrying any mutation in each population, and (iii) the frequency of individuals with homozygous mutations in each population. Regressions were fitted as linear models using the lm function from the R package stats and visualized with ggplot2. We did not adjust associations for relatedness among populations because the evolutionary events leading to resistance appear to be largely independent.

4.10 Establishment of near isofemale lines

Single immature TSSM females were picked from resistant populations and reared individually in a petri dish. When unmated, females will only produce haploid male eggs. To obtain the genotype of males for the sdh genes, the DNA of each female was extracted and genotyped as mentioned above after it had laid more than 10 eggs, while the eggs (> 10) were allowed to develop into male adults (F₁ male adults). To obtain the genotype of virgin female, the first six eggs were genotyped, while the female were kept (F₀ female). Then, a pair of F₀ female and F₁ male adults with the same homozygous genotypes for the three mutation sites were placed into the same petri dish to allow them to mate and produce male and female progeny, with the subsequent line defined as a near isofemale line. The mated male and female pairs were from the same population to reduce differences in their genetic backgrounds. The near isofemale lines with homozygous genotypes at the mutation sites were tested for susceptibility to cyetpyrafen using the bioassay method.

Data availability statement

Raw reads of whole genome resequencing data were submitted to NCBI Sequence Read Archive, with accession numbers SRR30666541 - SRR30666561 and SRR28000465. All scripts and data used in the study were deposited to GitHub: https://github.com/caolijun/TU_pool-seq and Figshare: 10.6084/m9.figshare.27231003 (Wei and Cao, 2025).

Acknowledgements

We thank Ming-Liang Li, Hong-Ping Tian, Liang-Bin Zeng, for their assistance with sample collections. This study is supported by National Key R&D Program of China (2023YFD1401200) and Program of Beijing Academy of Agriculture and Forestry Sciences (QNJJ202006, JKZX202208). Ary Hoffmann and Joshua Thia are supported by the Grains Research and Development Corporation and Horticulture Innovation (Australia).

Additional information

Author contributions

S.J.W. designed research; J.C.C., L.J.C., L.X.Z., Y.Y, M.C.Y., J.Y.Z., X.Y.Z., and Y.J.G., performed research; L.J.C., S.J.W., J.A.T, T.L.S, Q.Y., and A.A.H analyzed data; S.J.W., A.A.H., L.J.C., J.A.T., T.L.S., H.L., R.F.C., and X.X.C. wrote the paper.