# A model species for agricultural pest genomics: the genome of the Colorado potato beetle, Leptinotarsa decemlineata (Coleoptera: Chrysomelidae)

Sean Ma, Yolanda Chen, Martin Ma, Joshua Benoit, Anita Bhandari, Julia Bowsher, Kristian Brevik, Kaat Cappelle, Mei-Ju Li, Anna Childers, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Sean Ma, Yolanda Chen, Martin Ma, Joshua Benoit, Anita Bhandari, et al.. A model species for agricultural pest genomics: the genome of the Colorado potato beetle, Leptinotarsa decemlineata (Coleoptera: Chrysomelidae). Scientific Reports, 2018, 8 (1), pp.1931. 10.1038/s41598-018-20154-1. hal-02094614

HAL Id: hal-02094614

https://amu.hal.science/hal-02094614

Submitted on 26 May 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
SCIENTIFIC REPORTS

<figure><img src="image_3.png" /></figure>

OPEN

# A model species for agricultural pest genomics: the genome of the Colorado potato beetle, *Leptinotarsa decemlineata* (Coleoptera: Chrysomelidae)

Received: 17 October 2017

Accepted: 13 January 2018

Published online: 31 January 2018

Sean D. Schoville<sup>1</sup>, Yolanda H. Chen<sup>2</sup>, Martin N. Andersson<sup>3</sup>, Joshua B. Benoit<sup>4</sup>, Anita Bhandari<sup>5</sup>, Julia H. Bowsher<sup>6</sup>, Kristian Brevik<sup>2</sup>, Kaat Cappelle<sup>7</sup>, Mei-Ju M. Chen<sup>8</sup>, Anna K. Childers<sup>9,10</sup>, Christopher Childers<sup>10</sup>, Olivier Christiaens<sup>7</sup>, Justin Clements<sup>1</sup>, Elise M. Didion<sup>4</sup>, Elena N. Elpidina<sup>11</sup>, Patamarerk Engsontia<sup>12</sup>, Markus Friedrich<sup>13</sup>, Inmaculada García-Robles<sup>14</sup>, Richard A. Gibbs<sup>15</sup>, Chandan Goswami<sup>16</sup>, Alessandro Grapputo<sup>17</sup>, Kristina Gruden<sup>18</sup>, Marcin Grynberg<sup>19</sup>, Bernard Henrissat<sup>20,21,22</sup>, Emily C. Jennings<sup>4</sup>, Jeffery W. Jones<sup>13</sup>, Megha Kalsi<sup>23</sup>, Sher A. Khan<sup>24</sup>, Abhishek Kumar<sup>25,26</sup>, Fei Li<sup>27</sup>, Vincent Lombard<sup>20,21</sup>, Xingzhou Ma<sup>27</sup>, Alexander Martynov<sup>10,28</sup>, Nicholas J. Miller<sup>29</sup>, Robert F. Mitchell<sup>30</sup>, Monica Munoz-Torres<sup>31</sup>, Anna Muszewska<sup>19</sup>, Brenda Oppert<sup>32</sup>, Subba Reddy Palli<sup>23</sup>, Kristen A. Panfilio<sup>33,34</sup>, Yannick Pauchet<sup>35</sup>, Lindsey C. Perkin<sup>32</sup>, Marko Petek<sup>18</sup>, Monica F. Poelchau<sup>8</sup>, Éric Record<sup>36</sup>, Joseph P. Rinehart<sup>10</sup>, Hugh M. Robertson<sup>37</sup>, Andrew J. Rosendale<sup>4</sup>, Victor M. Ruiz-Arroyo<sup>14</sup>, Guy Smagghe<sup>10,7</sup>, Zsófia Szendrei<sup>38</sup>, Gregg W.C. Thomas<sup>39</sup>, Alex S. Torson<sup>6</sup>, Iris M. Vargas Jentzsch<sup>33</sup>, Matthew T. Weirauch<sup>10,40,41</sup>, Ashley D.™Yates<sup>42,43</sup>, George D. Yocum<sup>10</sup>, June-Sun Yoon<sup>23</sup> & Stephen Richards<sup>15</sup>

<footnote n="1"><sup>1</sup>Department of Entomology, University of Wisconsin-Madison, Madison, USA. <sup>2</sup>Department of Plant and Soil Sciences, University of Vermont, Burlington, USA. <sup>3</sup>Department of Biology, Lund University, Lund, Sweden. <sup>4</sup>Department of Biological Sciences, University of Cincinnati, Cincinnati, USA. <sup>5</sup>Department of Molecular Physiology, Christian-Albrechts-University at Kiel, Kiel, Germany. <sup>6</sup>Department of Biological Sciences, North Dakota State University, Fargo, USA. <sup>7</sup>Department of Crop Protection, Ghent University, Ghent, Belgium. <sup>8</sup>USDA-ARS National Agricultural Library, Beltsville, MD, USA. <sup>9</sup>USDA-ARS Bee Research Lab, Beltsville, MD, USA. <sup>10</sup>USDA-ARS Insect Genetics and Biochemistry Research Unit, Fargo, ND, USA. <sup>11</sup>A.N. Belozersky Institute of Physico-Chemical Biology, Lomonosov Moscow State University, Moskva, Russia. <sup>12</sup>Department of Biology, Faculty of Science, Prince of Songkla University, Amphoe Hat Yai, Thailand. <sup>13</sup>Department of Biological Sciences, Wayne State University, Detroit, USA. <sup>14</sup>Department of Genetics, University of Valencia, Valencia, Spain. <sup>15</sup>Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, Texas, USA. <sup>16</sup>National Institute of Science Education and Research, Bhubaneswar, India. <sup>17</sup>Department of Biology, University of Padova, Padova, Italy. <sup>18</sup>Department of Biotechnology and Systems Biology, National Institute of Biology, Ljubljana, Slovenia. <sup>19</sup>Institute of Biochemistry and Biophysics, Polish Academy of Sciences, Warsaw, Poland. <sup>20</sup>Architecture et Fonction des Macromolécules Biologiques, CNRS, Aix-Marseille Université, 13288, Marseille, France. <sup>21</sup>INRA, USC 1408 AFMB, F-13288, Marseille, France. <sup>22</sup>Department of Biological Sciences, King Abdulaziz University, King Abdulaziz, Saudi Arabia. <sup>23</sup>Department of Entomology, University of Kentucky, Lexington, USA. <sup>24</sup>Department of Entomology, Texas A&M University, College Station, USA. <sup>25</sup>Department of Genetics & Molecular Biology in Botany, Christian-Albrechts-University at Kiel, Kiel, Germany. <sup>26</sup>Division of Molecular Genetic Epidemiology, German Cancer Research Center, Heidelberg, Germany. <sup>27</sup>Department of Entomology, Nanjing Agricultural University, Nanjing, China. <sup>28</sup>Center for Data-Intensive Biomedicine and Biotechnology, Skolkovo Institute of Science and Technology, Moscow, Russia. <sup>29</sup>Department of Biology, Illinois Institute of Technology, Chicago, USA. <sup>30</sup>Department of Biology, University of Wisconsin-Oshkosh, Oshkosh, USA. <sup>31</sup>Environmental Genomics and Systems Biology Division, Lawrence Berkeley National Laboratory, Berkeley, USA. <sup>32</sup>USDA-ARS Center for Grain and Animal Health Research, New York, USA. <sup>33</sup>Institute for Developmental Biology, University of Cologne, Köln, Germany. <sup>34</sup>School of Life Sciences, University</footnote>
of Warwick, Gibbet Hill Campus, England, UK. ³⁵Department of Entomology, Max Planck Institute for Chemical Ecology, Jena, Germany. ³⁶INRA, Aix-Marseille Université, UMR1163, Biodiversité et Biotechnologie Fongiques, Marseille, France. ³⁷Department of Entomology, University of Illinois at Urbana-Champaign, Champaign, IL, USA. ³⁸Department of Entomology, Michigan State University, East Lansing, USA. ³⁹Department of Biology and School of Informatics and Computing, Indiana University, Bloomington, USA. ⁴⁰Center for Autoimmune Genomics and Etiology, Division of Biomedical Informatics and Division of Developmental Biology, Cincinnati Children's Hospital Medical Center, Cincinnati, USA. ⁴²Department of Pediatrics, University of Cincinnati College of Medicine, Cincinnati, USA. ⁴²Department of Entomology, The Ohio State University, Columbus, USA. ⁴³Center for Applied Plant Sciences, The Ohio State University, Columbus, USA. Correspondence and requests for materials should be addressed to S.D.S. (email: sean.schoville@wisc.edu)

The Colorado potato beetle is one of the most challenging agricultural pests to manage. It has shown a spectacular ability to adapt to a variety of solanaceae plants and variable climates during its global invasion, and, notably, to rapidly evolve insecticide resistance. To examine evidence of rapid evolutionary change, and to understand the genetic basis of herbivory and insecticide resistance, we tested for structural and functional genomic changes relative to other arthropod species using genome sequencing, transcriptomics, and community annotation. Two factors that might facilitate rapid evolutionary change include transposable elements, which comprise at least 17% of the genome and are rapidly evolving compared to other Coleoptera, and high levels of nucleotide diversity in rapidly growing pest populations. Adaptations to plant feeding are evident in gene expansions and differential expression of digestive enzymes in gut tissues, as well as expansions of gustatory receptors for bitter tasting. Surprisingly, the suite of genes involved in insecticide resistance is similar to other beetles. Finally, duplications in the RNAi pathway might explain why *Leptinotarsa decemlineata* has high sensitivity to dsRNA. The *L. decemlineata* genome provides opportunities to investigate a broad range of phenotypes and to develop sustainable methods to control this widely successful pest.

The Colorado potato beetle, *Leptinotarsa decemlineata* Say 1824 (Coleoptera: Chrysomelidae), is widely considered one of the world's most successful globally-invasive insect herbivores, with costs of ongoing management reaching tens of millions of dollars annually¹ and projected costs if unmanaged reaching billions of dollars². This beetle was first identified as a pest in 1859 in the Midwestern United States, after it expanded from its native host plant, *Solanum rostratum* (Solanaceae), onto potato (*S. tuberosum*)³. As testimony to the difficulty in controlling *L. decemlineata*, the species has the dubious honor of starting the pesticide industry, when Paris Green (copper (II) acetoarsenite) was first applied to control it in the United States in 1864⁴. *Leptinotarsa decemlineata* is now widely recognized for its ability to rapidly evolve resistance to insecticides, as well as a wide range of abiotic and biotic stresses⁵, and for its global expansion across 16 million km² to cover the entire Northern Hemisphere within the 20th century⁶. Over the course of 150 years of research, *L. decemlineata* has been the subject in more than 9,700 publications (according to the Web of Science™ Core Collection of databases) ranging from molecular to organismal biology from the fields of agriculture, entomology, molecular biology, ecology, and evolution.

In order to be successful, *L. decemlineata* evolved to exploit novel host plants, to inhabit colder climates at higher latitudes⁷⁻⁹, and to cope with a wide range of novel environmental conditions in agricultural landscapes¹⁰,¹¹. Genetic data suggest the potato-feeding pest lineage directly descended from populations that feed on *S. rostratum* in the U.S. Great Plains¹². This beetle subsequently expanded its range northwards, shifting its life history strategies to exploit even colder climates⁷,⁸,¹³, and steadily colonized potato crops despite substantial geographical barriers¹⁴. *Leptinotarsa decemlineata* is an excellent model for understanding pest evolution in agroecosystems because, despite its global spread, individuals disperse over short distances and populations often exhibit strong genetic differentiation¹⁵⁻¹⁷, providing an opportunity to track the spread of populations and the emergence of novel phenotypes. The development of genomic resources in *L. decemlineata* will provide an unparalleled opportunity to investigate the molecular basis of traits such as climate adaptation, herbivory, host expansion, and chemical detoxification. Perhaps most significantly, understanding its ability to evolve rapidly would be a major step towards developing sustainable methods to control this widely successful pest in agricultural settings.

Given that climate is thought to be the major factor in structuring the range limits of species¹⁸, the latitudinal expansion of *L. decemlineata*, spanning more than 40° latitude from Mexico to northern potato-producing countries such as Canada and Russia⁶, warrants further investigation. Harsh winter climates are thought to present a major barrier for insect range expansions, especially near the limits of a species' range⁷,¹⁹. To successfully overwinter in temperate climates, beetles need to build up body mass, develop greater amounts of lipid storage, have a low resting metabolism, and respond to photoperiodic keys by initiating diapause²⁰,²¹. Although the beetle has been in Europe for less than 100 years, local populations have demonstrating remarkably rapid evolution in life history traits linked to growth, diapause and metabolism⁸,¹³,²⁰. Understanding the genetic basis of these traits, particularly the role of specific genes associated with metabolism, fatty acid synthesis, and diapause induction, could provide important information about the mechanism of climate adaptation.

Although *Leptinotarsa decemlineata* has long-served as a model for the study of host expansion and herbivory due to its rapid ability to host switch¹⁷,²², a major outstanding question is what genes and biological pathways are associated with herbivory in this species? While >35,000 species of Chrysomelidae are well-known herbivores, most species feed on one or a few host species within the same plant family²³. Within *Leptinotarsa*, the majority of species feed on plants within Solanaceae and Asteraceae, while *L. decemlineata* feeds exclusively on solanaceous
species²⁴. It has achieved the broadest host range amongst its congeners (including, but not limited to: buffalobur (S. rostratum), potato (S. tuberosum), eggplant (S. melongena), silverleaf nightshade (S. elaeagnifolium), horsenettle (S. carolinense), bittersweet nightshade (S. dulcamara), tomato (S. lycopersicum), and tobacco (Nicotiana tabacum))¹⁷,²²,²⁵, and exhibits geographical variation in the use of locally abundant Solanum species²⁶.

Another major question is what are the genes that underlie the beetle's remarkable capacity to detoxify plant secondary compounds and are these the same biological pathways used to detoxify insecticidal compounds?²⁷? Solanaceous plants are considered highly toxic to a wide range of insect herbivore species²⁸, because they contain steroidal alkaloids and glycoalkaloids, nitrogen-containing compounds that are toxic to a wide range of organisms, including bacteria, fungi, humans, and insects²⁹, as well as glandular trichomes that contain additional toxic compounds³⁰. In response to beetle feeding, potato plants upregulate pathways associated with terpenoid, alkaloid, and phenylpropanoid biosynthesis, as well as a range of protease inhibitors³¹. A complex of digestive cysteine proteases is known to underlie *L. decemlineata*'s ability to respond to potato-induced defenses³²,³³. There is evidence that larvae excrete³⁴ and perhaps even sequester toxic plant-based compounds in the hemolymph³⁵,³⁶. Physiological mechanisms involved in detoxifying plant compounds, as well as other xenobiotics, have been proposed to underlie pesticide resistance²⁷. To date, while cornerstone of *L. decemlineata* management has been the use of insecticides, the beetle has evolved resistance to over 50 compounds and all of the major classes of insecticides. Some of these chemicals have even failed to control *L. decemlineata* within the first year of release¹⁰, and notably, regional populations of *L. decemlineata* have demonstrated the ability to independently evolve resistance to pesticides and to do so at different rates³⁷. Previous studies have identified target site mutations in resistance phenotypes and a wide range of genes involved in metabolic detoxification, including carboxylesterase genes, cytochrome P450s, and glutathione S-transferase genes³⁸⁻⁴².

To examine evidence of rapid evolutionary change underlying *L. decemlineata*'s extraordinary success utilizing novel host plants, climates, and detoxifying insecticides, we evaluated structural and functional genomic changes relative to other beetle species, using whole-genome sequencing, transcriptome sequencing, and a large community-driven biocuration effort to improve predicted gene annotations. We compared the size of gene families associated with particular traits against existing available genomes from related species, particularly those sequenced by the i5k project (http://i5k.github.io), an initiative to sequence 5,000 species of Arthropods. While efforts have been made to understand the genetic basis of phenotypes in *L. decemlineata* (for example, pesticide resistance)³²,⁴³,⁴⁴, previous work has been limited to candidate gene approaches rather than comparative genomics. Genomic data can not only illuminate the genetic architecture of a number of phenotypic traits that enable *L. decemlineata* to continue to be an agricultural pest, but can also be used to identify new gene targets for control measures. For example, recent efforts have been made to develop RNAi-based pesticides targeting critical metabolic pathways in *L. decemlineata*¹⁴,⁴⁵,⁴⁶. With the extensive wealth of biological knowledge and a newly-released genome, this beetle is well-positioned to be a model system for agricultural pest genomics and the study of rapid evolution.

## Results and Discussion

### Genome Assembly, Annotation and Assessment.
A single female *L. decemlineata* from Long Island, NY, USA, a population known to be resistant to a wide range of insecticides⁴⁷,⁴⁸, was sequenced to a depth of ~140x coverage and assembled with ALLPATHS⁴⁹ followed by assembly improvement with ATLAS (https://www.hgsc.bcm.edu/software/). The average coleopteran genome size is 760 Mb (ranging from 160–5,020 Mb⁵⁰), while most of the beetle genome assemblies have been smaller (mean assembly size 286 Mb, range 160–710 Mb)⁵¹⁻⁵⁵. The draft genome assembly of *L. decemlineata* is 1.17 Gb and consists of 24,393 scaffolds, with a N50 of 414 kb and a contig N50 of 4.9 kb. This assembly is more than twice the estimated genome size of 460 Mb⁵⁶, with the presence of gaps comprising 492 Mb, or 42%, of the assembly. As this size might be driven by underlying heterozygosity, we also performed scaffolding with REDUNDANS⁵⁷, which reduced the assembly size to 642 Mb, with gaps reduced to 1.3% of the assembly. However, the REDUNDANS assembly increased the contig N50 to 47.4 kb, the number of scaffolds increased to 90,205 and the N50 declined to 139 kb. By counting unique 19 bp kmers and adjusting for ploidy, we estimate the genome size as 816.9 Mb. Using just the small-insert (180 bp) 100 bp PE reads, average coverage was 24X for the ALLPATHS assembly and 26X for the REDUNDANS assembly. For all downstream analyses, the ALLPATHS assembly was used due to its increased scaffold length and reduced number of scaffolds.

The number of genes in the *L. decemlineata* genome predicted based on automated annotation using MAKER was 24,671 gene transcripts, with 93,782 predicted exons, which surpasses the 13,526–22,253 gene models reported in other beetle genome projects (Fig. 1)⁵¹⁻⁵⁵. This may be in part due to fragmentation of the genome, which is known to inflate gene number estimates⁵⁸. To improve our gene models, we manually annotated genes using expert opinion and additional mRNA resources (see Supplementary File for more details). A total of 1,364 genes were manually curated and merged with the unedited MAKER annotations to produce an official gene set (OGS v1.0) of 24,850 transcripts, comprised of 94,859 exons. A total of 12 models were curated as pseudogenes. A total of 1,237 putative transcription factors (TFs) were identified in the *L. decemlineata* predicted proteome (Fig. 2). The predicted number of TFs is similar to some beetles, such as *Anoplophora glabripennis* (1,397) and *Hypothenemus hampei* (1,148), but substantially greater than others, such as *Tribolium castaneum* (788), *Nicrophorus vespilloides* (744), and *Dendroctonus ponderosae* (683)⁵¹⁻⁵⁵.

We assessed the completeness of both the ALLPATHS and REDUNDANS assemblies, and the OGS separately, using benchmarking sets of universal single-copy orthologs (BUSCOs) based on 35 holometabolous insect genomes⁵⁹, as well as manually assessing the completeness and co-localization of the homeodomain transcription factor gene clusters in the ALLPATHS assembly. Using the reference set of 2,442 BUSCOs, the ALLPATHS genome assembly, REDUNDANS genome assembly, and OGS were 93.0%, 91.9%, and 71.8% complete, respectively. We found an additional 4.1%, 5.4%, and 17.9% of the BUSCOs present but fragmented in each dataset, respectively. For the highly conserved Hox and Iroquois Complex (Iro-C) clusters, we located and annotated
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Ultrametric tree with branch lengths in millions of years for *Leptinotarsa decemlineata* relative to five other Coleoptera genomes. The *L. decemlineata* lineage is shown in orange. Branches are labeled with their length in years (top) and with the number of gene family expansions (blue) and contractions (purple) that occurred on that lineage. Rapid changes for both types are in parentheses.</figcaption></figure>

complete gene models in the ALLPATHS genome assembly for all 12 expected orthologs, but these were split across six different scaffolds (Supplementary Figure 1S and Table 1S). All linked Hox genes occurred in the expected order and with the expected, shared transcriptional orientation, suggesting that the current draft assembly was correct but incomplete (see also Supplementary Figures 2S and 3S). Assuming direct concatenation of scaffolds, the Hox cluster would span a region of 3.7 Mb, similar to the estimated 3.5 Mb Hox cluster of *A. glabripennis*⁵⁴. While otherwise highly conserved with *A. glabripennis*, we found a tandem duplication for *Hox3/zen* and an Antennapedia-class (ANTP-class) homeobox gene with no clear ortholog in other arthropods. We also assessed the ALLPATHS genome assembly for evidence of contamination using a Blobplot (Supplementary Figure 4S), which identified a small proportion of the reads (3.3%) as putative contaminants (the largest single taxonomic group represented, 3.0%, was assigned to Platyhelminthes).

**Gene Annotation, Gene Family Evolution and Differential Expression.** We estimated a phylogeny among six coleopteran genomes (*A. glabripennis*, *Agrilus planipennis*, *D. ponderosae*, *L. decemlineata*, *Onthophagus taurus*, and *T. castaneum*; unpublished genomes available at http://i5k.github.io/)⁵¹,⁵²,⁵⁴ using a conserved set of single copy orthologs and compared the official gene set of each species to understand how gene families evolved the branch representing Chysomelidae. *Leptinotarsa decemlineata* and *A. glabripennis* (Cerambycidae) are sister taxa (Fig. 1), as expected for members of the same superfamily Chrysomeloidea. We found 166 rapidly evolving gene families along the *L. decemlineata* lineage (1.4% of 11,598), 142 of which are rapid expansions and the remaining 24 rapid contractions (Table 1). Among all branches of our coleopteran phylogeny, *L. decemlineata* has the highest average expansion rate (0.203 genes per million years), the highest number of genes gained, and the greatest number of rapidly evolving gene families.

Examination of the functional classification of rapidly evolving families in *L. decemlineata* (Supplementary Tables 2S and 3S) indicates that a subset of families are clearly associated with herbivory. The peptidases, comprising several gene families that play a major role in plant digestion³²,⁶⁰, displayed a significant expansion in genes (OrthoDB family Id: EOG8JDKNM, EOG8GTNV8, EOG8DRCSN, EOG8GTNT0, EOG8K0T47, EOG88973C, EOG854CDR, EOG8Z91BB, EOG80306V, EOG8CZF0X, EOG80P6ND,,EOG8BCH22, EOG8ZKN28, EOG8F4VSG, EOG80306V, EOG8BCH22, EOG8ZKN28, EOG8F4VSG). While olfactory receptor gene families have rapidly contracted (EOG8Q5C4Z, EOG8RZ1DX), subfamilies of odorant binding proteins and gustatory receptors have grown (see Sensory Ecology section below). The expansion of gustatory receptor subfamilies are associated with bitter receptors, likely reflecting host plant detection of nightshades (Solanaceae). Some gene families associated with plant detoxification and insecticide resistance have rapidly expanded (Glutathione S-transferases: EOG85TG3K, EOG85F05D, EOG8BCH22; UDP-glycosyltransferase: EOG8BCH22; cuticle proteins: EOG8QJV4S, EOG8DNHJQ; ABC transporters: EOG83N9TJ), whereas others have contracted (Cytochrome P450s: EOG83N9TJ). Finally, gene families associated with immune defense (fibrinogen: EOG8DNHJQ; Immunoglobulin: EOG8Q87B6, EOG8S7N55, EOG854CDX, EOG8KSRZK, EOG8CNT5M) exhibit expansions that may be linked to defense against pathogens and parasitoids that commonly attack exposed herbivores⁶¹. A substantial proportion of the rapidly evolving gene families include proteins with transposable element domains (25.3%), while other important functional groups include DNA and protein binding (including many transcription factors), nuclease activity, protein processing, and cellular transport.

Diversification of transcription factor (TF) families potentially signals greater complexity of gene regulation, including enhanced cell specificity and refined spatiotemporal signaling⁶². Notably, several TF families
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Heatmap distribution of the abundance of transcription factor families in *Leptinotarsa decemlineata* compared to other insects. Each entry indicates the number of TF genes for the given family per genome, based on presence of predicted DNA binding domains. Color scale is log (base 2) and the key is depicted at the top (light blue means the TF family is completely absent). Families discussed in the main text are indicated by arrows.</figcaption></figure>

are substantially expanded in *L. decemlineata*, including *HTH_psq* (194 genes vs. a mean of only 24 across the insects shown in Fig. 2), *MADF* (152 vs. 54), and *THAP* (65 vs. 41). Two of these TF families, *HTH_psq* and *THAP*, are DNA binding domains in transposons. Of the 1,237 *L. decemlineata* TFs, we could infer DNA binding motifs for 189 (15%) (Supplementary Table 4S), mostly based on DNA binding specificity data from *Drosophila*
<table><caption>**Table 1.** Summary of gene gain and loss events inferred after correcting for annotation and assembly error across six Coleoptera species. The number of rapidly evolving families is shown in parentheses for each type of change and the rate is genes per million years.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Expansions</th><th colspan="3">Contractions</th><th rowspan="2">No Change</th><th rowspan="2">Average Expansion</th></tr><tr><th>Families</th><th>Genes gained</th><th>Rate</th><th>Families</th><th>Genes lost</th><th>Rate</th></tr></thead><tbody><tr><td>Anoplophora glabripennis</td><td>865 (13)</td><td>1231</td><td>1.42</td><td>1988 (107)</td><td>3125</td><td>1.57</td><td>5850</td><td>-0.182341</td></tr><tr><td>Agrilus plannipennis</td><td>739 (100)</td><td>1991</td><td>2.69</td><td>707 (8)</td><td>769</td><td>1.09</td><td>7257</td><td>0.119108</td></tr><tr><td>Dendroctonus ponderosae</td><td>933 (72)</td><td>1982</td><td>2.12</td><td>1606 (21)</td><td>1887</td><td>1.17</td><td>6164</td><td>0.006438</td></tr><tr><td>Leptinotarsa decemlineata</td><td>426 (40)</td><td>855</td><td>2.01</td><td>1556 (48)</td><td>2116</td><td>1.36</td><td>6721</td><td>-0.127501</td></tr><tr><td>Onthophagus taurus</td><td>1299 (142)</td><td>2952</td><td>2.27</td><td>767 (24)</td><td>895</td><td>1.17</td><td>6637</td><td>0.203380</td></tr><tr><td>Tribolium castaneum</td><td>786 (51)</td><td>1428</td><td>1.82</td><td>516 (27)</td><td>909</td><td>1.76</td><td>7401</td><td>0.055645</td></tr></tbody></table>

melanogaster (124 TFs), but also from species as distant as human (45 TFs) and mouse (11 TFs). Motifs were inferred for a substantial proportion of the TFs in the largest TF families, including Homeodomain (59 of 90, 66%), bHLH (34 of 46, 74%), and Forkhead box (14 of 16, 88%). We could only infer a small number of C2H2 zinc finger motifs (21 of 439, ~5%), which is expected as these sequences evolve quickly by shuffling zinc finger arrays, resulting in largely dissimilar DNA-binding domains across metazoans⁶³. Collectively, the almost 200 inferred DNA binding motifs for *L. decemlineata* TFs provide a unique resource to begin unraveling gene regulatory networks in this organism.

To identify genes active in mid-gut tissues, life-stages, and sex differences, we examined differential transcript expression levels using RNA sequencing data. Comparison of significantly differentially expressed genes with >100-fold change, after Bonferroni correction, indicated higher expression of digestive enzymes (proteases, peptidases, dehydrogenases and transporters) in mid-gut versus whole larval tissues, while cuticular proteins were largely expressed at lower levels (Fig. 3A, Supplementary Table 5S). Comparison of an adult male and female showed higher expression of testes and sperm related genes in males, while genes involved in egg production and sterol biosynthesis are more highly expressed in females (Fig. 3B, Supplementary Table 6S). Comparisons of larvae to both an adult male (Fig. 3C, Supplementary Table 7S) and an adult female (Fig. 3D, Supplementary Table 8S) showed higher expression of larval-specific cuticle proteins, and lower expression of odorant binding and chemosensory proteins. The adults, both drawn from a pesticide resistant population, showed higher constitutive expression of cytochrome p450 genes compared to the larval population, which is consistent with the results from previous studies of neonicotinoid resistance in this population⁴⁸.

**Transposable Elements.** Transposable elements (TEs) are ubiquitous mobile elements within most eukaryotic genomes and play critical roles in both genome architecture and the generation of genetic variation⁶⁴. Through insertional mutagenesis and recombination, TEs are a major contributor to the generation of novel mutations (50–80% of all mutation events within the genome of *D. melanogaster*)⁶⁵⁻⁶⁷, and are increasingly thought to generate much of the genetic diversity that contributes to rapid evolution⁶⁸⁻⁷⁰. In addition to finding that genes with TE domains comprise 25% of the rapidly evolving gene families within *L. decemlineata*, we found that at least 17% of the genome consists of TEs (Supplementary Table 9S). This is substantially greater than the 6% found in *T. castaneum*⁵¹, but less than some Lepidoptera (35% in *Bombyx mori* and 25% in *Heliconius mel-pomene*)⁷¹,⁷². LINEs were the largest TE class, comprising ~10% of the genome, while SINEs were not detected. Curation of the TE models with intact protein domains resulted in 334 current models of potentially active TEs, meaning that these TEs are capable of transposition and excision. Within the group of active TEs, we found 191 LINEs, 99 DNA elements, 38 LTRs, and 5 Helitrons. Given that TEs have been associated with the ability of species to rapidly adapt to novel selection pressures⁷³⁻⁷⁵, particularly via alterations of gene expression patterns in neighboring genomic regions, we scanned gene rich regions (1 kb neighborhood size) for active TE elements. Genes with active neighboring TEs have functions that include transport, protein digestion, diapause, and metabolic detoxification (Supplementary Table 10S). Because TE elements have been implicated in conferring insecticide resistance in other insects⁷⁶, future work should investigate the role of these TE insertions on rapid evolutionary changes within pest populations of *L. decemlineata*.

**Population Genetic Variation and Invasion History.** To understand the propensity for *L. decemlineata* pest populations to rapidly evolve across a range of environmental conditions, we examined geographical patterns of genomic variability and the evolutionary history of *L. decemlineata* (to see the genomic analysis of environmental adaptation, refer to the Diapause and Environmental Stress section of Supplementary File). High levels of standing variation are one mechanism for rapid evolutionary change⁷⁷. We identified 1.34 million biallelic single nucleotide polymorphisms (SNPs) from pooled RNAseq datasets, or roughly 1 variable site for every 22 base pairs of coding DNA. This rate of polymorphism is exceptionally high when compared to vertebrates (e.g. ~1 per kb in humans, or ~1 per 500 bp in chickens)⁷⁸,⁷⁹, and is 8-fold higher than other beetles (1 in 168 for *D. ponderosae* and 1 in 176 bp for *O. taurus*)⁵²,⁸⁰ and 2 to 5-fold higher than some dipterans (1 in 54 bp for *D. melanogaster* and 1 in 125 bp for *Anopheles gambiae*)⁷⁸,⁸¹. It is likely that these values simply scale with effective population size, although the dipterans, with the largest known population sizes, have reduced variation due to widespread selective sweeps and genetic bottlenecks⁸².
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Volcano plots showing statistically significant gene expression differences after Bonferroni correction in *Leptinotarsa decemlineata*. (A) Mid-gut tissue versus whole larvae, (B) an adult male versus an adult female, (C) an adult male versus whole larvae, and (D) an adult female versus whole larvae. Points outside the gray area indicate >100-fold-differences in expression. Blue points indicate down-regulated genes and red points indicate up-regulated genes in each contrast.</figcaption></figure>

Evolutionary relationships and the amount of genetic drift among Midwestern USA, Northeastern USA, and European *L. decemlineata* populations were estimated based on genome-wide allele frequency differences using a population graph. A substantial amount of local genetic structure and high genetic drift is evident among all populations, although both the reference lab strain from New Jersey and European populations appear to have undergone more substantial drift, suggestive of strong inbreeding (Fig. 4). Population genetic divergence (*F*<sub>ST</sub>) values (Supplementary Table 11S) range from 0.035 (among Wisconsin populations) to 0.182 (Wisconsin compared to Europe). The allele frequency spectrum was calculated for populations in Wisconsin, Michigan, and Europe to estimate the population genetic parameter *θ*, or the product of the mutation rate and the ancestral effective population size, and the ratio of contemporary to ancestral population size in models that allowed for single or multiple episodes of population size change. Estimates of ancestral *θ* are much higher for Wisconsin (*θ* = 12595) and Michigan (*θ* = 93956) than Europe (*θ* = 3.07; Supplementary Table 12S), providing support for a single introduction into Europe following a large genetic bottleneck<sup>15</sup>. In all three populations, a model of population size growth is supported, in agreement with historical accounts of the beetles expanding from the Great Plains into the Midwestern U.S. and Europe<sup>3,15</sup>, but the dynamics of each population appear independent, with the population from Michigan apparently undergoing a very recent decline in contemporary population size (the ratio of contemporary to ancient population size is 0.066, compared to 3.3 and 2.1 in Wisconsin and Europe, respectively).

**Sensory Ecology and Host Plant Detection.** To interact with their environment, insects have evolved neurosensory organs to detect environmental signals, including tactile, auditory, chemical and visual cues<sup>83</sup>. We examined neural receptors, olfactory genes, and light sensory (opsin) genes to understand the sensory ecology and host-plant specializations of *L. decemlineata*.

We found high sequence similarity in the neuroreceptors of *L. decemlineata* compared to other Coleoptera. The transient receptor potential (TRP) channels are permeable transmembrane proteins that respond to
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Population genetic relationships and relative rates of genetic drift among *Leptinotarsa decemlineata* pest populations based on single nucleotide polymorphism data. Population codes: NJ- New Jersey lab strain, WIs- imidacloprid susceptible population from Arlington, Wisconsin, WIr- imidacloprid resistant population from Hancock, Wisconsin, MI- imidacloprid resistant population from Michigan, and EU- European samples combined from Italy and Russia.</figcaption></figure>

<table><caption>**Table 2.** Numbers of putatively functional proteins in four chemosensory families in *Leptinotarsa decemlineata* and other beetle species.</caption><thead><tr><td>Species</td><td>Odorant Binding Proteins (OBPs)</td><td>Olfactory Receptors (ORs)</td><td>Gustatory Receptors (GRs)</td><td>Ionotropic Receptors (IRs)</td></tr></thead><tbody><tr><td>Leptinotarsa decemlineata</td><td>58</td><td>75</td><td>144</td><td>27</td></tr><tr><td>Anoplophora glabripennis</td><td>52</td><td>120</td><td>234</td><td>72</td></tr><tr><td>Tribolium castaneum</td><td>49</td><td>264</td><td>219</td><td>72</td></tr></tbody></table>

temperature, touch, pain, osmolarity, pheromones, taste, hearing, smell and visual cues of insects<sup>84,85</sup>. In most insect genomes, there are typically 13–14 TRP genes located in insect stretch receptor cells and several are targeted by commercial insecticides<sup>86</sup>. We found 12 TRP genes present in the *L. decemlineata* genome, including the two TRPs (Nanchung and Inactive) that are targeted by commercial insecticides, representing a complete set of one-to-one orthologs with *T. castaneum*. Similarly, the 20 known amine neurotransmitter receptors in *T. castaneum* are present as one-to-one orthologs in *L. decemlineata*<sup>87,88</sup>. Amine receptors are G-protein-coupled receptors that interact with biogenic amines, such as octopamine, dopamine and serotonin. These neuroactive substances regulate behavioral and physiological traits in animals by acting as neurotransmitters, neuromodulators and neurohormones in the central and peripheral nervous systems<sup>89</sup>.

The majority of phytophagous insects are restricted to feeding on several plant species within a genus, or at least restricted to a particular plant family<sup>90</sup>. Thus, to find their host plants within heterogeneous landscapes, insect herbivores detect volatile organic compounds through olfaction, which utilizes several families of chemosensory gene families, such as the odorant binding proteins (OBPs), odorant receptors (ORs), gustatory receptors (GRs), and ionotropic receptors (IRs)<sup>91</sup>. OBPs directly bind with volatile organic compounds emitted from host plants and transport the ligands across the sensillar lymph to activate the membrane-bound ORs in the dendrites of the olfactory sensory neurons<sup>92</sup>. The ORs and GRs are 7-transmembrane proteins related within a superfamily<sup>93</sup> of ligand-gated ion channels<sup>94</sup>. The ionotropic receptors are related to ionotropic glutamate receptors and function in both smell and taste<sup>95</sup>. These four gene families are commonly large in insect genomes, consisting of tens to hundreds of members.

We compared the number of genes found in *L. decemlineata* in the four chemosensory gene families to *T. castaneum* and *A. glabripennis* (Table 2, as well as Supplementary Tables 13S–16S for details). While the OBP family is slightly enlarged, the three receptor families are considerably smaller in *L. decemlineata* than in either *A. glabripennis* or *T. castaneum*, consistent with the specialization of this beetle on one genus of plants. However, each beetle species exhibits species-specific gene subfamily expansions (Supplementary Figures 5S–8S); in particular, some lineages of GRs related to bitter taste are expanded in *L. decemlineata* relative to *A. glabripennis* and other beetles. Among the OBPs, we identified a major *L. decemlineata*-specific expansion of proteins belonging to the Minus-C class (OBPs that have lost two of their six conserved cysteine residues) that appear unrelated to the 'traditional' Minus-C subfamily in Coleoptera, indicating that coleopteran OBPs have lost cysteines on at least two occasions.

To understand the visual acuity of *L. decemlineata*, we examined the G-protein-coupled transmembrane receptor opsin gene family. We found five opsins, three of which are members of rhabdomeric opsin (R-opsin) subfamilies expressed in the retina of insects<sup>51,96</sup>. Specifically, the *L. decemlineata* genome contains one member of the long wavelength-sensitive R-opsin and two short wavelength UV-sensitive R-opsins. The latter were found
to be closely linked in a range of less than 20,000 bp, suggestive of recent tandem gene duplication. Overall, the recovered repertoire of retinally-expressed opsins in *L. decemlineata*⁹⁷ is consistent with the beetle's attraction to yellow light and to the yellow flowers of its ancestral host plant, *S. rostratum*⁹⁸,⁹⁹, and is consistent with the beetle's sensitivity in the UV- and LW-range¹⁰⁰. In addition, we found a member of the Rh7 R-opsin subfamily, which is broadly conserved in insects including other beetle species (*A. glabripennis*), although it is missing from *T. castaneum*. Finally, *L. decemlineata* has a single ortholog of the c-opsin subfamily shared with *T. castaneum*, which is absent in *A. glabripennis* and has an unclear role in photoreception¹⁰¹.

**Host Plant Utilization.** *Protein digestion.* Insect herbivores are fundamentally limited by nitrogen availability¹⁰², and thus need to efficiently break down plant proteins in order to survive and develop on host plants¹⁰³. *Leptinotarsa decemlineata* has serine and cysteine digestive peptidases (coined “intestains”)³²,¹⁰⁴,¹⁰⁵, as well as aspartic and metallo peptidases³³, for protein digestion. For the vast majority of plant-eating beetles (the infraorder Cucujiformia, which includes Chrysomelidae), cysteine peptidases contribute most strongly to proteolytic activity in the gut¹⁰⁵,¹⁰⁶. In response to herbivory, plants produce a wide range of proteinase inhibitors to prevent insect herbivores from digesting plant proteins¹⁰⁵,¹⁰⁵,¹⁰⁷. Coleopteran peptidases are differentially susceptible to plant peptidase inhibitors, and our annotation results suggest that gene duplication and selection for inhibitor insensitive genotypes may have contributed to the success of leaf-feeding beetles (Chrysomelidae) on different plants. We found that gene expansion of cysteine cathepsins from the C1 family in *L. decemlineata* correlates with the acquisition of greater digestive function by this group of peptidases, which is supported by gene expression activity of these genes in mid-gut tissue (Fig. 3A, Fig. 5). The gene expansion may be explained by an evolutionary arms race between insects and plants that favors insects with a variety of digestive peptidases in order to resist plant peptidase inhibitors¹⁰⁸,¹⁰⁹ and allows for functional specialization¹¹⁰.

Cysteine peptidases of the C1 family were represented by more than 50 genes separated into four groups with different structure and functional characteristics (Supplementary Table 18S): cathepsin L subfamily, cathepsin B subfamily, TINAL-like genes, and cysteine peptidase inhibitor domains (CPIDs). Cathepsin L subfamily cysteine peptidases are endopeptidases¹¹¹ that can be distinguished by the cathepsin propeptide inhibitor domain I29 (pfam08246)¹¹²,¹¹³. Within the cathepsin L subfamily, we found sequences that were similar to classical cathepsin L, cathepsin F, and cathepsin O. However, there were 28 additional predicted peptidases of this subfamily that could not be assigned to any of the “classical” cathepsin types, and most of these were grouped into two gene expansions (uL1 and uL2) according to their phylogenetic and structural characteristics. Cathepsin B subfamily cysteine peptidases are distinguished by the specific peptidase family C1 propeptide domain (pfam08127). Within the cathepsin B subfamily, there was one gene corresponding to typical cathepsin B peptidases and 14 cathepsin B-like genes. According to the structure of the occluding loop, only the typical cathepsin B may have typical endo- and exopeptidase activities, while a large proportion of cathepsin B-like peptidases presumably possesses only endopeptidase activity due to the absence of a His-His active subsite in the occluding loop, which is responsible for exopeptidase activity¹¹¹. Only one gene corresponding to a TINAL-like-protein was present, which has a domain similar to cathepsin B in the C-terminus, but lacks peptidase activity due to the replacement of the active site Cys residue with Ser¹¹⁴. Cysteine peptidase inhibitor domain (CPID) genes encode the I29 domain of cysteine peptidases without the mature peptidase domain. Within the CPID group, there were seven short inhibitor genes that lack the enzymatic portion of the protein. A similar trend of “stand-alone inhibitors” has been observed in other insects, such as *B. mori*¹¹⁵. These CPID genes may be involved in the regulation of cysteine peptidases. We note that we found multiple fragments of cysteine peptidase genes, suggesting that the current list of *L. decemlineata* genes may be incomplete. Comparison of these findings with previous data on *L. decemlineata* cysteine peptidases¹¹⁶ demonstrates that intestains correspond to several peptidase genes from the uL1 and uL2 groups (Supplementary Table 18S). These data, as well as literature for Tenebrionidae beetles¹¹⁷, suggest that intensive gene expansion is typical for peptidases that are involved in digestion.

We also found a high number of digestion-related serine peptidase genes in the *L. decemlineata* genome (Supplementary Table 19S), but they contribute only a small proportion of the beetle's total gut proteolytic activity³². Of the 31 identified serine peptidase genes and fragments, we annotated 16 as trypsin-like peptidases and 15 as chymotrypsin-like peptidases. For four chymotrypsin-like and one trypsin-like peptidase, we identified only short fragments. All complete (and near-complete) sequences have distinctive S1A peptidase subfamily motifs, a conserved catalytic triad, conserved sequence residues such as the “CWC” sequence and cysteines that form disulfide bonds in the chymotrypsin protease fold. The number of serine peptidases was higher than expected based upon the number of previously identified EST clones³², but lower than the number of chymotrypsin and trypsin genes in the *T. castaneum* genome.

**Carbohydrate digestion.** Carbohydrates are the other category of essential nutrients for *L. decemlineata*. The enzymes that assemble and degrade oligo- and polysaccharides, collectively termed Carbohydrate active enzymes (CAZy), are categorized into five major classes: glycoside hydrolases (GH), polysaccharide lyases, carbohydrate esterases (CE), glycosyltransferases (GT) and various auxiliary oxidative enzymes¹¹⁸. Due to the many different roles of carbohydrates, the CAZy family profile of an organism can provide insight into “glycobiological potential” and, in particular, mechanisms of carbon acquisition¹¹⁹. We identified 182 GHs assigned to 25 families, 181 GTs assigned to 41 families, and two CEs assigned to two families in *L. decemlineata*; additionally, 99 carbohydrate-binding modules (which are non-catalytic modules associated with the above enzyme classes) were present and assigned to 9 families (Supplementary Table 20S; the list of CAZy genes is presented in Supplementary Table 21S). We found that *L. decemlineata* has three families of genes associated with plant cell wall carbohydrate digestion (GH28, GH45 and GH48) that commonly contain enzymes that target pectin (GH28) and cellulose (GH45 and GH48), the major structural components of leaves¹²⁰. We found evidence of massive gene duplications in the GH28 family (14 genes) and GH45 family (11 genes, plus one additional splicing variant), whereas GH48 is represented by only three genes in the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Phylogenetic relationships of the cysteine peptidase gene family in *Leptinotarsa decemlineata* compared to model insects. Species abbreviations are: *L. decemlineata* (Ld, green color), *Drosophila melanogaster* (Dm, blue color), *Apis mellifera* (Am, purple color), and *Tribolium castaneum* (Tc, red color). Mid-gut gene expression (TPM) of highly expressed *L. decemlineata* cysteine peptidases is shown as bar graphs across three replicate treatments.</figcaption></figure>

genome¹²¹. Overall, the genome of *L. decemlineata* shows a CAZy profile adapted to metabolize pectin and cellulose contained in leaf cell walls. The absence of specific members of the families GH43 (α-L-arabinofuranosidases) and GH78 (α-L-rhamnosidases) suggests that *L. decemlineata* can break down homogalacturonan, but not substituted galacturonans such as rhamnogalacturonan I or II¹²⁰. The acquisition of these plant cell wall degrading enzymes has been linked to horizontal transfer in the leaf beetles and other phytophagous beetles¹²⁰, with strong phylogenetic evidence supporting the transfer of GH28 genes from a fungal donor (Pezizomycotina) in *L. decemlineata*, as well as in the beetles *D. ponderosae* and *Hypothenemus hampei*, but a novel fungal donor in the more closely related cerambycid beetle *A. glabripennis*⁵⁴ and a bacterial donor in the weevil *Callosobruchus maculatus*¹²⁰.

**Insecticide Resistance.** To understand the functional genomic properties of insecticide resistance, we examined genes important to neuromuscular target site sensitivity, tissue penetration, and prominent gene families involved in Phase I, II, and III metabolic detoxification of xenobiotics¹²². These include the cation-gated
nicotinic acetylcholine receptors (nAChRs), the γ-amino butyric acid (GABA)-gated anion channels and the histamine-gated chloride channels (HisClS), cuticular proteins, cytochrome P450 monooxygenases (CYPs), and the Glutathione S-transferases (GSTs).

Many of the major classes of insecticides (organochlorides, organophosphates, carbamates, pyrethroids, neonicotinoids and ryanoids) disrupt the nervous system (particularly ion channel signaling), causing paralysis and death<sup>123</sup>. Resistance to insecticides can come from point mutations that reduce the affinity of insecticidal toxins to ligand-gated ion superfamily genes<sup>124</sup>. The cys-loop ligand-gated ion channel gene superfamily is comprised of receptors involved in mediating synaptic ion flow during neurotransmission<sup>125</sup>. A total of 22 cys-loop ligand-gated ion channels were identified in the *L. decemlineata* genome in numbers similar to those observed in other insects<sup>126</sup>, including 12 nAChRs, three GABA receptors, and two HisClS (Supplementary Figure 9S). The GABA-gated chloride channel homolog of the *Resistance to dieldrin (Rdl)* gene of *D. melanogaster* was examined due to its role in resistance to dieldrin and other cyclodienes in *Diptera*<sup>127</sup>. The coding sequence is organized into 10 exons (compared to nine in *D. melanogaster*) on a single scaffold, with duplications of the third and sixth exon (Supplementary Figure 10S). Alternative splicing of these two exons encodes for four different polypeptides in *D. melanogaster*<sup>128,129</sup>, and as the splice junctions are present in *L. decemlineata*, we expect the same diversity of *Rdl*. The point mutations in the transmembrane regions TM2 and TM3 of *Rdl* are known to cause insecticide resistance in *Diptera*<sup>124,127</sup>, but were not observed in *L. decemlineata*.

Cuticle genes have been implicated in imidacloprid resistant *L. decemlineata*<sup>130</sup> and at least one has been shown to have phenotypic effects on resistance traits following RNAi knockdown<sup>131</sup>. A total of 163 putative cuticle protein genes were identified and assigned to one of seven families (CPR, CPAP1, CPAP3, CPF, CPCFC, CPLCG, and TWDL) (Supplemental Figure 11S and Table 22S). Similar to other insects, the CPR family, with the RR-1 (soft cuticle), RR-2 (hard cuticle), and unclassifiable types, constituted the largest group of cuticle protein genes (132) in the *L. decemlineata* genome. While the number of genes in *L. decemlineata* is slightly higher than in *T. castaneum* (110), it is similar to *D. melanogaster* (137)<sup>132</sup>. Numbers in the CPAP1, CPAP3, CPF, and TWDL families were similar to other insects, and notably no genes with the conserved sequences for CPLCA were detected in *L. decemlineata*, although they are found in other Coleoptera.

A total of 89 CYP (P450) genes were identified in the *L. decemlineata* genome, an overall decrease relative to *T. castaneum* (143 genes). Due to their role in insecticide resistance in *L. decemlineata* and other insects<sup>38,48,133</sup>, we examined the CYP6 and CYP12 families in particular. Relative to *T. castaneum*, we observed reductions in the CYP6BQ, CYP4BN, and CYP4Q subfamilies. However, five new subfamilies (CYP6BJ, CYP6BU, CYP6E, CYP6F and CYP6K) were identified in *L. decemlineata* that were absent in *T. castaneum*, and the CYP12 family contains three genes as opposed to one gene in *T. castaneum* (CYP12h1). We found several additional CYP genes not present in *T. castaneum*, including CYP413A1, CYP421A1, CYP4V2, CYP12J and CYP12J4. Genes in CYP4, CYP6, and CYP9 are known to be involved in detoxification of plant allelochemicals as well as resistance to pesticides through their constitutive overexpression and/or inducible expression in imidacloprid resistant *L. decemlineata*<sup>48,130</sup>.

GSTs have been implicated in resistance to organophosphate, organochlorine and pyrethroid insecticides<sup>134</sup> and are responsive to insecticide treatments in *L. decemlineata*<sup>130,135</sup>. A total of 27 GSTs were present in the *L. decemlineata* genome, and while they represent an expansion relative to *A. glabripennis*, all have corresponding homologs in *T. castaneum*. The cytosolic GSTs include the epsilon (11 genes), delta (5 genes), omega (4 genes), theta (2 genes) and sigma (3 genes) families, while two GSTs are microsomal (Supplementary Figure 12S). Several GST-like genes present in the *L. decemlineata* genome represent the Z class previously identified using transcriptome data<sup>135</sup>.

**Pest Control via the RNAi Pathway.** RNA interference (RNAi) is the process by which small non-coding RNAs trigger sequence-specific gene silencing<sup>136</sup>, and is important in protecting against viruses and mobile genetic elements, as well as regulating gene expression during cellular development<sup>137</sup>. The application of exogenous double stranded RNA (dsRNA) has been exploited as a tool to suppress gene expression for functional genetic studies<sup>138,139</sup> and for pest control<sup>41,140</sup>.

We annotated a total of 49 genes associated with RNA interference, most of them were found on a single scaffold. All genes from the core RNAi machinery (from all three major RNAi classes) were present in *L. decemlineata*, including fifteen genes encoding components of the RNA Induced Silencing Complex (RISC) and genes known to be involved in double-stranded RNA uptake, transport, and degradation (Supplementary Table 22S). A complete gene model was annotated for *R2D2*, an essential component of the siRNA pathway that interacts with *dicer-2* to load siRNAs into the RISC, and not previously detected in the transcriptome of the *L. decemlineata* mid-gut<sup>141</sup>. The core components of the small interfering RNA (siRNA) pathway were duplicated, including *dicer-2*, an RNase III enzyme that cleaves dsRNAs and pre-miRNAs into siRNAs and miRNAs respectively<sup>142,143</sup>. The *dicer-2a* and *dicer-2b* CDS have 60% nucleotide identity to each other, and 56% and 54% identity to the *T. castaneum dicer-2* homolog, respectively. The *argonaute-2* gene, which plays a key role in RISC by binding small non-coding RNAs, was also duplicated. A detailed analysis of these genes will be necessary to determine if the duplications provide functional redundancy. The duplication of genes in the siRNA pathway may play a role in the high sensitivity of *L. decemlineata* to RNAi knockdown<sup>144</sup> and could benefit future efforts to develop RNAi as a pest management technology.

Conclusion

The whole-genome sequence of *L. decemlineata*, provides novel insights into one of the most diverse animal taxa, Chrysomelidae. It is amongst the largest beetle genomes sequenced to date, with a minimum assembly size of 640 Mb (ranging up to 1.17 Gb) and 24,740 genes. The genome size is driven in part by a large number of transposable element families, which comprise at least 17% of the genome and appear to be rapidly expanding relative to
other beetles. Population genetic analyses suggest high levels of nucleotide diversity, local geographic structure, and evidence of recent population growth, which helps to explain how *L. decemlineata* rapidly evolves to exploit novel host plants, climate space, and overcome a range of pest management practices (including a large and diverse number of insecticides). Digestive enzymes, in particular the cysteine peptidases and carbohydrate-active enzymes, show evidence of gene expansion and elevated expression in gut tissues, suggesting the diversity of the genes is a key trait in the beetle's phytophagous lifestyle. Additionally, expansions of the gustatory receptor sub-family for bitter tasting might be a key adaptation to exploiting hosts in the nightshade family, Solanaceae, while expansions of novel subfamilies of CYP and GST proteins are consistent with rapid, lineage-specific turnover of genes implicated in *L. decemlineata*'s capacity for insecticide resistance. Finally, *L. decemlineata* has interesting duplications in RNAi genes that might increase its sensitivity to RNAi and provide a promising new avenue for pesticide development. The *L. decemlineata* genome promises new opportunities to investigate the ecology, evolution, and management of this species, and to leverage genomic technologies in developing sustainable methods of pest control.

## Materials and Methods

**Genome Characteristics and Sampling of DNA and RNA.** Previous cytological work determined that *L. decemlineata* is diploid and consists of 34 autosomes plus an XO system in males, or an XX system in females¹⁴⁵. Twelve chromosomes are submetacentric, while three areacrocentric and two are metacentric, although one chromosome is heteromorphic (acrocentric and/or metacentric) in pest populations¹⁴⁶. The genome size has been estimated with Feulgen densitometry at 0.46 pg, or approximately 460 Mb⁴⁷. To generate a reference genome sequence, DNA was obtained from a single adult female, sampled from an imidacloprid resistant strain developed from insects collected from a potato field in Long Island, NY. Additionally, whole-body RNA was extracted for one male and one female from the same imidacloprid resistant strain. Raw RNAseq reads for 8 different populations were obtained from previous experiments: two Wisconsin populations (PRJNA297027)¹³⁰, a Michigan population (PRJNA400685), a lab strain originating from a New Jersey field population (PRJNA275431), and three samples from European populations (PRJNA79581 and PRJNA236637)¹⁴,¹²¹. All RNAseq data came from pooled populations or were combined into a population sample from individual reads. In addition, RNA samples of an adult male and female from the same New Jersey population were sequenced separately using Illumina HiSeq 2000 as 100 bp paired end reads (deposited in the GenBank/EMBL/DDBJ database, PRJNA275662), and three samples from the mid-gut of 4th-instar larvae were sequenced using SOLiD 5500 Genetic Analyzer as 50 bp single end reads (PRJNA400633).

**Genome Sequencing, Assembly, Annotation and Assessment.** Four Illumina sequencing libraries were prepared, with insert sizes of 180 bp, 500 bp, 3 kb, and 8 kb, and sequenced with 100 bp paired-end reads on the Illumina HiSeq 2000 platform at estimated 40x coverage, except for the 8kb library, which was sequenced at estimated 20x coverage. ALLPATHS-LG v35218⁴⁹ was used to assemble reads. Two approaches were used to scaffold contigs and close gaps in the genome assembly. The reference genome used in downstream analyses was generated with ATLAS-LINK v1.0 and ATLAS GAP-FILL v2.2 (https://www.hgsc.bcm.edu/software/). In the second approach, REDUNDANS was used⁵⁷, as it is optimized to deal with heterozygous samples. The raw sequence data and *L. decemlineata* genome have been deposited in the GenBank/EMBL/DDBJ database (Bioproject accession PRJNA171749, Genome assembly GCA_000500325.1, Sequence Read Archive accessions: SRX396450-SRX396453). This data can be visualized, along with gene models and supporting data, at the i5k Workspace@NAL: https://i5k.nal.usda.gov/Leptinotarsa_decemlineata and https://apollo.nal.usda.gov/lepdec/jbrowse/¹⁴⁷. We estimated the size of our genome using a kmer distribution plot in JELLYFISH¹⁴⁸, where we mapped the 100 bp paired-end reads from the 180 bp insert library, used a 19 bp kmer distribution plot, and corrected for ploidy.

Automated gene prediction and annotation were performed using MAKER v2.0¹⁴⁹, using RNAseq evidence and arthropod protein databases. The MAKER predictions formed the basis of the first official gene set (OGS v0.5.3). To improve the structural and functional annotation of genes, these gene predictions were manually and collaboratively edited using the interactive curation software Apollo¹⁵⁰. For a given gene family, known insect genes were obtained from model species, especially *T. castaneum*⁵¹ and *D. melanogaster*¹⁵¹, and the nucleotide or amino acid sequences were used in BLASTx or tBLASTn¹⁵²,¹⁵³ to search the *L. decemlineata* OGS v0.5.3 or genome assembly, respectively, on the i5k Workspace@NAL. All available evidence (AUGUSTUS, SNAP, RNA data, etc.), including additional RNAseq data not used in the MAKER predictions, were used to inspect and modify gene predictions. Changes were tracked to ensure quality control. Gene models were inspected for quality, incorrect splits and merges, internal stop codons, and gff3 formatting errors, and finally merged with the MAKER-predicted gene set to produce the official gene set (OGS v1.0; merging scripts are available upon request). For focal gene families (e.g. peptidase genes, odorant and gustatory receptors, RNAi genes, etc.), details on how genes were identified and assigned names based on functional predictions or evolutionary relationship to known reference genes are provided in the Supplementary Material (Supplementary File).

To assess the quality of our genome assemblies, we used BUSCO v2.0⁵⁹ to determine the completeness of each genome assembly and the official gene set (OGS v1.0), separately. We benchmarked our data against 35 insect species in the Endopterygota obd9 database, which consists of 2,442 single-copy orthologs (BUSCOs). Secondly, we annotated and examined the genomic architecture of the Hox and Iroquois Complex gene clusters. For this, tBLASTn searches were performed against the genome using orthologous Hox gene protein sequences from *T. castaneum* (Tcas3.0) and *A. glabripennis*. Provisional *L. decemlineata* models were refined, and potential gene duplications were identified, via iterative and reciprocal BLAST and by manual inspection and correction of protein alignments generated with ClustalW²¹⁵⁴, using RNAseq expression evidence when available. Finally, we used BlobTools v1.0¹⁵⁵ to assess the assembled genome for possible contamination by generating a Blobplot. This plot integrates guanine and cytosine (GC) content of sequences, read coverage and sequence similarity via
blast searches to assess genome contamination. Putative contaminants were identified using the combined results ('bestsumorder' command) of BLASTn against NCBI's nt database¹⁵⁶ and DIAMOND BLAST¹⁵⁷ against the UniProt reference database¹⁵⁸, and read coverage was assessed by mapping the 100 bp paired-end reads from the 180 bp insert library to the ALLPATHS genome with the Burrows-Wheeler Aligner (bwa-mem) version 0.7.5a¹⁵⁹. As these sequences are nested within the genome scaffolds, and could represent horizontally transferred DNA, we left the scaffolds intact and provide a supplementary dataset of the contaminant sequences.

**Gene Family Evolution.** In order to identify rapidly evolving gene families along the *L. decemlineata* lineage, we obtained ~38,000 ortho-groups from 72 Arthropod species as part of the i5k pilot project (Thomas *et al.* in prep) from OrthoDB version 8¹⁶⁰. For each ortho-group, we took only those genes present in the order Coleoptera, which was represented by the following six species: *A. glabripennis*, *A. planipennis*, *D. ponderosae*, *L. decemlineata*, *O. taurus*, and *T. castaneum* (http://i5k.github.io/)⁵¹,⁵²,⁵⁴. Finally, in order to make accurate inferences of ancestral states, families that were present in only one of the six species were removed. This resulted in a final count of 11,598 gene families that, among these six species, form the comparative framework that allowed us to examine rapidly evolving gene families in the *L. decemlineata* lineage.

Aside from the gene family count data, an ultrametric phylogeny is also required to estimate gene gain and loss rates. To make the tree, we considered only gene families that were single copy in all six species and that had another arthropod species also represented with a single copy as an outgroup. Outgroup species were ranked based on the number of families in which they were also single copy along with the coleopteran species, and the highest ranking outgroup available was chosen for each family. For instance, *Pediculus humanus* was the most common outgroup species. For any gene family, we chose *P. humanus* as the outgroup if it was also single copy. If it was not, we chose the next highest ranking species as the outgroup for that family. This process resulted in 3,932 single copy orthologs that we subsequently aligned with PASTA¹⁶¹. We used RAxML¹⁶² with the PROTGAMMAJTTF model to make gene trees from the alignments and ASTRAL¹⁶³ to make the species tree. ASTRAL does not give branch lengths on its trees, a necessity for gene family analysis, so the species tree was again given to RAxML along with a concatenated alignment of all one-to-one orthologs for branch length estimation. Finally, to generate an ultrametric species tree with branch lengths in millions of years (my) we used the software r8s¹⁶⁴, with a calibration range based on age estimates of a crown Coleopteran fossil at 208.5–411 my¹⁶⁵. This calibration point itself was estimated in a similar fashion in a larger phylogenetic analysis of all 72 Arthropod species (Thomas *et al.* in prep).

With the gene family data and ultrametric phylogeny as the input data, gene gain and loss rates ($\lambda$) were estimated with CAFE v3.0¹⁶⁶. CAFE is able to estimate the amount of assembly and annotation error ($\varepsilon$) present in the input data using a distribution across the observed gene family counts and a pseudo-likelihood search, and then is able to correct for this error and obtain a more accurate estimate of $\lambda$. Our analysis had an $\varepsilon$ value of about 0.02, which implies that 3% of gene families have observed counts that are not equal to their true counts. After correcting for this error rate, $\lambda = 0.0010$ is on par with those previously found for other Arthropod orders (Thomas *et al.* in prep). Using the estimated $\lambda$ value, CAFE infers ancestral gene counts and calculates p-values across the tree for each family to assess the significance of any gene family changes along a given branch. Those branches with low p-values are considered rapidly evolving.

**Gene Expression Analysis.** RNAseq analyses were conducted to (1) establish male, female, and larva-enriched gene sets and (2) identify specific genes that are enriched within the digestive tract (mid-gut) compared to entire larva. RNAseq datasets were trimmed with CLC Genomics v.9 (Qiagen) and quality was assessed with FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Each dataset was mapped to the predicted gene set (OGS v1.0) using CLC Genomics. Reads were mapped with >90% similarity over 60% of length, with two mismatches allowed. The number of reads were corrected to reads per million mapped to allow comparison between RNAseq datasets that having varying coverage. Transcripts per million (TPM) was used as a proxy for gene expression and fold changes were determined as the TPM in one sample relative to the TPM of another dataset¹⁶⁷. The Baggerly's test (t-type test statistic) followed by Bonferroni correction was used to identify genes with significant enrichment in a specific sample¹⁶⁸. Statistical values for Bonferroni correction were reported as the number of genes $x \alpha$ value. This stringent statistical analysis was used as only a single replicate was available for each treatment. Enriched genes were removed, and mapping and expression analyses were repeated to ensure low expressed genes were not missed. Genes were identified by BLASTx searching against the NCBI non-redundant protein databases for arthropods with an expectation value (E-value) <0.001.

**Transposable Elements.** We investigated the identities (family membership), diversity, and genomic distribution of active transposable elements within *L. decemlineata* in order to understand their contribution to genome structure and to determine their potential positional effect on genes of interest (focusing on a TE neighborhood size of 1 kb). To identify TEs and analyze their distribution within the genome, we developed three repeat databases using: (1) RepeatMasker¹⁶⁹, which uses the library repeats within Repbase (http://www.girinst.org/repbase/), (2) the program RepeatModeler¹⁷⁰, which identifies de-novo repeat elements, and (3) literature searches to identify beetle transposons that were not found within Repbase. The three databases were used within RepeatMasker to determine the overall TE content in the genome.

To eliminate false positives and examine the genome neighborhood surrounding active TEs, all TE candidate models were translated in 6 frames and scanned for protein domains from the Pfam and CDD database (using the software transeq from Emboss, hmmer3 and rps-blast). The protein domain annotations were manually curated in order to remove: (a) clear false positives, (b) old highly degraded copies of TEs without identifiable coding potential, and (c) the correct annotation when improper labels were given. The TE models that contained protein
domains were mapped onto the genome and used for our neighborhood analysis: we extracted the 1 kb flanking regions for each gene and scanned these regions for TEs with intact protein coding domains.

**Population Genetic Variation and Demographic Analysis.** Population genetic diversity of pooled RNAseq samples was used to examine genetic structure of pest populations and past population demography. For Wisconsin, Michigan and the lab strains from New Jersey, we aligned the RNAseq data to the genomic scaffolds, using Bowtie2 version 2.1.0¹⁷¹ to index the genome and generate aligned SAM files. We used bwa¹⁵⁰ to align the RNAseq from the three populations from Europe. SAMtools/BCFtools version 0.1.19¹⁷² was used to produce BAM and VCF files. All calls were filtered with VCFtools version 0.1.11¹⁷³ using a minimum quality score of 30 and minimum depth of 10. All indels were removed from this study. Population specific VCF files were sorted and merged using VCFtools, and the allele counts were extracted for each SNP. These allele frequency data were then used to infer population splits and relative rates of genetic drift using Treemix version 1.12¹⁷⁴. We ran Treemix with SNPs in groups of 1000, choosing to root the tree with the Wisconsin population. In addition, for each pair of populations, we estimated the average genetic divergence by using F-statistics ($F_{ST}$) in VCFtools to calculate the ratio of among- to within-population genetic divergence across SNP loci.

To infer patterns of demographic change in the Midwestern USA (Wisconsin and Michigan) and European populations, the genome-wide allele frequency spectrum was used in *dadi* version 1.6.3¹⁷⁵ to infer demographic parameters under several alternative models of population history. The history of *L. decemlineata* as a pest is relatively well-documented. The introduction of *L. decemlineata* into Europe in 1914¹⁷⁶ is thought to have involved a strong bottleneck¹⁵ followed by rapid expansion. Similarly, an outbreak of *L. decemlineata* in Nebraska in 1859 is thought to have preceded population expansion into the Midwest reaching Wisconsin in 1865³. For each population, a constant-size model, a two-epoch model of instantaneous population size change at a time point $\tau$, a bottle-growth model of instantaneous size change followed by exponential growth, and a three-epoch model with a population size change of fixed duration followed by exponential growth, was fit to infer $\theta$, the product of the ancestral effective population size and mutation rate, and relative population size changes.

**Availability of data and material.** All data generated or analyzed during this study have been made publicly available (see Methods for NCBI accession numbers), or included in this published article and its supplementary information. The genome assembly and official gene sets can be accessed at: https://data.nal.usda.gov/dataset/leptinotarsa-decemlineata-genome-assembly-10_5667, https://data.nal.usda.gov/dataset/leptinotarsa-decemlineata-genome-annotations-v053_5668 and https://data.nal.usda.gov/dataset/leptinotarsa-decemlineata-official-gene-set-v11.

**References**

1. Grafius, E. Economic impact of insecticide resistance in the Colorado potato beetle (Coleoptera: Chrysomelidae) on the Michigan potato industry. *J. Econ. Entomol.* **90**, 1144–1151 (1997).

2. Skryabin, K. Do Russia and Eastern Europe need GM plants? *Nat. Biotechnol.* **27**, 593–595 (2010).

3. Walsh, B. D. The new potato bug and its natural history. *Practical Entomol.* **1**, 1–4 (1865).

4. Gauthier, N. L., Hofmaster, R. N. & Semel, M. History of Colorado potato beetle control in *Advances in potato pest management* (ed. Lashomb, J. H. & Casagrande, R.) 13–33 (Hutchinson Ross Publishing Co., 1981).

5. Hare, J. D. Ecology and management of the Colorado potato beetle. *Annu. Rev. Entomol.* **41**, 81–100 (1990).

6. Weber, D. C. Colorado potato beetle, *Leptinotarsa decemlineata* (Say) (Coleoptera: Chrysomelidae) in *Encyclopedia of entomology* (ed. Capinera, J. L.) 1008–1013 (Springer Netherlands, 2008).

7. Lyytinen, A., Boman, S., Grapputo, A., Lindström, L. & Mappes, J. Cold tolerance during larval development: effects on the thermal distribution limits of *Leptinotarsa decemlineata*. *Entomol. Exp. Appl.* **133**, 92–99 (2009).

8. Hiiesaar, K. et al. Factors affecting development and overwintering of second generation Colorado potato beetle (Coleoptera: Chrysomelidae) in Estonia in 2010. *Acta Agric. Scand. Sect. B - Soil Plant Sci.* **63**, 506–515 (2013).

9. Hiiesaar, K. et al. Phenology and overwintering of the Colorado potato beetle *Leptinotarsa decemlineata* Say in 2008–2015 in Estonia. *Acta Agric. Scand. Sect. B - Soil Plant Sci.* **66**, 1–8 (2016).

10. Alyokhin, A., Baker, M., Mota-Sanchez, D., Dively, G. & Grafius, E. Colorado potato beetle resistance to insecticides. *Am J Potato Res.* **85**, 395–413 (2008).

11. Sparks, T. C. & Nauen, R. IRAC: Mode of action classification and insecticide resistance management. *Pestic. Biochem. Physiol.* **121**, 122–128 (2015).

12. Izzo, V., Chen, Y. H., Schoville, S. D., Wang, C. & Hawthorne, D. J. Origin of pest lineages of the Colorado potato beetle, *Leptinotarsa decemlineata*. *bioRxiv* 156612, https://doi.org/10.1101/156612 (2017).

13. Lehmann, P., Lyytinen, A., Piroinen, S. & Lindström, L. Northward range expansion requires synchronization of both overwintering behaviour and physiology with photoperiod in the invasive Colorado potato beetle (*Leptinotarsa decemlineata*). *Oecologia* **176**, 57–68 (2014).

14. Liu, N., Li, Y. & Zhang, R. Invasion of Colorado potato beetle, *Leptinotarsa decemlineata*, in China: dispersal, occurrence, and economic impact. *Entomol. Exp. Appl.* **143**, 207–217 (2012).

15. Grapputo, A., Boman, S., Lindström, L., Lyytinen, A. & Mappes, J. The voyage of an invasive species across continents: genetic diversity of North American and European Colorado potato beetle populations. *Mol. Ecol.* **14**, 4207–4219 (2005).

16. Azeredo-Espin, A., Schroder, R., Huettel, M. & Sheppard, W. Mitochondrial DNA variation in geographic populations of Colorado potato beetle, *Leptinotarsa decemlineata* (Coleoptera: Chrysomelidae). *Experientia* **47**, 483–485 (1991).

17. Lu, W., Kennedy, G. G. & Gould, F. Genetic analysis of larval survival and larval growth of two populations of *Leptinotarsa decemlineata* on tomato. *Entomol. Exp. Appl.* **99**, 143–155 (2001).

18. Gaston, K. J. The how and why of biodiversity. *Nature* **421**, 900–901 (2003).

19. Turnock, W. J. & Fields, P. G. Winter climates and cold hardiness in terrestrial insects. *Eur. J. Entomol.* **102**, 561–576 (2005).

20. Piroinen, S., Ketola, T., Lyytinen, A. & Lindström, L. Energy use, diapause behaviour and northern range expansion potential in the invasive Colorado potato beetle. *Funct. Ecol.* **25**, 527–536 (2011).

21. Lehmann, P., Lyytinen, A., Piroinen, S. & Lindström, L. Latitudinal differences in diapause related photoperiodic responses of European Colorado potato beetles (*Leptinotarsa decemlineata*). *Evol. Ecol.* **29**, 269–282 (2015).

22. Hsiao, T. H. Host plant adaptations among geographic populations of the Colorado potato beetle. *Entomol. Exp. Appl.* **24**, 237–247 (1978).
23. Jolivet, P. Food habits and food selection of Chrysomelidae. Bionomic and evolutionary perspectives in *Biology of the chrysomelidae* (ed. Jolivet, P., Petitpierre, E. & Hsiao, T. H.) 1-24 (Kluwer Acad. Press, 1988).

24. Hsiao, T. H. Host specificity, seasonality and bionomics of *Leptinotarsa* beetles in *Biology of the chrysomelidae* (ed. Jolivet, P., Petitpierre, E. & Hsiao, T. H.) 581-599 (Kluwer Acad. Press, 1988).

25. Hare, J. D. & Kennedy, G. G. Genetic variation in plant-insect associations: survival of *Leptinotarsa decemlineata* populations on *Solanum carolinense*. *Evolution* **40**, 1031-1043 (1986).

26. Jacques, R. L. Jr. & Fasulo, T. R. Colorado potato beetle, *Leptinotarsa decemlineata* (Say) and False potato beetle, *Leptinotarsa juncta* (Germar) (Insecta: Coleoptera: Chrysomelidae). *Univ. Florida IFAS Ext.* EENY146, 1-5 (2000).

27. Alyokhin, A. & Chen, Y. H. Adaptation to toxic hosts as a factor in the evolution of insecticide resistance. *Curr. Opin. Insect Sci.* **2**, 33-38 (2017).

28. Cárdenas, P. D. et al. The bitter side of the nightshades: Genomics drives discovery in Solanaceae steroidal alkaloid metabolism. *Phytochemistry* **113**, 24-32 (2015).

29. Milner, S. E. et al. Bioactivities of glycoalkaloids and their aglycones from *Solanum* species. *J. Agric. Food Chem.* **59**, 3454-3484 (2011).

30. Dimock, M. B. & Tingey, W. M. Host acceptance behaviour of Colorado potato beetle larvae influenced by potato glandular trichomes. *Physiol. Entomol.* **13**, 399-406 (1988).

31. Lawrence, S. D., Novak, N. G., Ju, C. J.-T. & Cooke, J. E. K. Examining the molecular interaction between potato (*Solanum tuberosum*) and Colorado potato beetle *Leptinotarsa decemlineata*. *Botany* **86**, 1080-1091 (2008).

32. Petek, M. et al. A complex of genes involved in adaptation of *Leptinotarsa decemlineata* larvae to induced potato defense. *Arch. Insect Biochem. Physiol.* **79**, 153-181 (2012).

33. Novillo, C., Castañera, P. & Ortego, F. Characterization and distribution of chymotrypsin-like and other digestive proteases in Colorado potato beetle larvae. *Arch. Insect Biochem. Physiol.* **36**, 181-201 (1997).

34. Armer, C. A. Colorado potato beetle toxins revisited: Evidence the beetle does not sequester host plant glycoalkaloids. *J. Chem. Ecol.* **30**, 883-888 (2004).

35. Deroe, C. & Pasteels, J. M. Defensive mechanisms against predation in the Colorado beetle (*Leptinotarsa decemlineata*, Say). *Arch. Biol. Sci.* **88**, 289-304 (1977).

36. Hsiao, T. H. & Fraenkel, G. Properties of Leptinotarsin: a toxic hemolymph protein from the Colorado potato beetle. *Toxicon* **7**, 119-130 (1969).

37. Forgash, A. J. Insecticide resistance in the Colorado potato beetle in *Proceedings of the symposium on the Colorado potato beetle*, XVII International Congress of Entomology, August 1984 (ed. Ferro, D. N. & Voss, R. H.). 33-52 (Massachusetts Agricultural Experimental Station, 1985).

38. Wan, P.-J. *et al*. Identification of cytochrome P450 monooxygenase genes and their expression profiles in cyhalothrin-treated Colorado potato beetle, *Leptinotarsa decemlineata*. *Pestic. Biochem. Physiol.* **107**, 360-368 (2013).

39. Zhang, J., Pelletier, Y. & Goyer, C. Identification of potential detoxification enzyme genes in *Leptinotarsa decemlineata* (Say) and study of their expression in insects reared on different plants. *J. Plant Sci.* **88**, 621-629 (2008).

40. Lü, F.-G. et al. Identification of carboxylesterase genes and their expression profiles in the Colorado potato beetle *Leptinotarsa decemlineata* treated with fipronil and cyhalothrin. *Pestic. Biochem. Physiol.* **122**, 86-95 (2015).

41. Zhu, F., Xu, J., Palli, R., Ferguson, J. & Palli, S. R. Ingested RNA interference for managing the populations of the Colorado potato beetle. *Leptinotarsa decemlineata*. *Pest Manag. Sci.* **67**, 175-182 (2011).

42. Clements, J., Schoville, S., Peterson, N., Lan, Q. & Groves, R. L. Characterizing molecular mechanisms of imidacloprid resistance in select populations of *Leptinotarsa decemlineata* in the Central Sands region of Wisconsin. *PLoS One* **11**, e0147844 (2016).

43. Kumar, A. et al. Sequencing, de novo assembly and annotation of the Colorado potato beetle, *Leptinotarsa decemlineata*, transcriptome. *PLoS One* **9**, e86012 (2014).

44. Rinkevich, F. D. et al. Multiple evolutionary origins of knockdown resistance (kdr) in pyrethroid-resistant Colorado potato beetle. *Leptinotarsa decemlineata*. *Pestic. Biochem. Physiol.* **104**, 192-200 (2012).

45. Zhou, L. T. et al. RNA interference of a putative S-adenosyl-L-homocysteine hydrolase gene affects larval performance in *Leptinotarsa decemlineata* (Say). *J. Insect Physiol.* **59**, 1049-1056 (2013).

46. Zhang, J. et al. Full crop protection from an insect pest by expression of long double-stranded RNAs in plastids. *Science* **347**, 991-994 (2015).

47. Mota-Sanchez, D., Hollingworth, R. M., Grafius, E. J. & Moyer, D. D. Resistance and cross-resistance to neonicotinoid insecticides and spinosad in the Colorado potato beetle, *Leptinotarsa decemlineata* (Say) (Coleoptera: Chrysomelidae). *Pest Manag. Sci.* **62**, 30-37 (2006).

48. Zhu, F., Moural, T. W., Nelson, D. R. & Palli, S. R. A specialist herbivore pest adaptation to xenobiotics through up-regulation of multiple cytochrome P450s. *Sci. Rep.* **6**, 20421 (2016).

49. Butler, J. et al. ALLPATHS: de novo assembly of whole-genome shotgun microreads. *Genome Res.* **18**, 810-20 (2008).

50. Gregory, T. R. Animal genome size database. http://www.genomesize.com (2017).

51. Tribolium Genome Sequencing Consortium. The genome of the model beetle and pest *Tribolium castaneum*. *Nature* **452**, 949-55 (2008).

52. Keeling, C. I. et al. Draft genome of the mountain pine beetle, *Dendroctonus ponderosae* Hopkins, a major forest pest. *Genome Biol.* **14**, R27 (2013).

53. Cunningham, C. B. et al. The genome and methylome of a beetle with complex social behavior, *Nicrophorus vespilloides* (Coleoptera: Silphidae). *Genome Biol. Evol.* **7**, 3383-3396 (2015).

54. McKenna, D. D. et al. Genome of the Asian longhorned beetle (*Anoplophora glabripennis*), a globally significant invasive species, reveals key functional and evolutionary innovations at the beetle-plant interface. *Genome Biol.* **17**, 227 (2016).

55. Vega, F. E. et al. Draft genome of the most devastating insect pest of coffee worldwide: the coffee berry borer, *Hypothenemus hampei*. *Sci. Rep.* **5**, 12525 (2015).

56. Petitpierre, E., Segarra, C. & Juan, C. Genome size and chromosomal evolution in leaf beetles (Coleoptera, Chrysomelidae). *Hereditas* **6**, 1-6 (1993).

57. Pryszcz, L. P. & Gabaldón, T. Redundans: an assembly pipeline for highly heterozygous genomes. *Nucleic Acids Res.* **44**, e113-e113 (2016).

58. Denton, J. F. et al. Extensive error in the number of genes inferred from draft genome assemblies. *PLoS Comput. Biol.* **10**, e1003998 (2014).

59. Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. *Bioinformatics* **31**, 3210-3212 (2015).

60. Gruden, K. et al. The cysteine protease activity of Colorado potato beetle (*Leptinotarsa decemlineata* Say) guts, which is insensitive to potato protease inhibitors, is inhibited by thyroglobulin type-1 domain inhibitors. *Insect Biochem. Mol. Biol.* **28**, 549-560 (1998).

61. Smilanich, A. M., Dyer, L. A., Chambers, J. Q. & Bowers, M. D. Immunological cost of chemical defence and the evolution of herbivore diet breadth. *Ecol. Lett.* **12**, 612-621 (2009).

62. Vidal, N. M., Grazziotin, A. L., Iyer, L. M., Aravind, L. & Venancio, T. M. Transcription factors, chromatin proteins and the diversification of Hemiptera. *Insect Biochem. Mol. Biol.* **69**, 1-13 (2016).

63. Najafabadi, H. S. et al. C2H2 zinc finger proteins greatly expand the human regulatory lexicon. *Nat. Biotech.* **33**, 555-562 (2015).
64. Chénais, B., Caruso, A., Hiard, S. & Casse, N. The impact of transposable elements on eukaryotic genomes: from genome size increase to genetic adaptation to stressful environments. *Gene* **509**, 7–15 (2012).

65. Tu, Z. Insect transposable elements in *Insect molecular biology and biochemistry* (ed. Gilbert, L. I.) 57–89 (Academic Press, 2012).

66. Biémont, C. & Vieira, C. Genetics: junk DNA as an evolutionary force. *Nature* **443**, 521–524 (2006).

67. González, J. & Petrov, D. A. The adaptive role of transposable elements in the *Drosophila* genome. *Gene* **448**, 124–133 (2009).

68. Bourque, G. et al. Evolution of the mammalian transcription factor binding repertoire via transposable elements. *Genome Res.* **18**, 1752–1762 (2008).

69. Cordaux, R. & Batzer, M. A. The impact of retrotransposons on human genome evolution. *Nat. Rev. Genet.* **10**, 691–703 (2009).

70. van de Lagemaat, L. N., Landry, J. R., Mager, D. L. & Medstrand, P. Transposable elements in mammals promote regulatory variation and diversification of genes with specialized functions. *Trends Genet.* **19**, 530–536 (2003).

71. Lavoie, C. A., Platt, R. N., Novick, P. A., Counterman, B. A. & Ray, D. A. Transposable element evolution in *Heliconius* suggests genome diversity within Lepidoptera. *Mob DNA* **4**, 21 (2013).

72. Osanai-Futahashi, M., Suetsugu, Y., Mita, K. & Fujiwara, H. Genome-wide screening and characterization of transposable elements and their distribution analysis in the silkworm, *Bombyx mori*. *Insect Biochem. Mol. Biol.* **38**, 1046–1057 (2008).

73. Schrader, L. et al. Transposable element islands facilitate adaptation to novel environments in an invasive species. *Nat. Commun.* **5**, 5495 (2014).

74. González, J., Karasov, T. L., Messer, P. W. & Petrov, D. A. Genome-wide patterns of adaptation to temperate environments associated with transposable elements in *Drosophila*. *PLoS Genet.* **6**, e1000905 (2010).

75. Cridland, J. M., Thornton, K. R. & Long, A. D. Gene expression variation in *Drosophila melanogaster* due to rare transposable element insertion alleles of large effect. *Genetics* **199**, 85–93 (2015).

76. Daborn, P. J. et al. A single P450 allele associated with insecticide resistance in *Drosophila*. *Science* **297**, 2253–2256 (2002).

77. Messer, P. W. & Petrov, D. A. Population genomics of rapid adaptation by soft selective sweeps. *Trends Ecol. Evol.* **28**, 659–669 (2013).

78. Aquadro, C. F., DuMont, V. B. & Reed, F. A. Genome-wide variation in the human and fruitfly: a comparison. *Curr. Opin. Genet. Dev.* **11**, 627–634 (2001).

79. International Chicken Polymorphism Map Consortium. A genetic variation map for chicken with 2.8 million single-nucleotide polymorphisms. *Nature* **432**, 717–722 (2004).

80. Choi, J.-H. et al. Gene discovery in the horned beetle *Onthophagus taurus*. *BMC Genomics* **11**, 703 (2010).

81. Morlais, I., Ponçon, N., Simard, F., Cohuet, A. & Fontenille, D. Intraspecific nucleotide variation in *Anopheles gambiae*: new insights into the biology of malaria vectors. *Am. J. Trop. Med. Hyg.* **71**, 795–802 (2004).

82. Campos, J. L., Zhao, L. & Charlesworth, B. Estimating the parameters of background selection and selective sweeps in *Drosophila* in the presence of gene conversion. *Proc. Natl Acad. Sci. USA* **114**, E4762–E4771 (2017).

83. Dangles, O., Irschick, D., Chittka, L. & Casas, J. Variability in sensory ecology: expanding the bridge between physiology and evolutionary biology. *Q. Rev. Biol.* **84**, 51–74 (2009).

84. Matsuura, H., Sokabe, T., Kohno, K., Tominaga, M. & Kadowaki, T. Evolutionary conservation and changes in insect TRP channels. *BMC Evol Biol.* **9**, 228 (2009).

85. Venkatachalam, K. & Montell, C. TRP channels. *Annu. Rev. Biochem.* **76**, 387–417 (2007).

86. Nesterov, A. et al. TRP channels in insect stretch receptors as insecticide targets. *Neuron* **86**, 665–671 (2015).

87. Hauser, F. et al. A genome-wide inventory of neurohormone GPCRs in the red flour beetle *Tribolium castaneum*. *Front. Neuroendocrinol.* **29**, 142–165 (2008).

88. Nishi, Y., Sasaki, K. & Miyatake, T. Biogenic amines, caffeine and tonic immobility in *Tribolium castaneum*. *J. Insect Physiol.* **56**, 622–628 (2010).

89. Roeder, T. Biochemistry and molecular biology of receptors for biogenic amines in locusts. *Microsc. Res. Tech.* **56**, 237–247 (2002).

90. Schoonhoven, L., Van Loon, J. J. A. & Dicke, M. *Insect-Plant Biology* 2nd edn (Oxford University Press, 2005).

91. Leal, W. S. Odorant reception in insects: roles of receptors, binding proteins, and degrading enzymes. *Annu. Rev. Entomol.* **58**, 373–391 (2013).

92. Pelosi, P., Zhou, J.-J., Ban, L. P. & Calvello, M. Soluble proteins in insect chemical communication. *Cell. Mol. Life Sci.* **63**, 1658–1676 (2006).

93. Robertson, H. M., Warr, C. G. & Carlson, J. R. Molecular evolution of the insect chemoreceptor gene superfamily in *Drosophila melanogaster*. *Proc. Natl. Acad. Sci. USA* **100**, 14537–14542 (2003).

94. Sato, K., Tanaka, K. & Touhara, K. Sugar-regulated cation channel formed by an insect gustatory receptor. *Proc. Natl. Acad. Sci. USA* **108**, 11680–11685 (2011).

95. Rytz, R., Croset, V. & Benton, R. Ionotropic Receptors (IRs): Chemosensory ionotropic glutamate receptors in *Drosophila* and beyond. *Insect Biochem. Mol. Biol.* **43**, 888–897 (2013).

96. Jackowska, M. et al. Genomic and gene regulatory signatures of cryptozoic adaptation: Loss of blue sensitive photoreceptors through expansion of long wavelength-opsin expression in the red flour beetle *Tribolium castaneum*. *Front. Zool.* **4**, 24 (2007).

97. Sharkey, C. R. et al. Overcoming the loss of blue sensitivity through opsin duplication in the largest animal group, beetles. *Sci. Rep.* **7**, 8 (2017).

98. Izzo, V. M., Mercer, N., Armstrong, J. & Chen, Y. H. Variation in host usage among geographic populations of *Leptinotarsa decemlineata*, the Colorado potato beetle. *J. Pest. Sci.* **87**, 597–608 (2014).

99. Otálora-Luna, F. & Dickens, J. C. Multimodal stimulation of Colorado potato beetle reveals modulation of pheromone response by yellow light. *PLoS One* **6**, e20990 (2011).

100. Briscoe, A. D. & Chittka, L. The evolution of color vision in insects. *Annu. Rev. Entomol.* **46**, 471–510 (2001).

101. Leung, N. Y. & Montell, C. Unconventional roles of opsins. *Annu. Rev. Cell. Dev. Biol.* **33**, 241–264 (2017).

102. Price, P. W., Denno, R. F., Eubanks, M. D., Finke, D. L. & Kaplan, I. *Insect ecology: behavior, populations and communities* (Cambridge University Press, 2011).

103. Martinez, M. et al. Phytocystatins: defense proteins against phytophagous insects and Acari. *Int. J. Mol. Sci.* **17**, 1747 (2016).

104. Wolfson, J. L. & Murdock, L. L. Suppression of larval Colorado potato beetle growth and development by digestive proteinase inhibitors. *Entomol. Exp. Appl.* **44**, 235–240 (1987).

105. Bolter, C. J. & Jongsma, M. A. Colorado potato beetles (*Leptinotarsa decemlineata*) adapt to proteinase inhibitors induced in potato leaves by methyl jasmonate. *J. Insect Physiol.* **41**, 1071–1078 (1995).

106. Murdock, L. L. et al. Cysteine digestive proteinases in Coleoptera. *Biochem. Physiol.* **87**, 783–787 (1987).

107. Chye, M.-L., Sin, S.-F., Xu, Z.-F. & Yeung, E. C. Serine proteinase inhibitor proteins: Exogenous and endogenous functions. *Vitr. Cell. Dev. Biol.-Plant* **42**, 100–108 (2006).

108. Gruden, K., Popovič, T., Cimerman, N., Križaj, I. & Štrukelj, B. Diverse enzymatic specificities of digestive proteases, “intestains”, enable Colorado potato beetle larvae to counteract the potato defence mechanism. *Biol. Chem.* **384**, 305–310 (2003).

109. Gruden, K. et al. Molecular basis of Colorado potato beetle adaptation to potato plant defence at the level of digestive cysteine proteinases. *Insect Biochem. Mol. Biol.* **34**, 365–375 (2004).

110. Goptar, I. A. et al. Cysteine digestive peptidases function as post-glutamine cleaving enzymes in tenebrionid stored-product pests. *Comp. Biochem. Physiol. Part B* **161**, 148–154 (2012).

111. Turk, V. et al. Cysteine cathepsins: From structure, function and regulation to new frontiers. *Biochim. Biophys. Acta - Proteins Proteomics* **1824**, 68–88 (2012).
112. Marchler-Bauer, A. et al. CDD: NCBI's conserved domain database. *Nucleic Acids Res.* **43**, D222–D226 (2014).

113. Finn, R. D. et al. Pfam: the protein families database. *Nucleic Acids Res.* **42**, D222–D230 (2013).

114. Wex, T. et al. TIN-ag-RP, a novel catalytically inactive cathepsin B-related protein with EGF domains, is predominantly expressed in vascular smooth muscle cells. *Biochemistry* **40**, 1350–1357 (2001).

115. Yamamoto, Y., Kurata, M., Watabe, S., Murakami, R. & Takahashi, S. Y. Novel cysteine proteinase inhibitors homologous to the proregions of cysteine proteinases. *Curr. Protein Pept. Sci.* **3**, 231–238 (2002).

116. Sainsbury, F., Rhéaume, A.-J., Goulet, M.-C., Vorster, J. & Michaud, D. Discrimination of differentially inhibited cysteine proteases by activity-based profiling using cystatin variants with tailored specificities. *J. Proteome Res.* **11**, 5983–5993 (2012).

117. Martynov, A. G., Elpidina, E. N., Perkin, L. & Oppert, B. Functional analysis of C1 family cysteine peptidases in the larval gut of *Tenebrio molitor* and *Tribolium castaneum*. *BMC Genomics* **16**, 75 (2015).

118. Levasseur, A., Drula, E., Lombard, V., Coutinho, P. M. & Henriissat, B. Expansion of the enzymatic repertoire of the CAZy database to integrate auxiliary redox enzymes. *Biotechnol. Biofuels* **6**, 41 (2013).

119. Scully, E. D., Hoover, K., Carlson, J. E., Tien, M. & Geib, S. M. Mid-gut transcriptome profiling of *Anoplophora glabripennis*, a lignocellulose degrading cerambycid beetle. *BMC Genomics* **14**, 850 (2013).

120. Kirsch, R. et al. Horizontal gene transfer and functional diversification of plant cell wall degrading polygalacturonases: key events in the evolution of herbivory in beetles. *Insect Biochem. Mol. Biol.* **52**, 33–50 (2014).

121. Pauchet, Y., Wilkinson, P., Chauhan, R. & Ffrench-Constant, R. H. Diversity of beetle genes encoding novel plant cell wall degrading enzymes. *PLoS One* **5**, e15635 (2010).

122. Meyer, U. A. Overview of enzymes of drug metabolism. *J. Pharmacokinet. Pharmacodyn.* **24**, 449–459 (1996).

123. Bloomquist, J. R. Chloride channels as tools for developing selective insecticides. *Arch. Insect Biochem. Physiol.* **54**, 145–156 (2003).

124. Le Goff, G., Hamon, A., Bergé, J.-B. & Amichot, M. Resistance to fipronil in *Drosophila simulans*: influence of two point mutations in the RDL GABA receptor subunit. *J. Neurochem.* **92**, 1295–1305 (2005).

125. Jones, A. K. & Sattelle, D. B. The cys-loop ligand-gated ion channel gene superfamily of the nematode, *Caenorhabditis elegans*. *Invert. Neurosci.* **8**, 41–47 (2008).

126. Jones, A. K. & Sattelle, D. B. Diversity of insect nicotinic acetylcholine receptor subunits in *Insect nicotinic acetylcholine receptors*. (ed. Thany, S. H.) 25–43 (Springer, 2010).

127. Ffrench-Constant, R. H., Mortlock, D. P., Shaffer, C. D., MacIntyre, R. J. & Roush, R. T. Molecular cloning and transformation of cyclodiene resistance in *Drosophila*: an invertebrate gamma-aminobutyric acid subtype A receptor locus. *Proc. Natl. Acad. Sci. USA* **88**, 7209–7213 (1991).

128. Ffrench-Constant, R. H. & Rocheleau, T. A. *Drosophila* γ-aminobutyric acid receptor gene Rdl shows extensive alternative splicing. *J. Neurochem.* **60**, 2323–2326 (1993).

129. Buckingham, S. D., Biggin, P. C., Sattelle, B. M., Brown, L. A. & Sattelle, D. B. Insect GABA receptors: splicing, editing, and targeting by antiparasitics and insecticides. *Mol. Pharmacol.* **68**, 942–951 (2005).

130. Clements, J., Schoville, S., Peterson, N., Lan, Q. & Groves, R. L. Characterizing molecular mechanisms of Imidacloprid resistance in select populations of *Leptinotarsa decemlineata* in the Central Sands Region of Wisconsin. *PLoS One* **11**, e0147844 (2016).

131. Clements, J. et al. RNA interference of three up-regulated transcripts associated with insecticide resistance in an imidacloprid resistant population of *Leptinotarsa decemlineata*. *Pestic. Biochem. Physiol.* **135**, 35–40 (2017).

132. Willis, J. H. Structural cuticular proteins from arthropods: Annotation, nomenclature, and sequence characteristics in the genomics era. *Insect Biochem. Mol. Biol.* **40**, 189–204 (2010).

133. McDonnell, C. M. et al. Evolutionary toxicogenomics: diversification of the Cyp12d1 and Cyp12d3 Genes in *Drosophila* species. *J. Mol. Evol.* **74**, 281–96 (2012).

134. Hu, F., Dou, W., Wang, J.-J., Jia, F.-X. & Wang, J.-J. Multiple glutathione S-transferase genes: identification and expression in oriental fruit fly, *Bactrocera dorsalis*. *Pest Manag. Sci.* **70**, 295–303 (2014).

135. Han, J.-B., Li, G.-Q., Wan, P.-J., Zhu, T.-T. & Meng, Q.-W. Identification of glutathione S-transferase genes in *Leptinotarsa decemlineata* and their expression patterns under stress of three insecticides. *Pestic. Biochem. Physiol.* **133**, 26–34 (2016).

136. Fire, A. et al. Potent and specific genetic interference by double-stranded RNA in *Caenorhabditis elegans*. *Nature* **391**, 806–811 (1998).

137. Siomi, H. & Siomi, M. C. On the road to reading the RNA-interference code. *Nature* **457**, 396–404 (2009).

138. Revuelta, L. et al. Contribution of Idace1 gene to acetylcholinesterase activity in Colorado potato beetle. *Insect Biochem. Mol. Biol.* **41**, 795–803 (2011).

139. Wan, P.-J., Fu, K.-Y., Lü, F.-G., Guo, W.-C. & Li, G.-Q. Knockdown of a putative alanine aminotransferase gene affects amino acid content and flight capacity in the Colorado potato beetle *Leptinotarsa decemlineata*. *Amino Acids* **47**, 1445–1454 (2015).

140. Baum, J. A. et al. Control of coleopteran insect pests through RNA interference. *Nat. Biotechnol.* **25**, 1322 (2007).

141. Swevers, L. et al. Colorado potato beetle (Coleoptera) gut transcriptome analysis: expression of RNA interference-related genes. *Insect Mol. Biol.* **22**, 668–684 (2013).

142. Bernstein, E., Caudy, A. A., Hammond, S. M. & Hannon, G. J. Role for a bidentate ribonuclease in the initiation step of RNA interference. *Nature* **409**, 363–366 (2001).

143. Grishok, A. et al. Genes and mechanisms related to RNA interference regulate expression of the small temporal RNAs that control *C. elegans* developmental timing. *Cell* **106**, 23–34 (2001).

144. Yoon, J.-S., Shukla, J. N., Gong, Z. J., Mogilicherla, K. & Palli, S. R. RNA interference in the Colorado potato beetle, *Leptinotarsa decemlineata*: Identification of key contributors. *Insect Biochem. Mol. Biol.* **78**, 78–88 (2016).

145. Guénin, H.-A. & Scherler, M. La formule chromosomale du doryphore *Leptinotarsa decemlineata* Stål. *Rev. Suisse Zool.* **58**, 359–370 (1951).

146. Hsiao, T. H. & Hsiao, C. Chromosomal analysis of *Leptinotarsa* and *Labidomera* species. *Genetica* **60**, 139–150 (1983).

147. Poelchau, M. et al. The i5k Workspace@NAL—enabling genomic data access, visualization and curation of arthropod genomes. *Nucleic Acids Res.* **43**, D714–D719 (2015).

148. Marcais, G. & Kingsford, C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. *Bioinformatics* **27**, 764–770 (2011).

149. Holt, C. & Yandell, M. MAKER2: an annotation pipeline and genome-database management tool for second-generation genome projects. *BMC Bioinformatics* **12**, 491 (2011).

150. Lee, E. et al. Web Apollo: a web-based genomic annotation editing platform. *Genome Biol.* **14**, R93 (2013).

151. Adams, M. D. et al. The genome sequence of *Drosophila melanogaster*. *Science* **287**, 2185–2195 (2000).

152. Camacho, C. et al. BLAST+: architecture and applications. *BMC Bioinformatics* **10**, 421 (2009).

153. Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. *J. Mol. Biol.* **215**, 403–410 (1990).

154. Larkin, M. A. et al. Clustal W and Clustal X version 2.0. *Bioinformatics* **23**, 2947–2948 (2007).

155. Laetsch, D. R. & Blaxter, M. L. BlobTools: Interrogation of genome assemblies. *F1000Research* **6**, 1287 (2017).

156. National Center for Biotechnology Information Nucleotide. Nucleotide (nt) database. Available from: https://www.ncbi.nlm.nih.gov/nucleotide/(accessed January 05, 2018).

157. Buchfink, B., Xie, C. & Huson, D. Fast and sensitive protein alignment using DIAMOND. *Nature Methods* **12**, 59–60 (2015).

158. The UniProt Consortium. UniProt: the universal protein knowledgebase. *Nucleic Acids Res.* **45**, D158–D169 (2017).

159. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**, 1754–1760 (2009).
160. Kriventseva, E. V. et al. OrthoDB v8: update of the hierarchical catalog of orthologs and the underlying free software. *Nucleic Acids Res.* **43**, D250–D256 (2015).

161. Mirarab, S. et al. PASTA: ultra-large multiple sequence alignment for nucleotide and amino-acid sequences. *J. Comput. Biol.* **22**, 377–386 (2015).

162. Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. *Bioinformatics* **30**, 2010–2011 (2014).

163. Mirarab, S. & Warnow, T. ASTRAL-II: coalescent-based species tree estimation with many hundreds of taxa and thousands of genes. *Bioinformatics* **31**, i44–i52 (2015).

164. Sanderson, M. J. r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock. *Bioinformatics* **19**, 301–302 (2003).

165. Wolfe, J. M., Daley, A. C., Legg, D. A. & Edgecombe, G. D. Fossil calibrations for the arthropod Tree of Life. *Earth-Sci. Rev.* **160**, 43–110 (2016).

166. Han, M. V., Thomas, G. W. C., Lugo-Martinez, J. & Hahn, M. W. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. *Mol. Biol. Evol.* **30**, 1987–1997 (2013).

167. Rosendale, A. J., Romick-Rosendale, L. E., Watanabe, M., Dunlevy, M. E. & Benoit, J. B. Mechanistic underpinnings of dehydration stress in the American dog tick revealed through RNA-seq and metabolomics. *J. Exp. Biol.* **219**, 1808–1819 (2016).

168. Baggerly, K. A., Deng, L., Morris, J. S. & Aldaz, C. M. Differential expression in SAGE: accounting for normal between-library variation. *Bioinformatics* **19**, 1477–1483 (2003).

169. Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0.

170. Smit, A. & Hubley, R. RepeatModeler Open-1.0.

171. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. *Nat. Methods* **9**, 357 (2012).

172. Li, H. et al. The sequence alignment/map format and SAMtools. *Bioinformatics* **25**, 2078–2079 (2009).

173. Danecek, P. et al. The variant call format and VCFtools. *Bioinformatics* **27**, 2156–2158 (2011).

174. Pickrell, J. K. & Pritchard, J. K. Inference of population splits and mixtures from genome-wide allele frequency data. *PLoS Genet.* **8**, e1002967 (2012).

175. Gutenkunst, R. N., Hernandez, R. D., Williamson, S. H. & Bustamante, C. D. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. *PLoS Genet.* **5**, e1000695 (2009).

176. Jacques, R. L. Jr. *The potato beetles* (No. 3) (Taylor & Francis, 1988).

## Acknowledgements

We sincerely thank the sequencing, assembly and annotation teams at the Baylor College of Medicine Human Genome Sequencing Center for their efforts. We would like to acknowledge the following funding sources: sequencing, assembly and automated annotation was supported by NIH grant NHGRI U54 HG003273 to RAG; the UVM Agricultural Experiment Station Hatch grant to YHC (VT-H02010); the NIH postdoctoral training grant to RFM (K12 GM000708); MMT's work with Apollo was supported by NIH grants (5R01GM080203 from NIGMS, and 5R01HG004483 from NHGRI) and by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy (contract No. DE-AC02-05CH11231); the National Science Centre (2012/07/D/NZ2/04286) and Ministry of Science and Higher Education scholarship to AM. Mention of trade names or commercial products in this publication is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the U.S. Department of Agriculture. USDA is an equal opportunity provider and employer.

## Author Contributions

M.N.A., A.B., J.H.B., K.C., A.K.C., O.C., E.M.D., E.N.E., P.E., M.F., I.G.-R., C.G., A.G., M.G., B.H., E.C.J., J.W.J., M.K., S.A.K., A.K., F.L., V.L., X.M., A.M., N.J.M., R.F.M., B.O., S.R.P., K.A.P., Y.P., L.C.P., M.P., É.R., J.P.R., H.M.R., A.J.R., V.M.R.-A., G.S., A.S.T., I.M.V.J., A.D.Y., G.D.Y., and J.-S.Y. contributed to the manual annotation effort, data analysis, and data interpretation, in addition to reading and approving the final manuscript. S.D.S. and Y.H.C. coordinated the project and drafted the manuscript. S.R. and R.A.G. coordinated genome sequencing, assembly and automated annotation at the Baylor College of Medicine Human Genome Sequencing Center. M.F.P., C.C., M.M.T., and M.-J.M.C. coordinated the biocuration of the genome. G.W.C.T. generated the phylogeny and conducted the gene family analysis. M.T.W. conducted the transcription factor analysis. A.G., K.G., M.P., and S.Z. contributed additional RNAseq data. J.B.B. coordinated the RNAseq data analysis. K.B., A.M., and Y.H.C. conducted the transposable element analysis. S.D.S. and J.C. conducted the population genetics analysis.

## Additional Information

Supplementary information accompanies this paper at https://doi.org/10.1038/s41598-018-20154-1.

**Competing Interests:** The authors declare that they have no competing interests.

**Publisher's note:** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2018