Sequence variation and regulatory variation in acetylcholinesterase genes contribute to insecticide resistance in different populations of Leptinotarsa decemlineata

Abstract Although insect herbivores are known to evolve resistance to insecticides through multiple genetic mechanisms, resistance in individual species has been assumed to follow the same mechanism. While both mutations in the target site insensitivity and increased amplification are known to contribute to insecticide resistance, little is known about the degree to which geographic populations of the same species differ at the target site in a response to insecticides. We tested structural (e.g., mutation profiles) and regulatory (e.g., the gene expression of Ldace1 and Ldace2, AChE activity) differences between two populations (Vermont, USA and Belchow, Poland) of the Colorado potato beetle, Leptinotarsa decemlineata in their resistance to two commonly used groups of insecticides, organophosphates, and carbamates. We established that Vermont beetles were more resistant to azinphos‐methyl and carbaryl insecticides than Belchow beetles, despite a similar frequency of resistance‐associated alleles (i.e., S291G) in the Ldace2 gene. However, the Vermont population had two additional amino acid replacements (G192S and F402Y) in the Ldace1 gene, which were absent in the Belchow population. Moreover, the Vermont population showed higher expression of Ldace1 and was less sensitive to AChE inhibition by azinphos‐methyl oxon than the Belchow population. Therefore, the two populations have evolved different genetic mechanisms to adapt to organophosphate and carbamate insecticides.

insecticides (IRAC, 2020). Although insecticides can have a specific selective pressure, where selection acts on a specific target, there is still considerable variation in the magnitude of insecticide resistance among the populations of the same insect species (Dively et al., 2020;Ryan et al., 2019;Wu et al., 1999). These differences have often been linked to the genetic differences at the target site (Ilias et al., 2014;Weill et al., 2003). Despite increasing knowledge of the biochemical basis of insecticide resistance, relatively little is known about the geographic difference in insecticide target site gene sequences or their regulation (Hawkins et al., 2019;Ryan et al., 2019).
Pest species can become resistant to insecticides via various structural and regulatory mechanisms (Feyereisen et al., 2015).
Structural mechanisms are better characterized and involve target site mutations in an enzyme that make the insecticide ineffective.
These can be produced with genetic differences such as nonsynonymous nucleotide variation at target site genes (Li & Han, 2004;Malekmohammadi & Galehdari, 2016;Weill et al., 2004;Zhu et al., 1996). Regulatory mechanisms have been less studied and can be related, for example, to biochemical processes such as the overexpression of the target site or enhanced metabolism or excretion of the insecticide (Barres et al., 2016;Ffrench-Constant, 2013). In general, we know less about the variation among populations at both the nucleotide and regulatory levels of the target site genes. By understanding population-level variation also at the regulatory level, we may be able to explain better why the same insecticide might cause different outcomes among populations. This is because our predictions based on gene sequence variation alone might bias our estimates of the overall pesticide resistance (see discussion in Hawkins et al., 2019).
Organophosphates (OPs) and carbamates are both widely known as acetylcholinesterase (AChE; EC 3.1.1.7) inhibitors, which have been heavily used as insecticides since the 1970s. The AChE enzyme functions at the synapses of cholinergic neurons in the central and peripheral nervous system (Fukuto, 1990;Taylor, 2011;Taylor et al., 2009). The enzyme terminates neurotransmission at cholinergic synapses in the synaptic cleft by hydrolyzing the neurotransmitter acetylcholine (Taylor et al., 2009). OP and carbamate insecticides inhibit AChE, thus interfering with neurotransmission, leading to paralysis and death (Colovic et al., 2013). In many invertebrates, there are two distinct acetylcholinesterases (AChE1 and AChE2), which are encoded by two ace genes (i.e., ace1 and ace2; Kim & Lee, 2013).
The two AChEs are probably homologous, derived from an ancient duplication event that occurred long before the differentiation of insects (Weill et al., 2002). In 67 insect species out of 100, AChE1 accounts for most of the AChE activity and thus is considered as the main catalytic enzyme (Kim & Lee, 2013) responsible for neuronal functions (Weill et al., 2002). However, the role of AChE2 differs between species. It can act as the main catalytic activity, be equally active to AChE1, or show little catalytic activity (Kim & Lee, 2013).
For example, in Bombyx mori, Bm-ace2 is more highly expressed and thus likely the main catalytic enzyme, while Bm-ace1 may contribute to metamorphosis (Chen et al., 2009). However, more functional studies are needed to demonstrate the exact functions of these two genes (Jiang et al., 2018).
The Colorado potato beetle (CPB), Leptinotarsa decemlineata Say (Coleoptera: Chrysomelidae), is a model organism for studying the evolution of insecticide resistance . CPB has played an important role in the modern pesticide industry. Since it was first targeted by insecticides in 1864, it has been heavily managed with insecticides (Alyokhin et al., 2008;Gauthier et al., 1981).
CPB has currently developed resistance to more than 50 different active ingredients used in insecticides, including OPs (16 active ingredients) and carbamates (five active ingredients; Alyokhin et al., 2008;Brevik et al., 2018;Mota-Sanchez & Wise, 2020 (Zhu & Clark, 1995), has been associated with the resistance to OP and carbamate insecticides in CPB (Zhu et al., 1996). Three mutations (S291G, R30K, and Y45H) in Ldace2 gene have been described and linked to OP and/or carbamate resistance (Kim et al., 2007;Zhu & Clark, 1997;Zhu et al., 1996), with the S291G mutation being the main target site conferring resistance (Zhu et al., 1996). Revuelta et al. (2011) described Ldace1 for the CPB which is an orthologous to the Anopheles gambiae ace1 gene. The authors suggested that the 2-to 11-fold higher expression of Ldace1 compared with Ldace2 indicated that Ldace1 rather than Ldace2 was the main contributor to the AChE activity and may have been the primary target of the OP insecticides (Revuelta et al., 2011).
We set to study whether the differences among populations in mortality are linked to sequence variation and/or regulatory resistance mechanisms. We hypothesized that sequence and regulatory variation in Ldace1 and Ldace2 genes (but see Piiroinen et al., 2013) might contribute to differences in OP and carbamate resistance development in different geographic populations. We tested how the Ldace1 and Ldace2 genes contribute to resistance to OP and carbamate insecticides by comparing CPB populations that differ in their resistance to these insecticides. Using bioassays, we started by estimating the resistance level of six different populations to the organophosphate azinphos-methyl (AZ) and the carbamate carbaryl (CAR) insecticides. Based on the resistance status, we selected two populations, Vermont (from the USA) as the most resistant, and Belchow (from Poland) as the less resistant population, for further investigation of Ldace1 and Ldace2 genes. To test for populationlevel differences, we first measured the efficacy of AChE inhibition by azinphos-methyl oxon (AZoxon) and CAR insecticides. Then, we measured Ldace1 and Ldace2 gene expression, both before exposure and after exposure to both insecticides. Finally, we sequenced the Ldace1 and Ldace2 genes and investigated the frequency of amino acid replacements (mutations) in the two populations. We predicted that the more resistant population (i) would be less sensitive to AChE insecticide inhibition, (ii) demonstrate higher target site expression levels, and (iii) show higher frequency of mutations at target sites.

| Study species and rearing conditions
Colorado potato beetles used in this study were descendants of Adult beetles of the summer generation overwintered individually at 5°C in environmental chambers.

| Insecticide bioassays
We used bioassays to determine whether the different geographic populations differ in the degree of insecticide resistance to azinphos-methyl (AZ; organophosphate) and the carbaryl (CAR; carbamate) insecticides (Ovčarenko et al., 2014). For the bioassays, we obtained both AZ and CAR as Pestanal analytical standards from Sigma-Aldrich, which were dissolved in acetone.
Bioassays were performed by applying insecticide treatment to early third instar (5-7 days old, Bointeau & Le Blanc, 1992) larvae using a pipette. Larvae were randomly divided into petri dishes (five individuals/family/petri dish/) and were randomly assigned to an insecticide treatment. Then, we applied a 3 µl drop of AZ, CAR, or acetone (as control) topically to the fifth and sixth dorsal abdominal segments. In order to calculate the median lethal dose (LD 50 ), we treated each of the six populations (N = 160-414 larvae) with 4-8 different insecticide doses ranging from 0.0075 to 30 µg/larvae for AZ and 0.06 to 90 µg/larvae for CAR (3-6 larvae per family, 20-74 larvae per insecticide dose, see Figure S1). After the insecticide application, we provided larvae a standardized potato leaf after 2 h and assessed survival after 22 h (24 h of survival from the application). A larva was considered dead if it was not able to move after being placed on its back.
We analyzed larval survival separately for both insecticides using generalized linear models (GLM; Binary logistic, logit link function) in SPSS (IBM SPSS Statistics 24.0.0). Population and insecticide dose were set as explanatory variables and survival as the dependent variable. We determined the resistance levels (LD 50 , i.e., lethal dose) and generated the 95% confidence limits for each population and insecticide using Probit analysis. We considered differences in LD 50 values between populations and insecticides significant if their 95% confidence limits did not overlap (Ovčarenko et al., 2014). No mortality was observed in any of the populations within the acetone control.

| Sampling AChE enzyme activity and gene expression studies
Based on bioassay results, the Vermont and Belchow populations were selected for further studies to investigate the AChE activity, gene expression, and nonsynonymous point mutations in the two genes (Tables 1-2; Figure 1). We selected the Vermont population because it showed the highest tolerance to both insecticides, while the Belchow population was less resistant. Since the LD 10 dose of Vermont population would have killed all the beetles from Belchow population, we applied population-specific LD 10 doses to cause 10% Note: LD 50 (i.e., lethal dose that kills 50% of exposed individuals within 24 h since exposure) values for the six Colorado potato beetle populations exposed to the AZ insecticide. a Resistance ratio showing the increase in resistance compared to the least (Ufa) resistant population. mortality (LD 10 ) and induce a similar level of insecticidal stress. We used the same protocol described above for applying insecticide doses of AZ (Vermont 0.375 µg/larvae and Belchow 0.0075 µg/larvae), CAR (Vermont 1.5 µg/larvae and Belchow 0.15 µg/larvae), and acetone as control. After the insecticide application, we collected the surviving larvae 2, 4, and 24 h from the insecticide application, froze them in liquid nitrogen, and stored at −80°C until further analysis. In total, we collected 713 larvae from 15 Vermont families and 417 larvae from 15 Belchow families.

| Measurement of AChE inhibition efficiency
We tested in vitro whether populations that can tolerate higher insecticide doses (see Tables 1-2

| Expression of AChE genes
We tested for differences in the gene expression level of Ldace1 and Ldace2 between Vermont and Belchow populations. We We designed primers used for qPCR (Table S1) amplification of Ldace1 (JF343436.1; (Revuelta et al., 2011)) and Ldace2 (L41180.1; (Zhu & Clark, 1995)) using Primer 3 (http://frodo.wi.mit.edu/prime r3/, v. 0.4.0) according to . The primers were designed based on sequences from the annotated transcriptome of CPB (Kumar et al., 2014) and sequences available in GenBank. The qPCRs were run on a Bio-Rad CFX96™ instrument with an initial denaturation step of 95°C for 3 min, followed by 39 cycles of 10 s at 95°C, 10 s at 56°C, and 30 s at 72°C. The qPCR was followed by melting curve analysis (65-95°C) to check the purity of qPCR. We analyzed ten biological replicates (individual beetles) with three technical replicates for each treatment group (control, AZ, and CAR) for each population. In order to calibrate the expression across samples, we used two positive control samples with two replicates that were added to each plate.
The efficiency of qPCR amplification was calculated for each gene using twofold serial dilutions of pooled cDNA. The amplification efficiencies were between 96% and 113%.
We calculated expression values (mean Cq) for all the samples using the normalized expression (ΔΔCq) method with default threshold values by using the CFX Manager 3.0 software (Bio-Rad Laboratories Inc.). We used forkhead transcription factor (FOXO) and ribosomal protein L13e (L13e) as reference genes (Table S1; Kumar et al., 2014;Yocum et al., 2009).
We analyzed the relative expression data with the REST program (http://rest.gene-quant ifica tion.info/; Pfaffl et al., 2002), according to . The REST program (with 10,000 iterations) was used for pairwise comparisons within and between population and treatment groups.

| Identification of sequence variation in the Ldace1 and Ldace2 genes
To identify and compare mutations in Ldace1 and Ldace2 between populations, we sequenced the genes from the experimental larvae. We extracted total RNA from 38 (Belchow n = 19, Vermont n = 19) whole larvae and synthesized cDNA using the same procedures described earlier. Primers for Ldace1 and Ldace2 (Table   S2) were designed using the CPB sequences available in GenBank
In contrast, there were no differences between the two populations in AChE activity after inhibition by CAR (F 1,23 = 0.034, p = .855; Figure 2).

| Expression of AChE genes
The Ldace1

| Identification of sequence variation in the Ldace1 and Ldace2 genes
Sequencing analyses revealed five nonsynonymous mutations in the AChE genes: G192S and Y402F in the Ldace1 gene and R30K, Y54H, and S291G in the Ldace2 gene. In the Ldace1 gene, the G192S and Y402F mutations were only present in the Vermont population ( Figure 4a). Allele frequencies differed significantly (G192S: χ 2 = 7.8, df = 2, p = .020; Y402F: χ 2 = 10.7, df = 2, p = .005; Figure 4b) between populations. Homozygotes for the G192S and Y402F mutations represented 31% and 8% of the individuals from the Vermont population, respectively. The two alleles never co-occurred in the same individual.
In the Ldace2 gene, the R30K mutation was only present in the individuals from the Vermont population, whereas the Y54H mutation was only present in the Belchow population. The S291G mutation was identified in both populations (Table 2; Figure 4b). Allele frequencies of the R30K and Y54H mutations differed significantly between the populations (R30K: χ 2 = 14.5, df = 2, p = .001; Y54H: χ 2 = 13.0, df = 2, p = .002; Figure 4b). In the Vermont population, the R30K and S291G mutations were present in 39% and 67% of the individuals, respectively. In the Belchow population, the frequency of the Y54H mutation was 42%, whereas that of the S291G mutation was 63% (Figure 4b). The frequency of the S291G mutation was similar in both populations (S291G: χ 2 = 1.9, df = 2, p = .382).
Both the R30K (in Vermont) and the Y54H (in Belchow) alleles occurred together with the S291G in 39% and 42% of individuals, respectively.

| D ISCUSS I ON
We investigated the population-level sequence variation of resistanceassociated Ldace1 and Ldace2 genes together with their responses to two commonly used pesticides against the Colorado potato beetle to understand the role of these genes on insecticide resistance. We found that geographic populations differ in their resistance to two commonly used insecticides due to multiple differences in the insecticide target site. We demonstrated that the North American population (Vermont, USA) had higher resistance compared with a European population (Belchow, Poland), which was associated with higher occurrence of mutations (Figure 4), higher baseline expression of the Ldace1 gene (Figure 3a), and an AChE enzyme less sensitive to AZoxon insecticide inhibition (Figure 2). In addition, the Vermont population had two mutations in the Ldace1 gene that were absent from Belchow (Figure 4a). Therefore, it is likely that repeated application of insecticides (Dively et al., 2020) or different insecticide intensities (Crossley et al., 2018) together with invasion history (Grapputo et al., 2005) has resulted in the evolution of sequence variation and regulatory changes in the target genes that probably contributed to the overall higher resistance of the Vermont population ( Figure 1).

We identified two novel mutations (G192S and F402Y) in the
Ldace1 gene in the Vermont population. The absence of these mutations in the Belchow population could be due to the loss of genetic variation when the beetle invaded Europe (Grapputo et al., 2005).
Although the mutation in the Ldace2 S291G site has been previously described as the main mutation contributing to organophosphate insecticide resistance in the CPB (Zhu & Clark, 1997;Zhu et al., 1996)  the fact that it was equally common in both Vermont and Belchow populations suggests that it is unlikely to be the main factor explaining explains the 107-time difference in survival between the populations (Figure 1). It may be possible that the ancestral state of the resistance mutation in the Ldace2 gene has been G291S rather than S291G. This is because individuals with the resistance mutation are more susceptible to the host plant alkaloids (Wierenga & Hollingworth, 1992). This mutation might be related to the beetles' adaptation to the lower concentrations of steroidal alkaloids in agricultural solanaceous plants rather than insecticide resistance (Piiroinen et al., 2013;Wierenga & Hollingworth, 1992;Zhu & Clark, 1995). Therefore, it is possible that the novel mutations in the Ldace1 gene might play a more important role in the insecticide resistance than the resistance-associated mutation S291G in the Ldace2 gene.
The two novel mutations (see Figure 4a) described here could contribute more to the insecticide resistance in the Vermont population than the previously identified mutations in the Ldace2 gene.
Similar substitutions in the ace1 gene (glycine to serine (G192S) and phenylalanine to tyrosine (F402Y)) have been associated with insecticide resistance (e.g., insensitivity) in other insects. Different substitutions in the ace-1 gene occur frequently across multiple species, for example, the glycine to serine substitution (i.e., G119S) in Anopheles gambiae (Weetman et al., 2015) and Culex pipiens (Weill et al., 2003(Weill et al., , 2004, phenylalanine to tyrosine (i.e., F327Y, F331Y, and F445Y) in Musca domestica, Bemisia tabaci, and Culex tritaeniorhnchus (Alon et al., 2008;Nabeshima et al., 2004;Oh et al., 2006;Walsh et al., 2001). Some of the mutations in the ace1 seemed to confer high levels of resistance in combination with other mutations in other species (Mutero et al., 1994;Vontas et al., 2002) but interestingly, in our sample, the two mutations (G192S and F402Y) never co-occurred in Ldace1. We would need more functional studies to confirm the role these mutations play the insecticide resistance of the Colorado potato beetle.
In addition to substitution differences, we also identified gene expression differences between Ldace1 and Ldace2 in the two populations (see also Dively et al., 2020). The Ldace1 gene was 49-to 145-fold more expressed than Ldace2 under control conditions, suggesting that ace1 encoding ACHE1 is the major catalytic enzyme also in the CPB. These results are consistent with previous data, which indicated that in 66 insect species out of 100, AChE1 is more highly expressed than AChE2 and therefore the major catalytic enzyme of acetylcholine (Kim & Lee, 2013). Compared to a previous study (Revuelta et al., 2011), the difference between the Ldace1 and Ldace2 gene expression within Belchow (49-fold difference) and Vermont (145-fold difference) populations was higher than the reported difference between the developmental stages (i.e., from embryos to adults, 2-to 11-fold difference). The increased expression in the Vermont beetles suggests that the Ldace1 gene could also play a role in insecticide resistance. Indeed, regulatory changes in the target gene (i.e., increased expression of acetylcholinesterase gene) have been previously shown to increase OP resistance in the greenbug (Shizaphis graminum; Gao & Zhu, 2002). The lack of differences in the Ldace2 expression between populations in the control groups further suggests a greater role for Ldace1 than Ldace2 in conferring resistance to insecticides in the CPB.
The Vermont population (see Figure 1) was less sensitive to AZoxon inhibition than the Belchow population. These differences could be explained either by regulatory resistance, that is, the lower gene expression of the Belchow populations compared with the Vermont population, or alternatively by structural resistance, that is, the mutation profile differences between the populations.
Unfortunately, our sampling does not allow us to separate these two hypotheses. Previously, AChE sensitivity in the CPB has been associated with S291G and R30K mutations in the Ldace2 gene (Kim et al., 2006), but they were not aware of the presence of the Ldace1 gene. Therefore, it would be interesting to test the sensitivity differences related to the two mutations in the Ldace1 (G192S and F402Y) instead. Although the CAR insecticide inhibited AChE, there were no differences between populations. This suggests that alternative resistance mechanisms are present since there was still 20-time survival difference between the two populations. Alternatively, the differences between populations could be due to differences in OP and carbamate insecticide-induced inhibitory actions (Colovic et al., 2013). Finally, we cannot entirely exclude an insufficient dose in the enzyme inhibition assay.
Our results demonstrate that CPB resistance to commonly used OP and carbamate insecticides is due to multiple mechanisms acting at one target site that are not mutually exclusive. Besides changes at these target sites, it is likely that other resistance mechanisms have also been under selection (see Barres et al., 2016;Mutero et al., 1994), so we cannot expect that geographically separated populations have similar mechanisms of resistance. Therefore, the same management strategy in vast areas may not be equally successful across different populations and may potentially select for further differences among populations. Differences in the resistance mechanism can also challenge the development of new application strategies if the resistance of a population differs more in structural (has higher genetic variation) or regulatory (is better able to deal with xenobiotics) mechanisms. Further research is clearly needed to uncover which mechanisms confer the resistance in different populations. Understanding the mechanisms behind the rapid evolution of insecticide resistance will help us in making better insecticide risk assessment and management strategies, in this and other pest species.

| CON CLUS IONS
Our results indicate that the differences we observe in insecticide resistance among populations are a result of multiple factors acting at the same time. We demonstrated that within the target site, we need to incorporate both sequence variation and regulatory changes in AChE biochemistry and physiology to understand resistance evolution. The fact that populations differ in multiple levels even within one target site means that populations can respond to the same selection pressure in different ways. From a management perspective, it is important to understand that insecticides select not only for resistance genes but also for their function at the same time. Beetles from the Vermont population had (1) an AChE that is less sensitive to insecticide inhibition, (2) higher Ldace1 gene expression (which suggests that AChE1 is the major catalytic enzyme in the CPB), and (3) more target site mutations (G192S and F402Y) in the Ldace1 gene than the Belchow population. At the same time, our results underlined that studying changes at the target sites were not sufficient in explaining the observed resistance differences between populations and that there are other pathways to achieve the resistance besides changes at the target sites (see Barres et al., 2016). Therefore, to develop more efficient pest management strategies, we need more studies to lighter the geographical variation in resistance to insecticides.

ACK N OWLED G EM ENTS
The authors would like to thank Andreas Plischke and Maxim Udalov for the beetles collected from Holland and Ufa. We thank Sari Viinikainen and Emily Knott for their advice during the data sequencing and Sara Calhim for fruitful discussions on survival data analysis. We thank Kati Kivisaari and Joel Rahkonen for their help of rearing the beetles. This work was financed by the Academy of Finland (Grant to LL: 308302, 250248) and Centre of Excellence in Biological interactions research (252411). Since the CPB is a quarantine species in Finland, this research was carried out under a permission from Evira (3861/541/2007).

CO N FLI C T O F I NTE R E S T
The authors declare no competing interest.