Selective Sweeps in a Nutshell: The Genomic Footprint of Rapid Insecticide Resistance Evolution in the Almond Agroecosystem

Abstract Among the most familiar forms of human-driven evolution on ecological time scales is the rapid acquisition of resistance to pesticides by insects. Since the widespread adoption of synthetic organic insecticides in the mid-twentieth century, over 500 arthropod species have evolved resistance to at least one insecticide. Efforts to determine the genetic bases of insecticide resistance have historically focused on individual loci, but the availability of genomic tools has facilitated the screening of genome-wide characteristics. We resequenced three contemporary populations of the navel orangeworm (Amyelois transitella), the principal pest of almond orchards in California, differing in bifenthrin resistance status to examine insecticide-induced changes in the population genomic landscape of this species. We detected an exceptionally large region with virtually no polymorphisms, extending to up to 1.3 Mb in the resistant population. This selective sweep includes genes associated with pyrethroid and DDT resistance, including a cytochrome P450 gene cluster and the gene encoding the voltage-gated sodium channel para. Moreover, the sequence along the sweep is nearly identical in the genome assembled from a population founded in 1966, suggesting that the foundation for insecticide resistance may date back a half-century, when California’s Central Valley experienced massive area-wide applications of DDT for pest control.


Introduction
Among the most familiar forms of human-driven evolution on ecological time scales is the rapid acquisition of resistance to pesticides in agriculture (Hendry et al. 2017). Compared with selection exerted by most natural mortality factors, pesticide selection is generally extremely strong and specific, and resistance can develop exceptionally fast (Hawkins et al. 2019).
Since the widespread adoption of synthetic organic Significance Although it is widely known that insecticide resistance is rapidly acquired by target insects, comprehensive population genomic scans that assess the evolutionary trajectory of the development of resistance are scant. This is particularly true for insect pests that are nonmodel organisms, some of which are especially important due to the substantial damage they can cause to crops and as consequence to food supply. This study looks at the genome footprints and consequences of intensive insecticide usage in a pest of almonds in California. We detected a large selective sweep, an unequivocal signature of positive selection, in a region of the insect genome comprising several genes with potential roles in target insensitivity and insecticide detoxification.
insecticides in the mid-twentieth century, more than 500 arthropod species, over half of which are agricultural pests, are known to have evolved resistance to at least one insecticide (Whalon et al. 2008;Sparks and Nauen 2015). Resistance evolution remains an enduring challenge to the use of synthetic insecticides in agriculture and human disease vector control. Historically, efforts to determine the genetic and evolutionary bases of insecticide resistance have focused on allelic variation and gene expression changes at individual loci associated with metabolic resistance or target-site insensitivity. At the genomic and population scale, studies of selection for insecticide resistance have focused mainly on drosophilid model species (Steele et al. 2015;Schmidt et al. 2017;Duneau et al. 2018) or medically important mosquito species Kamdem et al. 2017;Kotsakiozi et al. 2017). Surprisingly few genome-wide population studies of field-acquired insecticide resistance have been carried out on agricultural pests, and little is known about the evolutionary origins and dynamics that govern this process (Oakeshott et al. 2003;P elissi e et al. 2018).
By querying the full genome and its variation among individuals and among populations, population genomic approaches can illuminate multiple specific targets of positive selection, as well as whole affected regions. A genome-wide approach can also help to distinguish this positive selection from the effect of the demographic history of the populations under examination (Akey et al. 2004). A classic hallmark of positive selection is a reduction in nucleotide diversity in genome regions flanking genes that have been the targets of such selection (Kaplan et al. 1989). This shift in allele frequencies, or "selective sweep" (Kim and Stephan 2002), depends on recombination rates in the genomic region and on the timing and strength of the underlying selective pressure. Classic selective sweeps, or "hard sweeps," where heterozygosity reaches near-zero values at or very close to the causative mutation, are rare in nature and are not readily reproducible in laboratory settings (i.e., experimental evolution studies) which often result in soft sweeps (Burke 2012;Messer and Petrov 2013).
Production of almonds (Prunus dulcis), a multibillion-dollar industry centered in the California Central Valley, can be severely affected by the navel orangeworm (Amyelois transitella, Lepidoptera: Pyralidae), the most important pest of this commodity. Insecticides are routinely and intensively used in almond orchards to manage this pest. Although insecticides approved for use include representatives from five structural classes, pyrethroids have been widely used in recent years due to their low cost and high efficacy. Prior to 2013, insecticide resistance was not known to occur in navel orangeworm populations; in 2013, however, 10-fold resistance to bifenthrin was documented in populations in Kern County (Demkovich et al. 2015). A reference genome for A. transitella was assembled in 2013 (https://i5k.nal.usda.gov, last accesed November 1, 2020). This reference was based on the "SPIRL-1966" line established in the USDA laboratory facility in the San Joaquin Valley in 1966, during a period preceding the use of pyrethroids and during which the Central Valley was subjected via crop-dusting to "thousands of tons of DDT alone" (Cory et al. 1971). DDT shares the same metabolic target as pyrethroids and was widely used from 1945 until 1972. Early use of DDT was associated to preadaptation for resistance to pyrethroid insecticides in several mosquito populations (Rongsriyam and Busvine 1975;Chadwick et al. 1977;Hemingway et al. 1989).
In this study, we took advantage of the availability of the A. transitella reference genome and of the detailed record of past insecticide usage in the California Central Valley to interrogate the origins and evolutionary trajectory of insecticide resistance in this pest at the genomic level. We evaluated nucleotide diversity and genomic differentiation in three modern, geographically separate A. transitella populations with different histories of bifenthrin use and levels of bifenthrin resistance. We also assessed gene expression changes in resistant and susceptible genotypes. We demonstrate the existence of a large selective sweep which includes, among others, genes encoding for cytochrome P450s and the para gene, the target of pyrethroids and DDT.

Results
Nucleotide Diversity and Tajima's D  (table 1). Nucleotide diversity was calculated with the p statistics and the deviation of variability from null expectations under a neutral model with Tajima's D, both metrics were calculated across nonoverlapping 5-kb-long windows for each of the populations (table 1). To identify regions of low nucleotide diversity (low p), we sorted the 5-kb regions according to their calculated p values. The lowest p values were compared with those across the genome and with averages in nearby regions. In addition, the windows with the lowest p values were assessed for their calculated Tajima's D values at the 90% CI in the same region. With this analysis, an unusually large region displaying reduced genetic diversity across the three populations was readily detected ( fig. 2). This region showed p values as low as zero in the R347 population and ranging from 0.00002 to 0. The distribution of Tajima's D values for 5-kb windows for each of the three populations ( fig. 3) shows heavy tails at both positive and negative sides of the range. The resistant genotype (R347) displayed higher genome-wide average Tajima's D (between 0.239 and 0.246, CI ¼ 90%). compared with ALM and FIG (between 0.0417 to 0.0487 and 0.0316 to 0.038, respectively, at the same confidence level) (table 1). Negative Tajima's D indicates an excess of low-frequency mutations that could be due to population expansion, background selection, or selective sweeps, whereas positive Tajima's D indicates balancing selection (the maintenance of alternative alleles that benefit the population). The top 1% and bottom 1% values of the Tajima's D, which fall in the extreme positive and negative part of the range, respectively, were used to evaluate regions departing from neutrality (supplementary table S2, Supplementary Material online). The lowest Tajima's D values range from À2.428 to À2.5016 in the three populations, and were found in the same region where the selective sweep was detected, confirming that this region of low polymorphism is indeed a region that was positively selected and which significantly departs from a neutral model of natural selection. Unexpectedly, an examination of the read alignments to the reference population SPIRL-1966 shows that the nucleotide sequence of the reference is nearly identical to that of the three other populations along the hard sweep region. The reference genome sequence is also nearly identical to that of R347 in the regions flanking the selective sweep (supplementary fig. S1, Supplementary Material online).
To further narrow down the origin of the selective sweep, the p metric was recalculated on a gene-by-gene basis on the predicted gene models across the entire NW_013535362 scaffold. Of the 190 annotated gene models, 170 had enough sequencing coverage to calculate nucleotide diversity at the set cut-offs described in the Materials and Methods section. In the ALM population, nucleotide diversity reached zero (p ¼ 0) at the protein kish-A gene (XM_013328236.1), with no SNPs detected along this coding sequence, followed by the voltage-gated sodium channel para (XM_013328250.1), with a single SNP. In the FIG population, the lowest p value (0.0001) was also found in para, followed by the gene encoding the cytochrome P450 CYP6B56 (XM_013328369.1). In the R347 population, p was zero at the gene encoding CYP6B56 (table 2). The coverage and number of SNPs for each of the genes in the scaffold are detailed in supplementary table S3, Supplementary Material online. There are approximately 15 genes in the flanking region to the left of the hard-sweep, where polymorphisms start to accumulate in all three populations, including cytochrome c oxidase, a cyclin-dependent kinase, a phosphomevalonate kinase, and Krü ppel-like-9 transcription factor, among others. To the right of the sweep, where only R347 exhibits low nucleotide diversity, we found a small conductance calciumactivated potassium channel (SkCa2), two uncharacterized proteins, and an extensin-like protein-coding gene (table 2).

Mutations in the Hard Sweep Region
We manually screened the few polymorphisms segregating along the region where all three populations have near-zero p values. Only two genes had nonsynonymous and nonconservative mutations in all or in a portion of the reads covering the position across the three populations (supplementary table S4, Supplementary Material online). One of these was the para gene (which also showed p values equal to zero), with the mutation L934F. All three populations carried the mutation in 100% of the sequenced reads covering the position, in contrast with the reference genome (supplementary fig. S2A, Supplementary Material online). This mutation corresponds to the known kdr (knock down resistance) mutation in the para gene that confers target-site resistance to DDT and pyrethroids (Kalra et al. 1967;Farnham 1977;Williamson et al. 1996). Although the SPIRL-1966 reference genome shows nearly identical nucleotide sequence along para, it does not have the kdr mutation. In the absence of population-level data for the reference genome, we confirmed that this    . R347). F ST were recalculated after eliminating the sweep region from the genome, and the resulting F ST values were almost identical to the ones obtained with the full genome, indicating that, other than the region surrounding the selective sweep, these populations are also differentiated due to geographically restricted selection that is not related to insecticide selection and could include population bottlenecks, genetic drift, or migration. The F ST method may not account for unexpected population structure and requires that each sequenced individual is assigned to a single population, an assumption that might not be valid. To verify our findings by overcoming the limitations of F ST , we utilized an approach based on principal component analysis (PCA) on the calculated major allele frequencies. With this test, we explored the 20 most-differentiated SNPs according to their FDR-corrected P values. Of those, only 2 SNPs were differentiated between the resistant and susceptible populations: one in scaffold NW_013535509.1, in an intergenic region near an ecdysteroid kinase gene cluster, and a second in the NW_013535362.

Cytochrome P450s, Gene Expression Analyses, and dN/dS Rates Ratios
To further probe the genes in the sweep region and surrounding areas, we carried out RNA-seq with midguts of ALM and R347 larvae that consumed bifenthrin or a control diet. We inspected the expression values of the (a priori selected) 42 genes that fall within the defined selective sweep. Of these, the CYP6B tandem showed the highest expression values across all the treatments (up to 15 times higher than the average FPKM in the sweep region, and up to 17 times more than the genome-wide average). In addition, a pattern consistent with elevated expression of CYP6B56 upon bifenthrin treatment in the R347 was supported with ANOVA tests (F ¼ 6.8, P 0.05) ( fig. 5). A gene coding for polyubiquitin-C also exhibited increased transcript levels upon bifenthrin exposure in the susceptible ALM population, but this difference was not significant at the set threshold ( fig. 5 and supplementary table S6, Supplementary Material online). For CYP6B54 and CYP6B55, we could not draw conclusions as the predicted gene model in the reference genome is missing a stop codon and the reads are counted jointly for both transcripts. To obtain additional confirmation of the findings of the RNA-seq analysis and to differentiate CYP6B54 from CYP6B55, we conducted a qPCR experiment with the same samples by designing unique primers specific for each of the three P450s in the cluster. The transcript coding for CYP6B54 showed constitutively higher expression in R347 compared with ALM, independent of the bifenthrin treatment (t- Material online); CYP6B55 did not show statistically significant differences between treatments or between populations and CYP6B56 failed to amplify. Because there is only a limited range of nucleotide sequence that we could use to design primers that will not cross-amplify these three closely related genes, we could not repeat the test for CY6B56. An additional transcript coding for the Krü ppel-like transcription factor 9 was added to the test, and it showed low expression across samples and treatments in agreement with the RNAseq data.
The analysis of nonsynonymous to synonymous substitutions rate ratios (x) in CYP6B-coding sequences from seven species of Lepidoptera showed that a model allowing for an increase in the rate of positive selection along the branches originating after the duplication of the most recent common  fig. S6, Supplementary Material online) explained the data better than a null hypothesis model of no change in selection pressure across the whole CYP6B phylogeny (X 2 ¼ 45.52, P corr 0.000) (same x in all branches of the tree), further supporting a role of these P450s in adaptation. We also found significant differences in selective pressure after the duplication that led to CYP6B54 and CYP6B55 (Branch No. 23) (X 2 ¼ 28.6, P corr 0.000) but not to CYP6B56

Bioassays and Analysis of Historical Data on Insecticide Usage
To determine the functional significance of the identified genomic signatures, we compared levels of resistance to bifenthrin and DDT in ALM and R347 with CPQ, a population derived from the original SPIRL-1996 line that lacks the kdr mutation. The CPQ population was the most susceptible to bifenthrin, with an LC 50 value of 0.38 ppm; in contrast, the R347 population displayed the highest level of resistance, with a LC 50 value of 24.27 ppm, confirming field observations of insecticide failure. The LC50 of ALM was 7.45, comparable to that of FIG reported by Bagchi et al. (2016). There was no difference in LC 50 for DDT between the ALM and the R347 populations, both displaying relatively high resistance to DDT (LC 50 ¼ 259.85 and 310.33 ppm, respectively) compared with CPQ (LC 50 ¼ 25.32 ppm) (table 3). These values are suggestive of a direct association between the presence of kdr and DDT resistance but not between kdr and bifenthrin resistance. Statewide in California, according to data from the California Department of Pesticide Regulation, we found that the use of bifenthrin in pounds increased 3.5-to 10.1fold, and by area 4.4-to 5.2-fold from 2007 to 2013. Applications per pound of product across that period were consistently larger in Kern County than in Madera County (supplementary fig. S7 and table S7, Supplementary Material online).

Discussion
Among the identifiable footprints of positive selection are selective sweeps. Footprints associated with evolved insecticide resistance could help understanding processes and events that ultimately shape the genomes of insects in response to human-driven selection. In this study, we utilized pooled DNA sequencing (Pool-seq) as an efficient and cost-effective method (Schlö tterer et al. 2014) to screen the genomes of three spatially distinct populations of the navel orangeworm, one of which recently evolved resistance to bifenthrin. We confirmed the existence of a large selective sweep with virtually no polymorphisms across a 0.5-Mb-long region in all these modern populations. Values of Tajima's D in this same region are below À2, a strong indication of a departure from neutrality. In addition, the regions flanking both sides of the selective sweep show very high F ST values between the resistant R347 and either of the susceptible populations. This difference could be due to differences in the point in time when the frequency of the sweep increased or when mutations accumulated after the sweep reached fixation, but they could also be due to the presence of an allele or alleles possibly linked to the difference in resistance in either or both of these regions.

Insecticide Resistance and the Selective Sweep
Reports of selective sweeps of the size identified in this study are vanishingly rare. Tian et al. (2009) reported a 1.1-Mb sweep in domesticated Zea mays, with similar evidence of multiple selection targets; this sweep likely resulted artificial selection for domestication exerted over thousands of years (Kistler et al. 2018). In humans, a >1.5-Mb-long hard sweep surrounding the lactase (LCT) gene in lactose-tolerant northern European populations (Bersaglieri et al. 2004) and in African pastoralist populations (Tishkoff et al. 2007), a textbook example of both selective sweeps and convergent evolution, is considered one of the strongest selection signatures so far reported in the human genome (Tishkoff et al. 2007). As in the case of the Z. mays selective sweep, the LCTassociated selective sweep emerged over the past $7,000 years coincident with the domestication of cows.
Within the context of insecticide selection, only a handful of studies have identified sweeps associated with resistance and, in most cases to date, the sweeps were in relatively small regions (<200 kb). Several of these studies identified genes encoding cytochrome P450 associated with the sweeps.

GBE
Genome Biol. Evol. 13(1) doi:10.1093/gbe/evaa234 Advance Access publication 4 November 2020 Schlenke and Begun (2004), reported a $100-kb selective sweep associated with a Doc transposable element inserted upstream of Cyp6g1 in Drosophila simulans, but this sweep was only marginally associated with DDT resistance. Soft, incomplete selective sweeps were found in the a-esterase gene cluster that contains the polymorphic LcaE7 gene encoding forms of the protein conferring organophosphate insecticide resistance in the Australian sheep blow fly Lucilia cuprina (Rose et al. 2011). A region containing the quantitative trait loci responsible for pyrethroid resistance in the malaria vector mosquito Anopheles funestus displayed characteristics of a selective sweep, which was then narrowed to the CYP6P9A and CYP6P9B P450 genes tandem (Barnes et al. 2017). Song et al. (2015) analyzed nine Z chromosome-linked loci in different populations of the Old World bollworm Helicoverpa armigera and detected a region possibly indicative of a selective sweep surrounding the Cyp303a1 locus. Among studies reporting selective sweeps associated with insecticide resistance, the one with findings most similar to ours is the study of Kamdem et al. (2017), who compared urban and rural populations of mosquitoes in the Anopheles gambiae complex in Cameroon, where DDT is still used for malaria control, and identified a selective sweep containing circa 80 genes inclusive of the kdr locus.

The Para Gene
The kdr locus is a single recessive point mutation resulting in an amino acid substitution from L to F in the S6 transmembrane segment of domain II of the para gene, a voltage-gated sodium channel, and the main target of both pyrethroids and DDT (Soderlund and Bloomquist 1989). The kdr locus has long been associated with resistance to DDT and to other neurotoxic insecticides, including pyrethrins and pyrethroids in many insect species (Kalra et al. 1967;Farnham 1977;Miyazaki et al. 1996;Pittendrigh et al. 1997;Bass et al. FIG. 5.-RNA-seq results showing the transcript-per-million-mapped-reads (TMM) normalized read counts in the NW_013535362.1 scaffold of the Amyelois transitella reference assembly. ALM, susceptible population; R347, resistant population; Ctrl, control diet without insecticide; Bfth, diet with bifenthrin. Black arrow in the left shows the scaffold and its direction. Dotted red line in the left shows the region of extreme low nucleotide diversity. Note that CYP6B54 and CYP6B55 are predicted in the genome as a joint transcript (missing a stop codon), for that reason, the read counts are merged for both, and we could not draw any conclusion regarding transcript levels for these two genes.
2007; Haddi et al. 2012;Dong et al. 2014). There are thus, limited ways in which para as a target of neurotoxic insecticides can evolve resistance, providing an example of "biochemically precise" (Hartley et al. 2006) convergent evolution. Similar cases of convergency were observed for the alpha-esterase coding genes that evolved resistance to organophophate insecticides in several insect species (Hartley et al. 2006;ffrench-Constant 2007), and in Na, K-ATPase coding genes in disparate insect species that evolved resistance to plant-derived cardenolides (Dobler et al. 2012).
In this study, we found that nucleotide diversity in the para gene reaches zero or near-zero values in the three A. trasnsitella populations, and the kdr mutation is present in ALM, FIG, and R347 populations and absent in the reference genome and in the CPQ population (sequenced only for the kdr locus). According to our bioassays, resistance to DDT coincides with the presence of the kdr mutation, but resistance to bifenthrin does not, indicating that other genes or loci are also responsible for the R347 resistance. It is of special interest that the kdr mutation has been linked to preadaptative resistance to pyrethroids in multiple mosquito species (Rongsriyam and Busvine 1975;Chadwick et al. 1977).

Cytochrome P450s
Assays with piperonyl butoxide, a P450 synergist, have implicated cytochrome P450s in pyrethroid detoxification in A. transitella (Demkovich et al. 2015). Of the genes in the sweep that could be involved in insecticide resistance, there is a cluster of three cytochrome P450 genes (CYP6B54, CYP6B55, and CYP6B56), where CYP6B56 has a p value of zero in R347. Although the precise substrate specificity of these A. transitella P450s has not been defined, they belong to a CYP subfamily for which there is both direct (Li et al. 2004) and indirect (Wang and Hobbs 1995;Ranasinghe and Hobbs 1998;Zhang et al. 2010;Qiu et al. 2012) evidence of involvement in pyrethroid resistance.
Tandem clusters of P450s are likely recent duplications that have persisted due to their adaptive value (Ohno 1970;Nei and Rooney 2005). The existence of P450s with the ability to detoxify phytochemicals likely served as a preadaptation for rapid evolution of resistance to insecticides, particularly those originally derived from or inspired by phytochemicals (Li et al. 2007;Omura 2013). Gene expression analyses showed increased levels of CYP6B56 transcripts in susceptible and resistant populations exposed to bifenthrin compared with control, and a comparison of nonsynonymous to synonymous rates of amino acid substitution indicates increased positive selection along the branches originating after the most recent common ancestor of the three CYP6Bs in the tandem.

Other Genes in the Sweep and Possible Regulatory Regions
Genes in the flanking region to the left of the hard-sweep, where polymorphisms start to accumulate in all three populations and where population differentiation (F ST ) is high between resistant and susceptible populations, include a cytochrome c oxidase, an enzyme that has been directly associated with pyrethroid resistance in at least one insect species (Pridgeon and Liu 2003). In addition, within this flanking region is the KL9 transcription factor, which has been implicated in regulation of P450 expression in mammals (Koh et al. 2014), and which showed segregating nonsilent mutations that differ between the resistant and both susceptible populations. To the right of the sweep, where only R347 exhibits low nucleotide diversity, we identified a gene encoding a small conductance calcium-activated potassium channel (SkCa2). Although calcium and chloride ion channels are known to interact with pyrethroids in mammals (Soderlund 2012), to date there is no evidence that insect SkCa2 channels are targeted by pyrethroids or that they are involved in resistance. In addition, PCA analyses identified two highly differentiated SNPs between susceptible and resistant populations, one of which is located just downstream of the selective sweep in the intron-exon boundary of a functionally uncharacterized locus. The other highly differentiated SNP was identified near a cluster of three ecdysteroid kinase coding genes (EcKLs). Evidence exists that EcKLs could be responsible detoxification in insects by direct or indirect phosphorylation of xenobiotics (Scanlan et al. 2020), but our transcriptomic data did not provide evidence of EcKLs constitutive transcript enrichment in the resistant line or of induction of any EcKL transcripts in either genotype upon application of bifenthrin.

Limitations
The ability of Pool-seq data to infer population allele frequencies relies heavily on sample size. We used 100 individuals from each population to ascertain our ability to infer signatures of selection. With that sample size Pool-seq outperform other techniques that rely on full genome resequencing or marker-based sequencing such as RAD-seq (Schlö tterer et al. 2014). Pool-seq, however, has some limitations, the most important of which is the loss of haplotype information, which is critical to infer additional population genomic parameters.  Bagchi et al. (2016).
In addition, because the economic impacts of A. transitella as a pest species is geographically circumscribed, the demand for the development of genomic research tools to facilitate research on this nonmodel species has been limited, so there are no available genetic maps or calculated recombination rates for this species, which limits our ability to infer the coalescence by simulations. Access to museum species would have also improved the study and facilitated the interpretations.
Based on our findings, a plausible scenario to account for the presence of the large hard sweep in contemporary field populations of A. transitella is a stacking of selective forces that prevailed at a time preceding the founding of the SPIRL-1966 colony and that persisted under continuous selection until a new reinforcing selective pressure arrived in the form of pyrethroids. DDT, which shares the same metabolic target as pyrethroids, was widely used from 1945 until 1972, when virtually all federal registrations for its use in the United States were canceled by the Environmental Protection Agency. Despite its curtailed use, DDT and its derivatives are known to remain in the environment, persisting in the soil as contaminants and moving from there across trophic levels and geographic regions (Mansouri et al. 2017). SPIRL-1966 was founded with individuals collected during an era of massive environmental exposure to DDT; although individuals from a laboratory colony in 2016 did not carry the kdr mutation, our findings are consistent with a loss of the mutation after four decades of relaxed selection under laboratory conditions (Hanai et al. 2018). Alternatively, the SPIRL-1966 never carried the kdr mutation, in which case one or more of the cytochrome P450s present in the sweep are the main candidate targets of selection and causative of the sweep. The persistence of resistance alleles (including kdr and Rdl that confer resistance to dieldrin) after decades of discontinuation of the insecticides that selected for them has been documented in mosquito species (Schechtman and Souza 2015;Grau-Bov e et al. 2020). Taken together, our findings provide an illustration of the ability of humans, in efforts to manage economically important pests, to alter the pace of evolution and the genomic composition of species over relatively short periods (Palumbi 2001).

Alignment, SNP Calling, and Identification of Genomic Regions under Positive Selection
Reads were trimmed for residual adapters and quality using Trimmomatic v. 0.32 (Bolger et al. 2014) with the following parameters: PE (for paired-end reads), ILLUMINACLIP: Truseq3_Nextera_PE.fa (adapters file), and MINLEN: 50 (minimum length to keep the read). The trimmed reads were aligned to the A. transitella reference genome (NCBI accession ASM118610v1) using bwa mem v. 0.7.12-r1039 with pairedend mode (Li 2013). The resulting SAM files were sorted and optical-and PCR-derived sequencing duplicates were marked and removed using Picard v.1.48 (Broad Institute, http:// broadinstitute.github.io/picard/, last accesed November 1, 2020). Low-quality alignments (mapping quality score <20), improper pairs, and unmapped mate reads were removed using SAMtools v.1.7 (Li et al. 2009). A pileup file was created for each separate library using the mpileup command from SAMtools. Indels were identified and separated from the main pileup files. The three pileup libraries were then subsampled for uniform SNP coverage of reads using "identify-genomicindel-regions.pl" and "subsample-pileup.pl" scripts from the Popoolation v.1.2.2 package with the following parameters:min-qual (minimum base quality) ¼ 20, -target-coverage ¼ 50, and -max-coverage (the maximum allowed coverage) ¼ 100 (Kofler, Orozco-terWengel, et al. 2011). The subsampling was carried out to minimize biases in population genomic metrics affected by sequencing errors and copy number variation that could skew coverage in affected regions.
Nucleotide diversity (p), and Tajima's D were calculated for each of the populations using 5-kb-long nonoverlapping windows with the "Variance-sliding.pl" script from Popoolation. The Tajima's D values were used to build an empirically derived distribution using 10,000 bootstrap replications with sampling with replacement with R v3.5.1. Confidence intervals of the mean (CI) were calculated on Tajima's D values. On a separate pipeline, a multiple-pileup file was created that incorporated the three trimmed and quality-filtered mapped libraries, indels were removed as described, and libraries were synchronized using "mpileup2sync.jar" from Pooplation2 software (v.1.201) (Kofler, Pandey, et al. 2011). The 3-populations pileup file was subsampled for uniform coverage and pairwise F ST values were calculated using the "fst-sliding.pl" script from Popoolation2 on a per-SNP basis.
For pooled DNA sequencing, F ST values might present bias if the sample is too small or if nonequimolar amounts of DNA for each of the individuals from each population were used. Bias can also be introduced during the PCR amplification step before sequencing (Cutler and Jensen 2010). Even though F ST calculation in Popoolation2 implements a bias correction, the estimates are still deemed biased according to Hivert et al. (2018). For that reason, to validate our F ST estimates, the mpileup synchronized file obtained from Pooplation2, was converted to a pooldata object for the "Poolfstat" v.1.0.0 R package (https://cran.r-project.org/web/packages/poolfstat). Both methods yielded similar results, but Popoolation2 does not report overall F ST between populations; accordingly, we report results from both packages. "Pcadapt" v.1.1 (Luu et al. 2017) was used to detect outlier SNPs based on PCA. In PCAdapt, z scores were calculated based on the original set of SNPs with K ¼ 2. Outliers were identified on the z scores vector using Mahalanobis distances. The distances were transformed into P values assuming a chi-squared distribution with K degrees of freedom (Luu et al. 2017).

Insecticide Bioassays
To establish the median-lethal concentrations (LC 50 ) for bifenthrin and DDT in the sequenced populations, we conducted feeding assays with semisynthetic artificial diet (Waldbauer et al. 1984). Bifenthrin (Chem Service Inc., West Chester, PA), or DDT (Sigma-Aldrich Co., St. Louis, MO) was stirred into the diet at different concentrations for each population and poured into separate 1-oz (28 ml) cups to set. Treatments and concentrations were: bifenthrin in methanol-ALM: 2 ppm, 5 ppm, 10 ppm, 12 ppm, 15 ppm, 24 ppm; DDT-ALM: 50 ppm, 100 ppm, 200 ppm, 300 ppm, 400 ppm; bifenthrin-R347: 8 ppm, 16 ppm, 24 ppm, 48 ppm, 75 ppm; DDT-R347: 50 ppm, 100 ppm, 200 ppm, 300 ppm, 400 ppm; DDT-CPQ: 10 ppm, 20 ppm, 35 ppm, 50 ppm, 75 ppm, 100 ppm. Four neonates were transferred with a soft brush into each plastic cup containing bifenthrin or methanol as the solvent control. Twenty larvae from each population were exposed to their respective bifenthrin or DDT concentrations and each assay was replicated three times per concentration. Neonate mortality on diets was assessed after 48 h and scored according to a movement response after being touched by a soft brush. Probit analysis (SPSS version 22, SPSS Inc., Chicago, IL) was applied to identify medianlethal concentrations (LC 50 ). Differences between populations were considered significant if their respective 95% confidence intervals in the Probit analysis did not overlap. We were unable to perform these assays with the FIG population because we did not establish a colony from the population used for sequencing; however, for the purposes of comparison, we incorporated the bifenthrin LC 50 established in our laboratory for this population by Bagchi et al. (2016).

RNA-Seq Analyses
We used RNA extracted from midguts from A. transitella samples from a concurrent study (in review). Briefly, fifth instar caterpillars from the ALM and the R347 populations were placed on artificial diet containing either 0.5 ppm bifenthrin or methanol as the control solvent. After 48 h, the larval midguts were dissected and flash-frozen in liquid nitrogen. There were three replicates per treatment and population combinations for a total of 12 samples, each consisting of a single midgut. Total RNA from each of the 12 samples was extracted using a Nucleospin RNA kit (Macherey-Nagel, Dü ren, Germany) according to the manufacturer's protocol. The RNA was sequenced on an Illumina HiSeq 4000 (Illumina, San Diego, CA). After sequencing, reads were preprocessed, mapped to the reference genome, and quantified for differential expression following previously described methods (Calla et al. 2017;supplementary methods, Supplementary Material online). The data were deposited in the NCBI SRA repository under accession number PRJNA548705.
Quantitative Real-Time PCR About 1 mg of each of the RNA samples that were used for RNA-seq were reverse transcribed to cDNA with a Protoscript II kit (New England Biolab, Ipswich, MA). The resulting cDNA was diluted 10-fold with nuclease-free water. qPCR reactions were assembled with SYBR green fast-qPCR mix (Thermo Fisher Scientific, Walthasm, MA) and primers specific for either CYP6B541, CYP6B55, CYP6B56, or Krü ppel-like factor (KL-9). The GADPH (Glyceraldehyde 3-phosphate dehydrogenase) was used as housekeeping gene. Primer sequences are in supplementary table S8, Supplementary Material online.

Molecular Evolution (dN/dS) Analyses
Nonsynonymous to synonymous substitution rate ratios (x) were calculated across the phylogeny of a representative set of 7 species of Lepidoptera with fully curated sets of cytochrome P450s, for a total of 22 CYP6B sequences. Nucleotide coding sequences were codon aligned using MUSCLE (Edgar 2004), and the alignment was trimmed from the ends to remove misaligned regions. A maximum-likelihood gene tree was constructed with RaXML (Stamatakis 2014), and a species tree was used to create a gene-species tree reconciliation with Notung (Chen et al. 2000). The tree with the most parsimonious topology and the codon aligned nucleotide sequences were used to run the CodeML software from PAML v.4.8 (Yang 1996(Yang , 2007, applying branch models. Detailed methods for each model, trees, and statistics are in supplementary methods, Supplementary Material online.

Sequencing the Para Locus in Laboratory Populations
DNA was extracted from ten frozen whole-body fifth instar larvae from the reference genome population SPIRL-1966, and ten midguts of fifth instar larvae of the CPQ population that fed on semisynthetic artificial diet (Waldbauer et al. 1984). Extractions were carried out using an E.Z.N.A. insect DNA kit (Omega Bio-Tek, Norcross, GA) according to manufacturer's instructions. Primers were designed using Primer3 (Untergasser et al. 2012) as implemented in Geneious software (v 9.1.8) (Kearse et al. 2012) to sequence the region of the para gene containing the kdr mutation (Forward 5 0 -ACCAAGGTGGAACTTCACAGAT-3 0 Reverse 5 0 -AGCAATTTCAAGAAGTCAGCAACA-3 0 ). Amplicons were Sanger sequenced and trace analysis and alignments were performed in Geneious.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.