Genome evolution in an agricultural pest following adoption of transgenic crops.

Replacing synthetic insecticides with transgenic crops for pest management has been economically and environmentally beneficial, but these benefits erode as pests evolve resistance. It has been proposed that novel genomic approaches could track molecular signals of emerging resistance to aid in resistance management. To test this, we quantified patterns of genomic change in Helicoverpa zea, a major lepidopteran pest and target of transgenic Bacillus thuringiensis (Bt) crops, between 2002 and 2017 as both Bt crop adoption and resistance increased in North America. Genomic scans of wild H. zea were paired with quantitative trait locus (QTL) analyses and showed the genomic architecture of field-evolved Cry1Ab resistance was polygenic, likely arising from standing genetic variation. Resistance to pyramided Cry1A.105 and Cry2Ab2 toxins was controlled by fewer loci. Of the 11 previously described Bt resistance genes, 9 showed no significant change over time or major effects on resistance. We were unable to rule out a contribution of aminopeptidases (apns), as a cluster of apn genes were found within a Cry-associated QTL. Molecular signals of emerging Bt resistance were detectable as early as 2012 in our samples, and we discuss the potential and pitfalls of whole-genome analysis for resistance monitoring based on our findings. This first study of Bt resistance evolution using whole-genome analysis of field-collected specimens demonstrates the need for a more holistic approach to examining rapid adaptation to novel selection pressures in agricultural ecosystems.

The study of rapid evolutionary responses to natural selection has gained momentum in the past decade, and along with it recognition of the relevance of pest adaptation to control measures (1). Melander published the first observation of insecticide resistance evolution in 1914 (2), and since then pests and weeds have adapted to myriad chemical, biological, and cultural control measures (3–6). While there has been an emphasis on finding alleles at single genetic loci that confer resistance to pesticides, recognition that surviving insecticide or herbicide applications could involve a suite of genomic changes has increased (7). A clear understanding of the timing and complexity of changes associated with adaptation to agricultural ecosystems is lacking, however (8).

The bacterium Bacillus thuringiensis (Bt) produces one or more crystalline (Cry) proteins that are each toxic to a narrow range of insect pests when ingested. Genes encoding Cry toxins were engineered into cultivated crops to reduce pest damage, producing so-called Bt crops. This major innovation in agricultural pest management has been adopted globally (9). The economic and environmental benefits of this technology are significant. They can replace synthetic insecticides, which have damaging off-target effects on pollinators and natural enemies, and contribute to area-wide pest suppression in non-Bt crops (9–11). Many pest species have evolved resistance to one or more Cry toxins, and there is concern about the durability of these transgenic crops (12). Resistance management strategies have been implemented effectively for a few major pest species with high susceptibility to the toxins (9). Inadequacy of these strategies for species with only moderate susceptibility to Cry toxins has provided opportunities for resistance evolution (12).

Discovery of resistance-conferring genes has been of long-standing interest, in part due to expectations that this knowledge could lead to molecular diagnostic approaches that facilitate detection of emerging resistance in pests (13, 14). Molecular signals of resistance could then serve as triggers for regulatory actions that would minimize resistance risk and preserve environmentally sustainable technologies, like Cry toxins (13). Genes discovered to have major impacts on Cry resistance are used to a limited extent for molecular monitoring (15–18). To date, there are at least 11 such genes in the Lepidoptera (19, 20), many of which were discovered in laboratory-selected populations and from screening wild populations for single genes of major effect (16–18, 21). It is unclear whether molecular monitoring of these genes over time would effectively track resistance evolution in wild populations.

Helicoverpa zea is a migratory pest species (22) which comprises a single panmictic population in North America (23). It has been controlled with Bt crops since 1996, when corn expressing Cry1Ab and cotton expressing Cry1Ac became commercially available (24). Corn expressing Cry1F, a toxin with lower protein structure similarity to Cry1Ab and Cry1Ac, was commercialized in 1998 (24), and in 2003 corn cultivars expressing multiple toxins (i.e., pyramids) were introduced. One pyramided corn cultivar relevant to this work expresses a hybrid Cry1Ac and Cry1F toxin called Cry1A.105, as well as Cry2Ab2. Over the past 20 y, North American H. zea have been exposed to all of these toxins with varying protein structure similarity (24) and which were initially protective to the plants that expressed them (25–27).

North American adoption rates of Bt crops grew from 8 and 15% of corn and cotton acreage in 1997 to 82 and 88% by 2020. This rapid expansion placed significant selection pressure on target pests, including H. zea, which showed widespread, damaging levels of resistance by 2016 (27–30). Existing Bt corn and cotton cultivars do not produce a high dose of Cry toxin relative to what is required to kill H. zea larvae that are heterozygous for resistance alleles (31). This likely contributed to their recent resistant status (12). Under these conditions, resistance to Cry toxins was predicted to be polygenic (32, 33) and to arise from standing genetic variation (34). Furthermore, as resistance evolves to one toxin, cross-resistance can occur to others, suggesting that some genetic mechanisms of resistance are shared between toxins. H. zea selected for resistance to one toxin are more likely to have heightened resistance to others if the amino acid sequences of the toxins are similar (35).

H. zea is an ideal pest for examining the timing and complexity of resistance evolution due to its high levels of gene flow, recently described Bt resistance in wild populations, and previous predictions of polygenic resistance. Between 2002 and 2017, we saved ca. 1,000 adult moths in ethanol each year, and we used these collections to track genomic changes in H. zea over time. Instead of simply monitoring for changes in previously identified candidate genes, we used a gene-agnostic approach to identify genomic regions associated with H. zea Cry resistance. Our study aimed to 1) identify contiguous windows of the H. zea genome with strong temporal allele frequency changes as an indicator of field-evolved Bt resistance, 2) separately determine the genomic architecture of Cry1Ab and pyramided Cry1A.105+Cry2Ab2 resistance by quantitative trait locus (QTL) analysis, and 3) compare these two datasets to explore the timing and complexity of resistance evolution. None of the 11 previously described Cry resistance candidate genes showed significant allele frequency changes when whole-genome sequenced H. zea collected in 2002 and 2017 were compared. Our QTL analysis also revealed no major role for nine of those candidate genes. However, a Cry-associated QTL contained a cluster of aminopeptidases, a gene family known for involvement in Bt resistance. We found that the genomic architecture of Cry1Ab resistance was polygenic and Cry1Ab resistance QTL were identified on several chromosomes. A major QTL spanning most of one chromosome strongly and significantly impacted weight gain in larvae exposed to the Cry1A.105+Cry2Ab2 pyramid. Cry-associated QTL directly overlapped with only a subset of genomic regions that underwent significant temporal allele frequency change, highlighting both potential and pitfalls of genomic monitoring for detection of emerging resistance.

Results

Superscaffolding of the H. zea Genome.

The H. zea genome [v. 1.0, GCA_002150865.1 (36)] consisted of 2,975 scaffolds, making it difficult to identify contiguous genomic regions. We superscaffolded this fragmented assembly by aligning it to the more contiguous Helicoverpa armigera genome (v. 1.0, GCF_002156985.1). Resulting superscaffolds were organized into chromosomes using a linkage map [Rqtl v. 1.47-9 (37)] generated from one of the F2 genetic crosses used for QTL analysis. The final assembly contained 2,732 of the original 2,975 H. zea scaffolds and consisted of 31 linkage groups (LGs; equivalent to the 31 H. zea chromosomes) plus 15 unplaced superscaffolds (Dataset S1). The first 28 H. zea LGs were named according to their homologous Bombyx mori chromosome (v. 1.0, GCF_014905235.1).

Temporal Allele Frequency Shifts Occurred throughout the Genomes of Wild H. zea.

We initially studied the changes that occurred in the years following Bt crop adoption, by sparsely sampling and analyzing the genomes of H. zea collected in 2002, 2007, 2012, and 2016. Archived males (n = 265) collected from pheromone-baited traps (38) in Bossier Parish, LA (), were used to generate double-digest restriction site-associated DNA sequencing (ddRADseq) libraries and sequenced according to Fritz et al. (39, 40). Following filtering and genome alignment, we identified 14,398 single nucleotide polymorphisms (SNPs) from 259 individuals (). Nucleotide diversity (π) and heterozygosity (F) did not change between years, indicating there had been no strong population contraction (). Overall genomic divergence was low according to Weir and Cockerham’s FST (FST < 0.004; ), and this low divergence was also supported by a sequential k-means analysis that found our archived H. zea were best described as a single genetic cluster ().

Genetic divergence between 2012 and 2016 was three times higher (FST = 0.0025) than for any other sequentially sampled period (). We further explored this accelerated temporal genetic change with a discriminant analysis of principal components (DAPC). This identified SNP combinations that clustered individuals into groups (k = 2 to 4) by maximizing between-group and minimizing within-group variation. Using these SNP combinations, we estimated posterior membership probabilities for each individual. When k was equal to 2, samples from 2002, 2007, and 2012 were assigned to both genetic clusters. Samples collected in 2016 were strongly biased toward a single genetic cluster (Fig. 1), and this trend persisted regardless of k-value (). During this period (2010 to 2016), the mean area consumed per ear of Cry1A.105+Cry2Ab2-expressing sweet corn increased from near 0 to ca. 3 cm2, and the proportion of two-toxin-expressing ears containing late instars grew from near 0 to 75% (27).

Fig. 1.

Posterior membership probabilities for H. zea individuals collected in Bossier Parish, LA, in 2002, 2007, 2012, and 2016. Genotypic clustering used ddRAD SNPs and assumed a prior number of clusters (k) equal to 2.

To identify SNPs with the greatest influence on this genomic shift between 2012 and 2016, we removed the top 1%, 2.5%, and 5% of variance-contributing SNPs from our dataset and reevaluated population membership probability. Removal of 5% (n = 720 SNPs) resulted in nearly complete assignment of all individuals to a single genetic cluster, indicating that all contributed to the observed 2012-to-2016 genomic shift (). We used a second approach, FST outlier analysis (41), to identify genomic regions that underwent significant temporal divergence and detected 53 (of 14,398) ddRAD SNPs as having greater than expected change (FST =0.068 to 0.278; ). Twenty-one of those 53 were among the top variance contributing SNPs from our DAPC. This conservative set of outlier SNPs, generated from both FST outlier and DAPC analyses, were found on 12 LGs: 1, 5, 8, 9, 10, 11, 13, 15, 17, 22, 24, and 28. Altogether, this demonstrated that among-year population genetic divergence was low for H. zea from LA, but allele frequency changes mildly increased after 2012. The SNPs that contributed most to these temporal changes were not concentrated in a single genomic region but spread across multiple chromosomes.

Whole-Genome Analysis of Temporal Allele Frequency Change in H. zea.

To further examine temporal changes, we sampled entire genomes of a smaller number of wild H. zea males from Bossier Parish, LA, using Illumina whole-genome sequencing (WGS). Thirteen males collected in 2002, along with 11 males each from 2012 and 2017, were sequenced and produced a genome-aligned, filtered dataset containing 4,999,128 SNPs (). Using this SNP dataset, we conducted pairwise analyses of genomic divergence between populations collected from 2002, 2012, and 2017 using 5-, 10-, 20-, and 40-kb sliding window-averaged Weir and Cockerham’s FST. We report the 10-kb sliding window analysis because it captured the major genomic changes that occurred for all sliding window sizes, but all analyses are in . Our empirical threshold for statistical significance in a pairwise comparison was 6 SDs (ZFST > 6) from the mean FST value in a 10-kb window (42). Because both genetic drift and selection influence the landscape of genomic divergence between populations, we quantified the frequency with which our empirical threshold, calculated from the comparison of 2002-to-2017 samples (FST = 0.078), could be reached due to drift alone using genotypic simulations (43). These simulations demonstrated the probability was very low (P = 0.0180 to 0.0003) and negatively correlated with present-day population size ().

In total, 69 H. zea scaffolds on 30 LGs or unplaced superscaffolds contained 10-kb windows with greater-than-expected allele frequency change in the pairwise comparison of 2002 and 2017 samples. For comparisons from 2002 to 2012 and 2012 to 2017, 91 scaffolds and 79 scaffolds, respectively, each on 31 LGs or unplaced superscaffolds, contained 10-kb windows with greater than expected change. The greatest changes between 2002 and 2017 occurred on LG13 at four linked scaffolds (scaffolds 173, 569, 1612, and 759), where windowed-FST values greatly exceeded those from any other region of the genome (scaffold 1612; max FST = 0.292). LG9 showed the second highest degree of change (scaffold 332; max FST = 0.198), followed by LGs including but not limited to 14, 19, and 23 (). Adjacent windows were merged together if their FST values were significant, and we identified 11 windows, including one on LG13, with elevated FST in both 2002-to-2012 and 2012-to-2017 comparisons. Additional LGs containing windows of elevated FST in multiple by-year comparisons were 4, 5, 14, 21, 28, 31, and superscaffold 44. Temporal changes at LG9, our second-highest FST peak, were complex and found within the same superscaffold (NW_018395585.1 and NW_018395392.1), but no windows were shared among by-year comparisons. Regions of genomic divergence were compared with our ddRAD outliers, and 3 (of 21) SNP outliers were within or adjacent to windows of elevated FST. Two were within windows on LG13, and one was on LG9 within 20 kb of a window (). This supported the signatures of temporal change at the two largest peaks in our WGS analysis, but with a much larger sample size.

Allele Frequency Changes at Candidate Bt Resistance Genes.

To determine whether these significant allele frequency changes occurred at the 11 previously reported candidate Bt resistance genes (Table 1), we examined 10-kb sliding window-averaged FST values for the length of each gene using specimens collected in 2002 and 2017. FST ranged from −0.018 to 0.052 (Table 1 and ), with the strongest changes occurring at apn4. No 10-kb window exceeded our genome-wide empirical threshold for significant divergence, however. These modest allele frequency changes indicated that de novo mutations in the cis-regulatory or coding sequences of these genes were unlikely to be the cause of field-evolved Bt resistance in H. zea. Two additional candidate genes for resistance to synthetic pyrethroid insecticides, the voltage-gated sodium channel gene (vgsc) (53) and cyp337b3 (54), were also examined for evidence of allele frequency change, but none was found.

Table 1.

Average FST values and additive effect sizes of SNPs near Bt resistance candidate genes

Gene	Ref.	Average FST (no. of SNPs)	Cry1Ab effect size	Cry1A.105 + Cry2Ab2 effect size	Providence control effect size	Obsession control effect size
alp	44	−0.011 (171)	−10.4	2.6	−14.0	1.2
apn1	45	0.014 (152)	18.1*	29.0*	16.2	14.7
apn4	46	0.052 (188)	18.1*	29.0*	16.2	14.7
abcA2	47	0.011 (155)	−3.5	11.1	−33.6*	20.6
abcC2	48	−0.005 (100)	−7.8	3.1	−4.8	7.4
white	49	0.014 (197)	11.2	−3.8	32.0*	−15.2
calp	16	0.012 (221)	2.4	NA	−18.2	NA
cad2	50	−0.001 (118)	5.3	2.1	−10.9	0.7
cad-86C	51	−0.008 (1079)	6.7	4.8	5.5	3.2
map4K4	52	−0.018 (65)	−8.4	−0.3	−4.4	1.8
tspan1	18	−0.007 (165)	11.4	−3.9	34.3*	−15.0

Additive effect sizes are the average value of β (LMM) for SNPs within 250 kb of either side of each gene. They represent differences from the mean 7 d weight (milligrams) associated with a single copy of a field parent allele. Averages were calculated from between 1 and 16 SNPs. provides their UniProt/GenBank IDs and genome positions for each gene. Field-collected WGS samples from 2002 and 2017 were used to calculate FST values across the number of SNPs in parentheses. Asterisks denote that the candidate gene is found within a QTL window. NA indicates that no SNPs were within 250 kb of the candidate gene.

QTL Analysis of Cry 1Ab, Cry1A.105, and Cry2Ab2 Resistance.

We conducted a QTL analysis to identify genomic regions that strongly contributed to H. zea Cry resistance. The average weight of field-collected Cry-resistant larvae was 5.4 times greater than that of the laboratory-reared Cry-susceptible larvae after 7 d on a diet containing a diagnostic dose of lyophilized Cry1Ab-expressing leaf tissue and 13.8 times greater on diet containing a diagnostic dose of Cry1A.105+Cry2Ab2-expressing leaf tissue (). Differences between resistant and susceptible populations were significant according to a Mann–Whitney U test (Cry1Ab: df = 1, W = 189, P < 0.001; Cry1A.105+Cry2Ab2: df = 1, W = 38.5, P < 0.001). First instar F2 offspring from each of two genetic crosses between resistant and susceptible grandparents (n for F2s = 153 and 160) were split between bioassays on a diagnostic dose of Cry-expressing leaf tissue and a control containing nonexpressing leaf tissue. The controls were used to identify QTL for larval growth unrelated to Cry resistance. The family exposed to Cry1Ab-expressing leaf tissue had a mean 7-d larval weight of 77.7 mg (SEM = 5.3) compared to 173.1 mg (SEM = 10.6) on the untreated control (W = 4805, P < 0.001). Mean 7-d weights of the family exposed to Cry1A.105+Cry2Ab2-expressing leaf tissue were 61.1 mg (SEM = 4.3) and 166.5 mg (SEM = 7.2) on treated diet and the untreated control (W = 5869, P < 0.001). Weight distributions of F2 offspring on Cry-treated diets were continuous, and few Cry-exposed larvae reached the resistant grandparent’s weight, suggesting that resistance was polygenic ().

Two hundred seventy-seven F2 offspring from the four treatments were used to generate ddRAD libraries with their parents and grandparents, and all samples were sequenced to an average coverage depth of 76.7× per locus. Following read filtering and genome alignment, we applied a Bayesian sparse linear mixed model [BSLMM; GEMMA v.0.98.4 (55, 56)] to 6,756 to 6,268 filtered SNPs to identify the genomic architecture of Cry resistance. For F2s exposed to Cry1Ab, the proportion of phenotypic variance explained by all ddRAD SNPs was 0.345 [0.035, 0.711], and less than half of this genetic variance could be explained by SNPs of major effect (genetic variance explained by sparse effects [PGE] = 0.415 [0.000, 0.964], n SNPs = 57.61 [0, 270]). For F2s exposed to Cry1A.105+Cry2Ab2, a greater proportion of variance in larval growth could be explained by our full SNP dataset (phenotypic variance explained [PVE] = 0.576 [0.363, 0.783]), and SNPs of large effect size explained more of the total genetic variance (PGE = 0.770 [0.315, 0.989], n SNPs = 7.87 [1, 58]). Results from these single-family crosses indicated that H. zea resistance to Cry1Ab was polygenic, but fewer loci controlled resistance to Cry1A.105+Cry2Ab2.

We identified windows of SNPs throughout the H. zea genome with statistically significant additive effects on 7-d larval weight gain using a univariate linear mixed model [LMM; GEMMA v.0.98.4 (55)]. QTL windows consisted of one or more neighboring SNPs with likelihood ratio test (LRT) P values < 0.01. Significant SNPs within a window were never interrupted by more than 20 nonassociated SNPs, and we required that each QTL window contain at least one SNP with a posterior inclusion probability >0.01 in the BSLMM to make our analysis more conservative. Fourteen windows on eight LGs (3, 7, 9, 10, 15, 18, 21, 24) and unplaced superscaffold 41 influenced 7-d larval weights on the Cry1Ab diet (). Of these, alleles from the resistant grandparent had positive effects on growth in all QTL windows on LGs 9 (number of windows = 2), 10 (n = 3), and 18 (n = 3), suggesting a role in Cry resistance for each LG. Yet, no QTL windows were significant when we lowered our LRT significance threshold (α = 1 × 10−6), as would be expected if Cry1Ab resistance resulted from many genes of small effect size. For Cry1A.105+Cry2Ab2, we identified one QTL window spanning most of LG9 (), and this remained the same at our more conservative LRT significance threshold (α = 1 × 10−6).

Detecting QTL of small effect requires large sample sizes (57), and our power to detect such QTL was unclear, particularly for the more polygenic Cry1Ab resistance. We conducted multiple power analyses because our true “alpha” value fell within the range of 0.01 and 1 × 10−6 due to our two test criteria. These analyses also allowed for different levels of polygenic contribution to trait variance (57). We estimated that our family sizes allowed us to detect QTL explaining a minimum of 10 to 30% of the trait variance, depending upon α value (0.01 or 1 × 10−6) and level of polygenic contribution to trait variance ().

We attempted to separate QTL which provided a general growth advantage from those contributing to Cry resistance by comparing QTL for growth on Cry-expressing leaf tissue with growth on the nonexpressing control for each family ( and Datasets S2–S5). QTL associated with growth on diet containing Cry1A.105+Cry2Ab2 never overlapped with those for its control diet. Two QTL windows were associated with growth on both diet containing Cry1Ab-expressing leaf tissue and its control: one on LG9 (scaffold 1124) and another on LG10 (scaffold 143). For both LGs, resistant grandparent alleles significantly increased 7-d larval weight on Cry1Ab-expressing and nonexpressing leaf tissue, indicating they provided a general growth advantage. Resistant grandparent alleles within the large Cry1A.105+Cry2Ab2 QTL on LG9 also increased larval weight after feeding on toxin-treated diet, and this QTL overlapped with the LG9 growth-associated QTLs in the presence of both Cry1Ab and its control. The breadth of the Cry1A.105+Cry2Ab2 QTL was greater, and its effect size was nearly twice that of the Cry1Ab-associated QTL, however. We further analyzed the proportion of PVE by all SNPs on LG9 (GEMMA, BSLMM) because a similar amount of weight gain represented different proportions of total body mass in the toxin-containing and control treatments. Our estimates for LG9 PVE were 0.31 [0.05, 0.73] for Cry1Ab, 0.43 [0.23, 0.70] for Cry1A.105+Cry2Ab2 and 0.18 [0.01, 0.58] and 0.09 [0.00, 0.36] for their respective controls. Altogether, this indicated that a portion of LG9 conveys a general growth advantage to H. zea and enhances growth in the presence of Cry1Ab, Cry1A.105, and Cry2Ab2. Yet, most of LG9 is strongly associated with Cry1A.105+Cry2Ab2 resistance.

Genome positions of candidate Bt resistance genes were compared with QTL positions, and the average additive effect sizes of SNPs near candidate genes were calculated, providing insight into whether any were likely to play a role in H. zea Cry tolerance (Table 1). All apn genes, including apn1 and apn4, were within Cry1Ab and Cry1A.105+Cry2Ab2 QTL, but not within the general growth-associated QTL on LG9. A field parent allele at the SNPs near apn1 and apn4 on average increased offspring weight by 18.1 mg for the Cry1Ab treatment, 29.0 mg for the Cry1A.105+Cry2Ab2 treatment, and 16.2 mg and 14.7 mg for their respective controls. Yet, some offspring that were homozygous for the susceptible grandparent allele at the SNP nearest apn1 and apn4 (within 20 kb) grew to be some of the largest larvae after feeding on toxin-containing diet (). While we could not rule out the possibility that apn1 and apn4 play a role in H. zea Cry resistance, wild resistance alleles at these genes were not necessary for larval weight gain on either Cry-treated diet. No other candidate genes were within Cry-associated QTL, although abcA2 was within 100 kb of a negative growth-associated QTL (LG17) on untreated diet. Genes white and tspan1 were within a large, positive growth-associated QTL (LG10) on untreated diet. As described above, this overlapped with a QTL for Cry1Ab, but this Cry1Ab QTL consisted of a single SNP which was more than 2Mb from white and tspan1, suggesting they are not involved in Cry1Ab resistance.

Resistance-Associated QTL Overlap with a Subset of FST Outliers.

Windows of temporal genomic divergence that do not overlap with Cry-associated QTL would refute a strong contribution of cis-regulatory or coding sequence mutations in these regions to resistance. This was true of many regions with elevated FST in our wild H. zea. Yet, overlapping Cry-associated QTL with significant temporal genomic change in wild H. zea should support involvement of those regions in resistance. We compared the results of our FST and QTL analyses and identified windows that overlapped by >1 bp. LGs 9 and 18 contained regions of strong temporal divergence in wild H. zea, as well as Cry-associated QTL (Fig. 2). Many windows of the H. zea genome showed elevated levels of divergence, however, leading us to quantify the probability that overlaps between the QTL and FST outlier windows could be expected by chance alone. We conducted a permutation test, which fixed the QTL positions on our superscaffolded genome, then randomly shuffled the positions of our FST outlier windows 10,000 times and counted their overlaps with fixed QTL (58). Observed overlap between Cry-associated QTL and FST outlier windows was greater than the average numbers of overlaps due to chance alone, but this difference was only significant for the Cry1A.105+Cry2Ab2 treatment (), likely due to the breadth of the QTL. Some of the FST outlier windows in the Cry1A.105+Cry2Ab2 may contain protein-coding or cis-regulatory mutations linked to field-evolved resistance, but we were unable to demonstrate that FST outlier windows which overlapped with Cry1Ab QTL were not due to chance alone. For the control treatments, overlap was not significantly greater than would be expected by chance ().

Fig. 2.

Genomic divergence in Cry-associated regions of the H. zea genome. SNP additive effect sizes, β (LMM), of the resistant parent allele are plotted for Cry1Ab (A) and Cry1A.105+Cry2Ab2 (B). All SNPs in significant QTL windows are in color. Genome-wide divergence in field-collected H. zea from Louisiana (2002-2017) are in C. Alternating light and dark gray points represent pairwise FST values for 10-kb windows with a 1-kb step. Points above the red line underwent significant temporal genomic divergence. Points in color were also associated with increased growth on Cry1Ab (red), Cry1A.105+Cry2Ab2 (blue), or both (purple) in our QTL analysis.

Notably, LG13 showed the greatest allele frequency changes between 2002 and 2017, yet it did not contribute to growth differences on either Cry1Ab- or Cry1A.105+Cry2Ab2-containing diet (Fig. 2 and ). These strong changes on LG13 provided us the opportunity to compare them with the strength of selection observed in Cry-associated regions of LGs 9 and 18. Linkage disequilibrium (LD) extended for 3 kb along most of the genome, including on Cry-associated scaffolds (). Correlations among alleles on LG13 (scaffolds 1612 and 569) were high relative to the genome-wide average in all years, however, and increased over time. This buildup of LD over time was not observed on Cry1Ab- or Cry1A.105+Cry2Ab2-associated scaffolds.

The magnitude of the SNP frequency changes between 2002 and 2017 was also greater on LG13 scaffolds 569 and 1612 compared to Cry-associated scaffolds. The maximum change was 0.69 for both scaffolds on LG13 but only 0.59 for Cry-associated scaffolds. We quantified changes in haplotype frequency within windows of elevated FST for Cry-associated LGs 9 and 18, as well as LG13. For Cry1Ab-associated regions of LG18, the dominant haplotype in 2017 was at higher frequency in 2002 and declined over time (Fig. 3). This contrasted with the dominant haplotypes on LG9 near both Cry1Ab and Cry1A.105+Cry2Ab2 QTL, whose frequencies were lowest in 2002 then increased (Fig. 3). By 2017, haplotypes at all Cry-associated regions were at intermediate frequencies: 0.38 to 0.46. Like LG9, haplotype frequencies at LG13 were low in 2002, but unlike for LG9, haplotypes at LG13 underwent strong and sustained directional change over time (Fig. 3).

Fig. 3.

Haplotype frequency changes at two Cry-associated LGs (9 and 18) and LG13, which was not associated with resistance to Cry toxins. Haplotype frequencies in red were associated with Cry1Ab resistance, blue with Cry1A.105+Cry2Ab2, and purple with both Cry1Ab and Cry1A.105+Cry2Ab2.

Such strong directional selection on LG13 haplotypes led us to examine H. zea collected from our Louisiana trapping sites in 1998, to determine whether allele frequencies were even lower in that year. Moths from 1998 had DNA quality that was too low for use in ddRADseq and WGS analyses, but short fragments could be PCR-amplified. Two genes on LG13, one carboxypeptidase q (cpq1) and one cytochrome p450 (cyp333b3), contained nonsynonymous mutations () and fell within regions of significant temporal divergence in our genome scan. Sanger sequencing of ca. 200-bp fragments spanning one nonsynonymous mutation in each of cpq1 and cyp333b3 indicated their frequencies in 1998 were 0.07 and 0.1, respectively. These same mutations had frequencies of 0.27 and 0.35 in 2002 and 0.86 and 0.82 in 2017 (). Such strong directional change and buildup of LD is consistent with a hard selective sweep in this region of the genome (59), but likely in response to some factor other than the Cry toxins studied here. Cry1Ab- and Cry1A.105+Cry2Ab2-associated regions of the genome did not show these same patterns, however.

Contemporary Patterns of Cry-Associated Selection in H. zea.

Ongoing selection in Cry-associated regions of the genome were examined by collecting fifth and sixth instar H. zea from ears of two sweet corn varieties grown in Central Maryland Research and Education Center (CMREC) plots in 2017: Cry1A.105+Cry2Ab2 expressing Obsession II (n = 8) and its non-Bt isoline, Obsession I (n = 8). The genomes of these individuals were sequenced with Illumina short reads. H. zea populations are considered panmictic due to their large effective population sizes and migratory behavior (60). We reasoned that the same LGs under selection by Cry1Ab, Cry1A.105, and Cry2Ab2 in H. zea from LA, would be involved in resistance in MD. We confirmed panmixia by comparing 2017 LA and MD samples using a filtered dataset of 14,697 SNPs. The pairwise genome-wide FST value was 0.0004, and the 95% confidence intervals around this value (−0.0002, 0.0008) overlapped with zero, demonstrating little divergence between populations. We then used our full SNP dataset to examine 10-kb window-averaged FST values, which identified differences in allele frequencies for larvae that survived on Cry1A.105+Cry2Ab2-expressing sweet corn relative to its non-Bt isoline. No differences between larvae collected from Cry-expressing and nonexpressing corn were observed at LG9, but significant divergence was found in the Cry1Ab QTL region on LG18 (). The Cry1Ab QTL on LG18 was very large (), and regions of divergence between MD field collections did not overlap with regions of temporal change in LA between 2002 and 2017. Potential hypotheses to explain this lack of overlap are below.

Discussion

We paired QTL analyses of Cry resistance with a whole-genome analysis of wild H. zea specimens collected over time and space to quantify resistance-associated changes that occurred following adoption of Bt crops in the United States. Our results supported previous predictions that resistance to Cry1Ab and, to a lesser extent, pyramided Cry1A.105+Cry2Ab2, was the result of polygenic adaptation from standing genetic variation. Substantial allele frequency changes occurred throughout the genomes of wild H. zea in concert with increasing resistance to Cry toxins, but none occurred at previously identified lepidopteran Cry resistance genes. While candidate genes apn1 and apn4 showed some potential for involvement in Cry resistance according to QTL analysis, our QTL results also indicated the possibility of a general growth advantage from wild H. zea alleles in a nearby region of the same chromosome. Combined QTL and whole genome analyses identified several major Cry-associated changes that occurred in novel genomic regions in wild H. zea over our 15-y sampling period, shedding light on the evolution of resistance in this pest. Results from our study also offer insight into how, when, and whether genomic approaches could be used for future resistance risk monitoring.

For four decades there has been an emphasis on identifying single genes that confer resistance to control measures in insect pests and weeds, and we now have a catalog of such genes (5, 61–63). There are cases where a single (64) or multigene diagnostic targeting previously established candidate genes (65) have been useful for resistance monitoring. Yet, recently there has been concern that these focused approaches could mislead, as resistance evolves due to alleles already at low to intermediate frequencies in unexpected (63) or even multiple genes (5, 8, 66). Our work demonstrates that these concerns are appropriate. Aminopeptidases, for example, have been identified as functional receptors of Cry1A toxins (67). Temporal genomic divergence at apn4 was highest of all candidate genes, making it the most promising for resistance diagnostics, but allele frequency changes at this gene were not statistically significant. Alleles near apn4 from field-collected resistant individuals provided a modest general growth advantage and a Cry tolerance advantage in our QTL analysis, yet some of the largest larvae in our growth assays were homozygous for the susceptible apn4 allele. These results indicate that use of candidate genes, even apn4 which showed strongest potential for involvement in Cry resistance, can be problematic. In our case, resistance to Cry1Ab and Cry1A.105+Cry2Ab2 pyramids may be aided by, but did not require, apn4, likely because it (or a neighboring apn) is one of many genes that impact growth on Cry toxins. Solely focusing on this most promising candidate gene for resistance diagnostics would have missed the coming problem.

Instead, our combined QTL analyses and genome scans identified novel regions of the H. zea genome likely under selection by Cry toxins expressed in Bt corn and cotton. According to our GEMMA analysis, Cry1Ab resistance was polygenic, where SNPs of major effect explained less than half of the genetic variance contributing to larval growth (PGE = 0.415). Our sample size was sufficient to detect QTL explaining > 10 to 30% of the variance in larval weight, and we identified multiple regions of LGs 9, 10, and 18 where alleles from the field-collected, resistant grandparent increased larval growth (Fig. 2). Yet, it is likely that many Cry1Ab resistance-conferring loci of minor effect remain to be discovered. Resistant grandparent alleles in some LG9 and 10 QTLs also increased larval growth on a diet containing leaf tissue without Cry toxins, indicating they provide a general growth advantage and have a low fitness cost in the absence of Cry toxins.

When we examined temporal allele frequency changes in wild H. zea within Cry1Ab QTL we were unable to show that more FST outliers occurred within these QTL than expected by chance alone. This highlights one major challenge for molecular resistance monitoring. When traits are polygenic and arise from extensive standing variation in nature, selection can act upon many variants, each of small effect size (68). All of the minor effect variants segregating in natural populations may not be represented by a single mapping family, like the one used here. Therefore, some of the temporal allele frequency changes observed in wild H. zea may represent changes at loci of small effect that were not detectable in our QTL analysis. Separating selection events that occur at minor-effect resistance loci from those caused by other environmental factors represents a challenge for molecular resistance monitoring, due to the large numbers of samples required to link loci of minor effect to a trait.

QTL results for larval growth on diet containing Cry1A.105+Cry2Ab2 revealed that most of LG9 strongly impacted larval growth. Due to its breadth, this QTL encompassed the Cry1Ab and control growth-associated regions of LG9, but portions were uniquely associated with growth on Cry1A.105+Cry2Ab2. One Cry1A.105+Cry2Ab2-associated region corresponded to the second greatest temporal change observed in field-collected samples (Fig. 2) and may be uniquely linked to Cry2Ab2 resistance. Previous studies of resistance mechanisms support this conclusion. Cry1A.105 competes with Cry1Ab and Cry1Ac for the same lepidopteran midgut receptors (69), and selection for resistance to one of these toxins results in cross-resistance for the others (35). Yet, major cross-resistance to Cry1A and Cry2A toxins is uncommon (70), including for H. zea (71), and is likely due to their distinct larval midgut receptors (72–74). This difference supports a genetic mechanism for Cry2Ab2 resistance that is distinct from Cry1A. Although we found a single, large QTL for Cry1A.105+Cry2Ab2 resistance, Yang et al. (75) found that phenotypic variance for Cry2Ab2 resistance fit a polygenic model. Perhaps loci of small effect size on other chromosomes were not detectable in our QTL analysis () and further study would reveal those loci.

Our data begin to shed light on the time course of field-evolved Cry resistance evolution in H. zea. Modest changes on a superscaffold in a Cry1Ab-associated region of LG18 (NW_018395428.1) and a Cry1A.105+Cry2Ab2-associated region of LG9 (NW_018395392.1) occurred between 2002 and 2012 (). These genomic changes were detectable as early as 2012, when Dively et al. (27) was detecting damaging levels of phenotypic resistance to Cry1Ab-expressing sweet corn and growing damage to Cry1A.105+Cry2Ab2-expressing sweet corn. This offers some hope that genomic monitoring could be useful for predicting emerging resistance, but future work should be directed at identifying the earliest time point these changes would have been detectable.

Further changes on Cry1Ab-associated LG18 (NW_018395783.1) and the Cry1A.105+Cry2Ab2-associated region of LG9 (NW_018395392.1) occurred between 2012 and 2017 in wild LA H. zea (). Notably, the same regions of Cry1A.105+Cry2Ab2-associated LG9 were under selection in both time periods (2002 to 2012 and 2012 to 2017), but patterns of selection were more complex for Cry1Ab-associated LG18. Changes on LG18 were observed on multiple superscaffolds over time: NW_018395428.1 showed changes between 2002 and 2017 and NW_018395783.1 between 2012 and 2017, while NW_018395637.1 and NW_018395512.1 differed between our 2017 Maryland populations, which indicated selection is ongoing on this LG (Fig. 2 and ). Such complexity is in line with selection for a polygenic trait from standing genetic variation, where multiple variants of small effect size are differentially selected among populations. The detectable Cry1Ab-associated QTL on LG18 covered hundreds of kilobases and may contain multiple variants of small effect, considering the polygenic nature of Cry1Ab resistance. Each of these variants may be seized upon by selection differently in independent populations over time, or by the distinct Cry1A toxins (Cry1Ab and Cry1A.105 toxins) examined here.

Although we detected genomic windows within Cry1Ab- and Cry1A.105+Cry2Ab2-associated QTL that underwent strong, significant changes over the 15-y sampling period, the strongest changes occurred in a broad region of LG13. Independence of this region from our Cry-associated QTL refutes its strong contribution to Cry1Ab, Cry1A.105, or Cry2Ab2 resistance in H. zea. Interestingly, three genes in this sweep region, cpq1, cpq2, and cyp333b3, have previously been linked to xenobiotic resistance (76–78), as well as binding of Cry1Ac toxins in H. armigera (79). Perhaps this region plays a role in resistance to Cry toxins that we did not study (e.g., Cry1Ac, Cry1Fa). Despite this finding, we compared the strong changes on LG13 to those on Cry-associated scaffolds, which offered insight into both the polygenic nature of field-evolved Cry resistance and challenges for genomic monitoring.

The strong directional changes observed for both SNP and haplotype frequencies (Fig. 3), combined with the buildup of LD over time (), indicated a hard selective sweep occurred on LG13 (59). Unlike patterns observed on LG13, the 2017 dominant haplotypes at two (of three) LG9 regions started at low frequency but rose to and remained at intermediate frequency in 2017, suggesting a partial selective sweep, or recombination linking resistance alleles to previously susceptible haplotypes over time. In contrast to both LGs 13 and 9, the dominant 2017 haplotype in a Cry1Ab-associated region of LG18 started at high frequency and declined over time, likely due to a decline of susceptible haplotypes and replacement by resistant ones. The decline was modest, and high haplotype variation was maintained in this region of LG18 in 2017 (Dataset S6). We speculate that putatively susceptible haplotypes were replaced by resistance alleles linked to multiple haplotypes, consistent with a sweep from standing genetic variation (80). Altogether, our results underscore that use of genomic approaches to study rapid adaptation for polygenic pesticide resistance is challenging but feasible under the right experimental conditions. Allele frequency changes were detected in Cry-associated regions of LGs 9 and 18 in our study, but they did not show the same strength and breadth as the hard sweep on LG13.

Use of genomic approaches to study rapid evolution in natural populations has recently gained momentum (40, 65, 81–84), yet few studies have applied them to study pest adaptation to agricultural systems (85). As these approaches are applied to a variety of agricultural systems, species- and method-specific issues will need to be addressed to ensure accuracy. Not all genomic changes identified through resistance monitoring will be the result of selection from the control tool being studied. Deployment of pest management innovations do not take place in the absence of other environmental selection pressures. When Bt crops were introduced into North America, agricultural use of synthetic insecticides was reduced, though certainly not eliminated. In a previous study, we found that another pest of cotton, Chloridea virescens, did not experience changes in the frequency of Bt resistance candidate genes. Rather, changes in the frequency of a pyrethroid resistance mutation were detected (40). Indeed, we found that the greatest allele frequency changes in H. zea were not associated with Cry1Ab, Cry1A.105, or Cry2Ab2 resistance. Therefore, confirming the strength of association between any major genomic change with a resistance phenotype is critical.

Confirmation testing can be conducted using the types of simple crosses and QTL analyses carried out in our study. Such assays can determine the importance of rapidly changing loci to resistance phenotypes but require large sample sizes to detect loci of minor effect. In cases similar to H. zea on Cry-expressing crops, where the control tactic does not kill all pest individuals or has sublethal effects, collections from paired sentinel plots could also point to genomic regions associated with resistance. Such paired sentinel plots are already being used for determining the current extent of resistance but typically monitor damage to the crop rather than the genetics of the pest (27, 86).

A second concern is that population demographic processes, including migration patterns and population expansion, can result in the appearance of genome-wide shifts in allele frequency not due to selection. When this occurs, there is a risk that genomic monitoring would falsely implicate a genomic change as a sign of evolving resistance. We were fortunate to be studying a highly migratory species (22, 23, 87) with little population structure over large geographic areas (23, 88), and this was supported by our low genome-wide pairwise FST for H. zea from Louisiana and Maryland. Additionally, H. zea feeds on many host plant species, resulting in relatively stable population densities even after release of Bt crops (89). For a species that interbreeds over such a large area, population demographic processes are unlikely to result in detection of false positives. Had we been studying an insect or weed with strong local population structure, or one in which the control measure caused severe reduction and resurgence of the pest, our experimental and computational approaches would have required additional careful consideration. Fortunately, new computational approaches which can better differentiate population demographic processes from selection, including soft and partial sweeps, are now available, even for temporal datasets (90, 91).

The benefits of genomic monitoring of pest species go beyond identifying the mechanisms of resistance and risk mitigation. Adoption of genomic approaches for resistance monitoring will be crucial for addressing the decades-long debate on exactly how pests rapidly adapt to management practices (92). Here we present evidence on the complex nature of genomic change during expansion in use of a major agricultural innovation. Extension of carefully designed genomic monitoring programs to additional pest species sampled over broad spatial and temporal scales could provide critical insight into the nature of adaptation to human-imposed selection.

Methods

Insect Collections for ddRADseq and WGS Analyses.

H. zea adults were collected by pheromone-baited trap twice weekly and pooled by collection date from May through September in Bossier Parish, LA () in 2002, 2007, 2012, 2016, and 2017. Late instar (fifth to sixth) H. zea were also collected in 2017 at the University of Maryland CMREC farms in Prince George’s County, MD, and reared to adulthood prior to genomic analysis.

Genetic Crosses for Linkage Mapping and QTL Analysis.

In 2019, we collected Cry resistant late instar (fifth to sixth) H. zea larvae from the ears of Cry1A.105+Cry2Ab2-expressing sweet corn planted in Prince George’s County, MD, and reared them to pupation on a meridic diet in the laboratory. Upon emergence, pairs of field-collected adults were mated and their larval progeny were bioassayed on a Cry-containing diet alongside susceptible Benzon larvae (Benzon Research Inc.) according to Dively et al. (27) (). After confirming differences in Cry susceptibility for Benzon and field strains, we single-pair-mated them to produce F2 progeny from two mapping families.

The H. zea reference genome was superscaffolded by alignment to the H. armigera genome using RagTag [v. 1.1.1 (93)]. A linkage map was produced using the Cry1A.105+Cry2Ab mapping family with R package qtl [v. 1.47-9 (37)]. Markers with a logarithm of the odds score of 16 or greater and a recombination fraction less than 0.2 were placed into LGs and reordered iteratively until order quality plateaued. H. zea superscaffolds were grouped into chromosomes and ordered with the linkage map. Chromosome homology to B. mori was determined by aligning the superscaffolded H. zea assembly to the B. mori assembly (v1.0, GCF_014905235.1) using RagTag. Throughout our study, all sequencing reads were aligned to the original H. zea assembly for variant calling, but SNPs were ordered using our superscaffolded genome for visualization and downstream analyses.

ddRADseq Library Preparation of Wild H. zea and Bioinformatic Analysis.

Specimens from 2002, 2007, 2012, and 2016 were prepared into ddRADseq libraries for population genomic analysis according to Fritz et al. (39, 40) and sequenced on eight runs of an Illumina MiSeq. Merged, filter-trimmed reads were aligned to the H. zea genome [v. 1.0, NCBI Bioproject PRJNA378438 (36)] with Bowtie2 (94). SNP genotypes were called with BCFtools (95), and SNP filtering criteria are provided in .

We calculated average nucleotide diversity (π) and heterozygosity (the inbreeding coefficient F) values for populations collected in each year and examined genome-wide divergence for each pair of years using pairwise FST (96). A k-means analysis in R (97) modeled the number of putative populations or clusters (k = 1:10) into which the LA H. zea fit according to genotype [adegenet v. 2.1.1 (98, 99)]. Using the three best-fitting numbers of clusters (k = 2:4) as priors, we applied a DAPC in adegenet to identify groups of genetically similar H. zea individuals. Weir and Cockerham’s FST was calculated for each SNP to quantify genetic divergence across years and identify those with higher than expected divergence over time [OutFLANK v. 0.2 (41)].

WGS Library Preparation of Wild H. zea and Bioinformatic Analysis.

DNA from additional males collected in 2002 (n = 13), 2012 (n = 11), and 2017 (n = 11) was prepared into two Illumina TruSeq LT libraries (Illumina, Inc.) and sequenced on two Illumina NextSeq500 runs. Following filter-trimming and genome alignment, SNP genotypes were called with BCFtools. Temporal change at Bt resistance candidate genes was measured by pairwise FST averaged across SNPs found within each gene. We applied a sliding-window analysis, where windows with an average FST value greater than 6 SDs from the mean (ZFST >6) were considered to have undergone statistically significant divergence (42).

Sliding window averaged FST values were also used to identify novel signatures of selection across the genome using 10-kb windows with a 1-kb step. Additional sliding-window and step sizes were also examined (). Genotype simulations under neutral expectations using realistic population demographic conditions generated a distribution of FST values, which were compared against empirical FST values from our genome scan as a secondary approach to examining statistical significance (43) ().

QTL Analysis of Cry1Ab, Cry1A.105, and Cry2Ab2 Resistance.

Progeny from each of the two mapping families were assayed for 7-d weight gain on diets containing conventional (nonexpressing) or Cry-expressing leaf tissue collected from sweet corn isolines, where one family was exposed to Cry1Ab and its nonexpressing isoline (Providence). The other family was exposed to Cry1A.105+Cry2Ab2 and its nonexpressing isoline (Obsession). Following phenotyping, both families were prepared into ddRADseq libraries () and sequenced on one lane of an Illumina NovaSeq. 6000. Genotypes were called from filtered and genome-aligned reads and then associated with 7-d growth phenotypes for all treatments using an LMM and a BSLMM in GEMMA 0.98.4 (55, 56). The intersection of our full list of 10-kb genomic windows showing significant temporal changes and our Cry-associated QTL gave the putative genomic regions under selection by Cry1Ab and Cry1A.105+Cry2Ab.

Genome Similarity between H. zea Collected from Maryland and Louisiana in 2017.

A global FST analysis compared H. zea collected at CMREC in Prince George’s County, MD, in 2017 to those collected in Bossier Parish, LA, in 2017. Nonoverlap of the bootstrapped 95% confidence intervals with zero determined statistical significance. Maryland samples collected from Cry1A.105+Cry2Ab2-expressing and nonexpressing sweet corn were examined by a pairwise FST sliding window analysis, as described above.

Frequency of Nonsynonymous Substitutions Prior to 2002.

Genomic DNA was isolated from H. zea collected in 1998 from Hartstack traps in Bossier Parish, LA, PCR-amplified using primers flanking nonsynonymous mutations in cpq1 and cyp333b3 (), and Sanger-sequenced to determine their frequency prior to 2002.

Analysis of LD.

Decay of LD was measured using R2 for H. zea scaffolds under Cry selection in LGs 9 and 18 and compared with genome-wide patterns, as well as those at LG13. Changes in haplotype frequency were measured on Cry-associated LGs 9 and 18, as well as regions of LG13. We included SNPs that were present in >90% of our samples after thinning to 1 per 500 bp. Haplotype frequencies were calculated from called genotypes using the R package haplo.stats (100).