The draft genome of the specialist flea beetle Altica viridicyanea (Coleoptera: Chrysomelidae)

doi:10.21203/rs.3.rs-47391/v2

Download PDF

Research article

The draft genome of the specialist flea beetle Altica viridicyanea (Coleoptera: Chrysomelidae)

https://doi.org/10.21203/rs.3.rs-47391/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 07 Apr, 2021

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background: Altica (Coleoptera: Chrysomelidae) is a highly diverse and taxonomically challenging flea beetle genus that has been used to address questions related to host plant specialization, reproductive isolation, and ecological speciation. To further evolutionary studies in this interesting group, here we present a draft genome of a representative specialist, Altica viridicyanea, the first Alticinae genome and the fifth chrysomelid genome reported thus far.

Results: The genome is 864.8 Mb and consists of 4,490 scaffolds with a N50 size of 557 kb, which covered 98.6% complete and 0.4% partial insect Benchmarking Universal Single-Copy Orthologs. Repetitive sequences accounted for 62.9% of the assembly, and a total of 17,730 protein-coding gene models and 2,462 non-coding RNA models were predicted. To provide insight into host plant specialization of this monophagous species, we examined the key gene families involved in chemosensation, detoxification of plant secondary chemistry, and plant cell wall-degradation.

Conclusions: The genome assembled in this work provides an important resource for further studies on host plant adaptation and functionally affiliated genes. Moreover, this work also opens the way for comparative genomics studies among closely related Altica species, which may provide insight into the molecular evolutionary processes that occur during ecological speciation.

Epigenetics & Genomics

Altica

genome

annotation

host plant adaption

chemosensory

detoxification

The high rate of diversification among host-specific herbivorous insects is thought to result from their shift and specialization to distinct host-plant species, creating conditions that promote reproductive isolation and contribute to the process of speciation [1-3]. One herbivore group that has been used to address questions about host plant specialization, reproductive isolation, and ecological speciation is the leaf beetle genus Altica (Coleoptera: Chrysomelidae) [4-10]. This group has undergone rapid divergence that is largely associated with host plant use. For example, studies of three closely related species, Altica viridicyanea, A. cirsicola, and A. fragariae, have demonstrated that although these species are broadly sympatric and quite similar in morphology, they feed on distantly related host plants from different plant families [6]. Consequently, their divergence is likely the result of dietary shifts to unrelated host plants [6]. Further, their close relationship is supported by crossing studies that show that interspecific hybrids can be generated under laboratory conditions [4-6, 8, 10], and phylogenetic analysis indicates that these species diverged over a relatively short period of time. Although Altica has been pivotal for understanding the linkages between host plant use and speciation [6, 11-13], the lack of a representative Altica genome hinders our ability to more thoroughly investigate the molecular mechanisms underlying the processes of ecological adaptation and diversification within this interesting group.

Key aspects of host plant adaptation and speciation among Altica beetles involve at least two major evolutionary changes. First, a host plant shift will require behavioral adaptations that alter both feeding preference and oviposition preference for the novel host [14-15]. Both the adults and larvae of Altica species feed on the host plant, and this has been posited as a mechanism of reproductive isolation in herbivorous insects as an insect is more likely to encounter mates while feeding [7, 16-17]. Second, the insects will also require physiological adaptations that allow them to feed on new plants with different secondary compounds. This may involve new detoxification mechanisms that allow insects to avoid the deleterious effects of defensive chemistry [15, 18]. One aspect of host plant adaptation, then, is the prediction that the gene families involved in the detection of host chemical cues and those involved in xenobiotic detoxification will be key in facilitating successful shifts onto new host plant species. As a result, if we are to understand how host plant adaptation has played a role in speciation of Altica beetles, a reference genome would be helpful in making comparisons of candidate gene families involved in the diversification process.

The wealth of genetic and behavioral studies that have been conducted on A. viridicyanea makes it an excellent starting point for genomic investigations of host plant adaptation and speciation among Altica species. This species is an extreme host specialist of the plant Geranium nepalense (Sweet) (Geraniaceae), and these beetles will not feed or lay eggs on alternative hosts even in no-choice laboratory trials [4, 5, 7]. Reproductive isolation of A. viridicyanea is driven by the presence of species-specific cuticular hydrocarbons that determine mating preferences [8, 9]. In fact, these chemical signals are so specific that they allow beetles to assess sexual maturity of conspecifics [9]. Studies of F₁ hybrids between A. viridicyanea and a close relative have shown that these cuticular hydrocarbon profiles can also be modified by the beetle’s diet [8, 10]. Consequently, host plant use and mating preferences are intrinsically linked through chemistry.

Here we provide the genome assembly of A. viridicyanea (Fig. 1), the first genome of the subfamily Alticinae, and the fourth for the Chrysomelidae. Chrysomelids are a large and highly diverse family of beetles, with approximately 40,000 described species in the world [19], many of which are economically important pests of agricultural crops [20]. To date, genomes have been assembled for four species: the Colorado potato beetle, Leptinotarsa decemlineata [21], the cowpea weevil, Callosobruchus maculatus [22], the ragweed leaf beetle, Ophraella communa [23], and the western corn rootworm, Diabrotica virgifera virgifera (https://www.ncbi.nlm.nih.gov/genome/?term=Diabrotica+virgifera+virgifera). All of these species exhibit intermediate preferences in host diet and are restricted to feeding on a single plant family (oligophagous). Consequently, the genome for A. viridicyanea will add to our genomic resources for a monophagous (restricted to a single host plant species) member of the Chrysomelidae. Furthermore, this assembly will also expand our knowledge of beetles in general as we currently have only 21 published beetle genome assemblies [22, 24-33] (Table S1) which is a comparably small number for a diverse insect group with nearly 400,000 described species [34-36]. Finally, the genome assembled here will provide an important resource for further studies on host plant adaptation and functionally affiliated genes.

De novo genome assembly

Flow cytometry revealed that the genome size of A. viridicyanea ranged from 836.3±13.8Mb in females to 795.7±8.3Mb in males. We generated and assembled 187.3× coverage from Illumina short reads and 72.7× coverage via PacBio long reads from 157 female adults, thus creating a draft genome reference assembly of 864.8 Mb consisting of contig and scaffold N50s of 92.8 kb and 557.2 kb, respectively. The GC content was 31.67%. The size of the A. viridicyanea genome was larger than 85% of the currently published beetle genomes (Table S1). The draft genome assembly of A. viridicyanea was contained within 17,580 contigs that were assembled into 4,490 scaffolds, with the longest scaffold size of 5.6 Mb. Using the reference set of 1,658 insect BUSCOs, the genome contains 98.6% complete single-copy orthologs and multi-copy orthologs; using the reference set of 2,442 Endopterygota BUSCOs, our genome contains 95.8% complete single-copy orthologs and multi-copy orthologs (Table 1). Together, the results of the above analyses indicate that the genome of A. viridicyanea is a robust assembly. The estimated heterozygosity in the Illumina reads was about 0.70% ~ 0.96%, depending on k-mer size (k-mer 17, 19, 21, 23, 25 and 27).

We used these PacBio RNA-seq data to evaluate the genome assembly. Of the 13,550 polished reads, 60.46% could be successfully mapped to the genome. Transcripts that were unmapped or mapped with coverage or identity below the minimum threshold were partitioned into 1,177 gene families based on k-mer similarity. Of the 1,177 gene families, 121 could be mapped to the assembled genome, and hits to sequences from other species were found for 13 gene families. Blastn revealed that 1,032 of the remaining gene families were similar to sequences from plants which may represent DNA contamination by plant material during the DNA isolation step.

Genome annotation

Prior to gene prediction using the assembled sequences, repeat sequences were identified in the genome of A. viridicyanea. The repetitive sequence content was about 62.91 % of the assembly, which was similar to that of the cowpea weevil Callosobruchus maculatus (64%), lower than that of the ladybird Propylea japonica (71.33%), and much higher than that of other beetle species (Table S1). Most of the repetitive sequences were transposable elements. According to a uniform classification system for eukaryotic transposable elements [37], retrotransposons (Class I) accounted for 41.27% whereas DNA transposons (Class II) accounted for 26.24% of the genome (Table 2).

To check whether the PacBio reads could span most of the repeats (transposons here), we aligned the PacBio reads to the assembled genome. Focusing on the primary alignment only, there were 89.01% (5,792,616/6,507,752) of reads successfully mapped to the genome. We found that 99.33% (1,913,673/1,926,492) annotated transposons are fully covered by at least one read. Of these regions fully spanned by PacBio reads, the longest one was 29,177 bp. The result shows the length distribution of repeats fully covered and repeats not fully covered by reads. It is clear that repeats fully covered by reads are significantly shorter than those repeats not fully covered. So, we suggested longer reads could help to resolve these regions.

The integration of de novo, RNA-seq-based and homology-based gene prediction methods identified 17,730 protein-coding genes in A. viridicyanea (Table 3, Fig. S1), a number slightly less than the average of beetle species with available genomes (~18,600 genes on average, Table S1). In total, 16,625 genes were assigned to putative functions, accounting for approximately 93.77% of the predicted genes (Table S2), and 750 putative pseudogenes were identified (Table S3). There were 2,462 non-coding RNA models identified, including 45 miRNAs, 1093 rRNAs, and 1324 tRNAs, corresponding to 32, four and 24 gene families, respectively.

Phylogenetic analysis

We estimated the phylogenetic relationships of A. viridicyanea samples and an additional nine representative beetle species (Anoplophora glabripennis, Aethina tumida, Agrilus planipennis, Dendroctonus ponderosae, Diabrotica virgifera virgifera, Leptinotarsa decemlineata, Nicrophorus vespilloides, Onthophagus taurus and Tribolium castaneum). In total, 14,854 orthologs in A. viridicyanea clustered with the other nine representative beetle species. We identified 1,321 A. viridicyanea specific genes, corresponding to 470 gene families, and with the exception of Diabrotica virgifera virgifera, this number was much greater than the other representative beetle species included in this analysis (Table S4). The phylogenetic relationships were consistent with the results inferred from large datasets [38-40] based on 1,751 conserved single copy orthologs. For example, A. viridicyanea, Diabrotica virgifera virgifera and Leptinotarsa decemlineata, all belonging to chrysomelids, formed a clade, and these species clustered with Anoplophora glabripennis, a member of the superfamily Chrysomeloidea (Fig. 2). The estimated divergence time between A. viridicyanea and Diabrotica virgifera virgifera was about 74.7 million years ago. From this analysis, we also identified 155 gene families that expanded and 27 gene families that contracted along the A. viridicyanea lineage (Fig. 2). Some of these gene families were related to chemosensory and detoxification functions.

Chemosensory gene families

In many herbivorous insects, feeding, mating and oviposition behaviors are mediated by chemical cues [41]. The chemosensory system may also play important roles in speciation of some insects [42-44]. This is likely the case in A. viridicyanea as previous work has shown that this highly specialized beetle primarily uses chemical cues to achieve sexual isolation from its sibling species [8]. Furthermore, these contact chemicals also act as a mating signal to discriminate intraspecific variation in sexual maturity [9]. In addition, chemical cues are modified by and likely involved in host plant choice [8, 10]. Consequently, we investigated A. viridicyanea gene families known to be involved in chemosensory signaling in insects.

There are at least five gene families involved in the detection of chemicals, including three receptor families, odorant receptors (ORs), gustatory receptors (GRs) and ionotropic receptors (IRs), and two protein binding families, odorant binding proteins (OBPs) and chemosensory proteins (CSPs). These receptor families are usually expressed in insect olfactory sensory neurons and are involved in the detection of a suite of chemicals. For instance, volatile chemicals are detected by ORs [45-47], contact chemicals or carbon dioxide are detected by GRs [48], and nitrogen-containing compounds, acids, and aromatics are identified by IRs [49]. In contrast, the binding protein gene families are highly abundant in the sensillar lymph of insects and usually function as carriers of hydrophobic scent molecules to the receptors [50, 51].

In the genome of A. viridicyanea, we identified 195 putative chemosensory genes and two pseudogenes. Perhaps not surprisingly, the gene repertoire of the monophagous A. viridicyanea was considerably reduced as compared to that of host generalist species such as T. castaneum (630 genes plus 103 pseudogenes) and A. glabripennis (451 genes plus 65 pseudogenes). Upholding this pattern, A. viridicyanea has fewer chemosensory genes than the oligophagous species such as Dendroctonus ponderosae (240 genes plus 10 pseudogenes) and L. decemlineata (>300 genes) that specialize on a single family of host plants (Table S5). Yet there is a single outlier to this trend; Agrilus planipennis (132 genes and two pseudogenes), a species that is intermediate in host range, has fewer chemosensory genes than A. viridicyanea. These findings are generally consistent with the hypothesis that chemosensory gene content and host specificity should correlate in phytophagous beetles [52], although there is clearly an exception to this rule.

Insect ORs are proteins with seven transmembrane domains that are involved in the detection of volatile chemicals [46, 53, 54]. The number of ORs in beetle species varies widely from 30 to hundreds of ORs [55]. When we examined the A. viridicyanea genome for the presence of ORs, we found a diversity of gene families. There were 64 ORs and one pseudogene (PseudoGene48) that were classified into eight subfamilies: Group 1, 2A, 2B, 3, 4, 5A, 5B and 7 (Fig. 3; Table S5). Following the new OR classification scheme [55], we also identified one highly conserved olfactory co-receptor, Orco, that has been found in other beetle species. Interestingly, we also found a large expansion in A. viridicyanea in Group 4 that contained 17 ORs. By comparison, no more than four Group 4 OR genes have been previously identified in any other surveyed beetle species [52, 55].

In addition to ORs, we also compared GRs across beetle taxa. Most GRs are expressed in gustatory receptor neurons in taste organs and are involved in contact chemoreception and detection of CO₂ [48]. We annotated 38 GRs in A. viridicyanea, including three conserved candidate CO₂ receptors, seven candidate sugar receptors, and the remaining were candidate bitter taste receptors. Simple orthology of GRs is generally rare in beetles [52], and not surprisingly, no single-copy orthologs were revealed in the species that we compared. The phylogenetic analysis showed that 2-7 GRs from each of the seven species grouped within the clade of conserved sugar receptors. Additionally, two or three genes from five of the seven species formed a clade of CO₂ receptors (Fig. S2).

The number of GRs varied from 10 to 245 in the ten surveyed beetles (Table S5). Comparisons with A. viridicyanea identified as many as 147 GRs in an oligophagous chrysomelid species Leptinotarsa decemlineata, whereas fewer than 20 GRs were annotated in four other chrysomelids (Colaphellus bowringi, Ophraella communa, Pyrrhalta aenescens and Pyrrhalta maculicollis). The extremely low numbers of GRs in the latter four species is likely the result of differences in data collection—those species only had transcriptomic data available, and that approach generally does not describe the full complement of chemosensory genes. For example, a study in the longhorn beetle Anoplophora glabripennis found 11 GRs when using transcriptomic data, however, genomic data revealed 234 GRs [36, 56].

The next chemosensory receptor group that we examined was the IRs, a conserved family that evolved from a family of synaptic ligand-gated ion channels, ionotropic glutamate receptors (iGluRs) [49, 57, 58]. In insects, the IRs include two groups: the conserved “antennal IRs” that have an olfactory function, and the species-specific “divergent IRs” which are candidate gustatory receptors [59]. Our genome annotations revealed 33 ionotropic receptors (IRs). Only the members of the conserved antennal IR21a group were identified in all seven of the beetle species that we surveyed, whereas the clades IR8a, IR25a and IR76b were formed by single-copy orthologs from five species, excluding P. aenescens and P. maculicollis (transcriptomic data are available for both of these species) (Fig. 4). Furthermore, IRs from all seven species fell within the well-supported non-single-copy IR75 clade (Fig. 4). Remarkably, for A. viridicyanea, 18 of 33 IRs were clustered into a “Galerucinae+Alticinae” specific lineage that consisted of genes from O. communa, P. aenescens, P. maculicollis and A. viridicyanea (Fig. 4). These four species belong to the sister groups, Galerucinae and Alticinae, two closely related subfamilies of Chrysomelidae [60].

Finally, we examined the protein binding gene families. OBPs and CSPs are generally regarded as carriers of pheromones and odorants in insect chemoreception, and a multitude of additional functions have also been suggested such as carrying semiochemicals and visual pigments, promoting development and regeneration, and digesting insoluble nutrients [61]. OBPs are small, soluble proteins with six conserved cysteines [50]. Although the detailed mechanisms remain unclear [62], it is believed that OBPs deliver hydrophobic molecules to the receptors [50]. In A. viridicyanea, we annotated 48 putative OBP genes and one pseudogene (PseudoGene855). Among these, 34 genes belonged to the Minus C OBPs. We found four clades of classic OBPs, which include single-copy orthologs from each of the seven species in the analysis, e.g., Classic I, IV, VIII and IX. In clades VII and X, two copies from Dendroctonus ponderosae were clustered with single-copy orthologs from other six species. The clades of Classic II, III, V and VI were formed by single-copy orthologs from 5-6 species (Fig. S3). Plus-C OBPs were not found in A. viridicyanea, and are also absent in the Pyrrhalta species that belong to the “Galerucinae+Alticinae” taxonomic group.

CSPs are characterized by the presence of four cysteines that form two disulfide bridges [63]. We annotated 12 CSP genes in A. viridicyanea. The phylogenetic analysis revealed that only one clade (clade 1) was formed by single-copy orthologs from the eight species surveyed. Clades 2-7 were formed by single-copy orthologs from 5-7 beetle species. In these lineages, the absence of IRs from transcriptomic sources (e.g., P. aenescens, P. maculicollis and O. communa) was more common whereas the orthologs of A. viridicyanea also lacked members of clade 5 (Fig. S4).

Similar to previous work on GRs, transcriptomes often fail to describe the full set of chemosensory genes due to low expression, spatiotemporal variation in expression, or shallow sequencing depth. For instance, 106 chemosensory genes were detected in Anoplophora glabripennis using transcriptomic sequencing [56] whereas more than 500 chemosensory genes (65 pseudogenes included) were annotated from its genome [52].

Detoxification supergene families

Novel plant secondary compounds often present a challenge for herbivorous insects, and physiological adaptation to novel plant secondary metabolites is a key problem. The detoxification and metabolism of most xenobiotics occurs via a common set of detoxification-related enzymes, all of which belong to multigene families [64]. The cytochrome P450s (P450s), carboxyl/cholinesterases (CCEs), and glutathione S-transferases (GSTs) are widely regarded as the major insect gene/enzyme families involved in xenobiotic detoxification [65-67]. In addition, the UDP-glucuronosyltransferases (UGTs) and ATP binding cassette transporters (ABCs) can also play a role in detoxification [68-70]. This diversity of detoxification enzymes is critical for many herbivorous insects [18, 71] as their diets often contain a suite of plant chemicals that can be toxic, reduce palatability, or slow development time.

The host plant of A. viridicyanea is Geranium nepalens which has a number of chemical defenses such as tannins, flavonoids and organic acids [72]. As a strict specialist, then, A. viridicyanea likely has adaptations that allow them to detoxify these chemicals. Indeed, we annotated 225 detoxification enzymes spanning all three families (101 P450s, 97 CCEs and 27 GSTs). Expansion and contraction of these gene families are considered important in adaptive phenotypic diversification [73]. Furthermore, meta-analyses have established that the size of the P450, CCE and GST gene families are correlated with insect diet breadth [66, 67, 74]. In contrast with these studies, we showed that although A. viridicyanea has a greater number of detoxifying genes than that of the closely-related oligophagous Leptinotarsa decemlineata (197 genes) [21, 67], it has fewer detoxifying genes than generalist T. castaneum (275 genes) [26, 75].

Insect cytochrome P450 proteins are important in both xenobiotic detoxification and synthesis and degradation of endogenous molecules such as ecdysteroids and juvenile hormone [76-78]. In insects, the cytochrome P450 family is divided into four major clades: the mitochondrial P450 clade, the CYP2 clade, the CYP3 clade, and the CYP4 clade [79]. We found 101 P450s in Altica viridicyanea spanning all four clades: five in the mitochondrial clade, seven in the CYP2 clade, 53 in CYP3 clade, and 36 in CYP4 clade (Fig. 5, Table S6). We found that a majority of these genes belonged to the CYP6 and CYP9 subfamilies of the CYP3 clade and the CYP4 subfamily of the CYP4 clade (Table S7). These P450 subfamilies are known to be involved in detoxification of plant allelochemicals as well as resistance to pesticides [21, 80-82].

In addition to the cytochrome P450s, the A. viridicyanea genome also contained 97 genes encoding putative CCEs (Fig. S5), which is slightly fewer than that of L. decemlineata (102), but more than that of the other eight beetle species that were included in the analysis (ranged from 44 to 82) [67]. The dietary/detoxification group included two clades: coleopteran xenobiotic metabolizing CCE (clade A) and ɑ-esterase type CCEs (clade B) [83]. In A. viridicyanea, there is a noteworthy expansion (62 genes) in clade A, whereas we did not identify any genes from Clade D (integument esterase), F (juvenile hormone esterase), or I (unknown function) (Fig. S5; Table S8).

Another group of detoxification enzymes that we examined are the GSTs. GSTs are involved in many cellular physiological activities, such as detoxification of endogenous and xenobiotic compounds, intracellular transport, biosynthesis of hormones and protection against oxidative stress [75, 84]. Insect GSTs are divided into two major groups, the cytosolic and the microsomal GST genes. The cytosolic group is further divided into six classes: Delta, Epsilon, Sigma, Omega, Theta, and Zeta [85]. The Delta and Epsilon classes are thought to be insect-specific [75, 84, 86], and members of the Epsilon subfamily are commonly involved in detoxification of xenobiotics [87]. We detected a total of 27 GST genes in A. viridicyanea (Fig. S6; Table S9). Both the total number and the number of detoxification-related Epsilon subfamily in A. viridicyanea were lower than that of most beetles [67]( Table S9).

UDP-glycosyltransferases (UGT) catalyze the conjugation of a range of diverse small lipophilic compounds with sugars to produce glycosides, playing an important role in the detoxification of xenobiotics and in the regulation of endobiotics in insects [68]. From 17 (Oryctes borbonicus) to 65 (Anoplophora glabripennis) UGTs were identified in the nine beetle species surveyed [3, 67]. Currently, the largest repertoire of UGTs in beetles was found in the polyphagous longhorn beetle Anoplophora glabripennis, with 65 putative UGT genes and 7 pseudogenes [36]. The expansion of UGTs in A. glabripennis is thought to be related to its ability to feed on a broad range of host plants [36]. In line with this, we annotated 32 UGTs in the A. viridicyanea genome. A number of UGT50s were identified in this species, which has been suggested as the most conserved UGT in insects [68], and we also observed a remarkable expansion in the UGT324 family (Fig. S7, Table 10).

Most ABC proteins engage in active transport of molecules across cell membranes. The ABC transporters are well-known components of various detoxification mechanisms across all phyla [88, 89]. In the present study, we identified 69 putative ABCs in A. viridicyanea, belonging to eight subfamilies (A to H). This is a similar number to two specialist species of chrysomelids, Chrysomela populi and Diabrotica virgifera virgifera, (65 in each based on transcriptomic data) (Table S11). The gene numbers of the conserved subfamilies D, E and F were consistent with other beetles analyzed (Table S11); however, the number of genes in subfamilies B and C in A. viridicyanea (46) are the highest among the five species with which we compared (Table S11; Fig. 6). These subfamilies are known to be involved in detoxification processes [64, 90].

Plant cell wall-degrading enzymes

Early views of insect digestion postulated that insects lack the endogenous enzymes required for plant cell wall (PCW) digestion, and that PCW digestion by insects depended on exogenous enzymes from symbiotic microorganisms [91]. Recent studies, however, have revealed that endogenous PCW degrading enzymes are present in many insects and are important in the digestion of cellulose, hemicelluloses, and pectin in PCW [40, 92]. In fact, these enzymes are likely a key innovation in the adaptive radiation of herbivorous beetles. Some insect PCW-degrading enzymes are also involved in immune-defense responses and detoxification [40, 92].

Beetle-encoded plant cell wall-degrading enzymes are carbohydrate esterases (CE), polysaccharide lyases (PL), and mainly glycoside hydrolases (GH) [40]. In A. viridicyanea, we identified 65 putative glycoside hydrolases, including 35 GH1 genes, 10 GH45 genes, two GH48 genes and 18 GH28 genes (Table S12). Genes of GH1 originated anciently in animals and are ubiquitous in beetles. The species of Phytophaga (i.e., Chrysomeloidea and Curculionoidea) examined thus far have a greater number of GH1 genes than A. viridicyanea, for example, 228 were found in Diabrotica undecimpunctata, 135 in Mastostethus salvini, and 136 in Rhynchitomacerinus kuscheli [21, 36, 40]. For another ancient and ubiquitous gene family, GH9, there are at least a dozen independent losses in beetles [40]. GH9 was not detected in A. viridicyanea, along with 4 of 7 other chrysomelid species (Callosobruchus maculatus, Donacia marginata, Diabrotica undecimpunctata and Leptinotarsa decemlineata) [21, 40].

The other plant cell wall-degrading gene families (CE8, PL4, GH32, GH5, GH10, GH43, GH44, GH45, GH48 and GH28) are suggested to be obtained from bacteria and fungi via horizontal gene transfer, and are mainly found in Buprestoidea and Phytophaga, with scattered genes in a few other taxa [40]. In A. viridicyanea, in addition to the ubiquitous GH1 genes, three families of PCW-degrading enzymes were identified, including cellulose degrading GH45 and GH48, and pectin degrading GH28. These observed gene numbers are similar to that of closely related species in the Chrysomelinae (Oreina cacaliae and Leptinotarsa decemlineata) and Galerucinae (Diabrotica undecimpunctata) [21, 40].

In the present study, we combined long reads of PacBio with the higher fidelity of the short reads generated with Illumina sequencing. From the genome annotation we found that A. viridicyanea, a host specialist herbivore, has a reduced number of chemosensory and detoxification genes as compared to more generalist herbivorous beetles, consistent with the idea that diet breadth should positively correlate with chemosensory and detoxification gene content. Although A. viridicyanea had fewer chemosensory and detoxification genes than more polyphagous beetles, we did observe expansions in some gene families that may be related to host plant adaptation. As a result, the genome assembled here provides an important resource for further studies on host plant adaptation and functionally affiliated genes. Additionally, this work will also open the opportunity for comparative genomics studies among closely related Altica species that may provide insights into the molecular evolutionary processes that occur during ecological speciation.

Beetles

Altica viridicyanea is a highly specialized herbivore of Geranium nepalense. They are elongate-ovate beetles with a length of 3-4 mm. The dorsal surface is black with a metallic blue reflection (Fig. 1). Equipped with dilated hind legs, these beetles typically jump in a flea-like fashion to escape predators.

To make genome sequencing easier, we created a laboratory colony by collecting adult A. viridicyanea in Changping (40.28’N, 116.05’E), Beijing, China. Adults were maintained in growth chambers held at 16:8 LD and 25ºC and fed with leaves of their host plant, Geranium nepalense (Sweet). A subset of these collected adults was used for genome size estimation (see below). We maintained the colony through successive single-pair sibling matings to create third generation lines for whole genome sequencing. When we were ready to sequence the beetles, beetles were starved for two days and then killed in liquid nitrogen. The samples were stored at -80°C until DNA extraction. In total, we used 157 virgin females for sequencing.

Genome size estimation

In preparation for whole genome sequencing, we first determined the genome size of A. viridicyanea. Genome size was estimated via flow cytometry [93] on four adult males and four adult females. The thoracic muscle of living beetles was dissected with sterilized fine forceps, and cut into small pieces in Galbraith buffer. The ground suspension was filtered through 40-mm nylon mesh (Easystrainer^TM) to remove cellular debris and the flow-through was collected in a 5-ml round-bottomed tube placed on ice. Propidium iodide was add to samples to a final concentration of 50 mg/ml and stained in the dark at 4°C for 2 hours. Then fluorescence intensity was estimated for each beetle using a superfluid cell sorting system (MoFol XDP, Beckman Coulter Life Sciences). Genome size was calculated by comparing samples to an internal reference subsisting of chicken blood.

DNA extraction and sequencing

Genomic DNA (gDNA) was extracted from whole beetles using DNeasy Blood & Tissue Kit (Qiagen, Germany) following the manufacturer’s instructions. Samples were first surface sterilized using 75% ethanol and sterile deionized water. The total amount of gDNA was measured using QubitFluorometer (Invitrogen), and the integrity of the gDNA was verified on an agarose gel that had reference lanes containing high molecular weight ladders (GeneRuler High Range DNA Ladder and D2000 DNA Ladder). From these extractions, seven Illumina sequencing libraries were prepared, with insert sizes and genome coverage of 270 bp (55.5×), 500 bp (29.4×), 800 bp (19.0×), 3 kb (12.5×), 5 kb (two libraries, 49.7×), to 10 kb (21.1×) (Table S13). The libraries were sequenced with 150 bp paired-end reads on the Illumina HiSeq 2500 platform. We also sequenced two PacBio libraries with 72.7× genome coverage, with N50 read lengths of 9.1 kb (Table S13). The first PacBio DNA library (20 kb) was constructed using the PacBio SMRTbell Template Prep Kit 3.0 (Pacific Biosciences, Menlo Park, CA, USA) and sequenced on a PacBio RS II sequencer with the P6 polymerase/C4 chemistry combination at 1GENE (Zhejiang, China). A total of 16 SMRT Cells were processed and the movie length was six hours. The second PacBio DNA library (20 kb) was constructed using the SMRTBell template preparation kit 1.0 (PacBio, USA), for which six SMRT Cells were run on the PacBio Sequel instrument at BGI (Guangdong, China) with SequelTM Sequencing Kit 2.1 (PacBio, USA). The movie length was 10 hours.

RNA-seq library construction and sequencing

We produced transcriptomic data to facilitate gene prediction analyses. To obtain more full-length or near full-length gene sequences, we generated sequence data by combining PacBio full-length transcriptome sequencing (8 Gb clean data) and Illumina sequencing (10 Gb clean data). Whole bodies of four males and four females with three libraries (size-selection: 1-2 kb, 2-3 kb and >3 kb lengths) were used for PacBio sequencing. The PacBio transcriptome libraries were constructed using PacBio SMRTbell Template Prep Kit 3.0 (Pacific Biosciences, Menlo Park, CA, USA) and sequenced on a PacBio RS II sequencer at 1GENE (Zhejiang, China). The sequencing chemistry was P6-C4, and the movie length was four hours. Five SMRT Cells were processed: one for the 1-2 kb library, two for the 2-3 kb library and two for the >3 kb library. To enrich for chemosensory genes, heads from 30 females and 30 males were also used for Illumina sequencing.

Data pre-processing

For Illumina Data, FastQC v0.11.6 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to check the quality of the raw reads under default settings, and then TrimGalore v0.4.5 (https://github.com/FelixKrueger/TrimGalore) was used to filter out the adapters and low-quality reads using the parameters “--length 36 -q 20 --trim-n”. Because A. viridicyanea are quite small, we used the whole body of adults for sequencing. Although we surface sterilized the insects, the likelihood of microbial contamination from gut contents was high; therefore, we used BBDuk in the BBMap v37.80 (https://sourceforge.net/projects/bbmap) with parameters “ordered=t k=31” to filter microbial contaminants. To do this, we first built a reference library which included all sequences of archaea, bacteria, fungi, protozoa, and viruses available in RefSeq (accessed on August 15, 2019). The mitochondrion sequence of A. viridicyanea [94] was also included to remove reads originating from the mitochondrial genome. After these filtering steps, the remaining reads were used for downstream genome assembly and analysis.

To simplify the genome assembly process, we also excluded reads originating from microbial sources that were present in the PacBio data. To accomplish this, we mapped all the raw PacBio reads to the above Illumina library using minimap2 v2.12 [95], using the parameters “-x map-pb -a -Q --split-prefix”. To improve the mapping specificity, we included 11 available coleopteran genomes [28] in the reference library. Then we removed the reads that mapped to the microbial sequences.

For the PacBio Iso-Seq data processing, we used the new Iso-Seq3 pipeline (https://github.com/PacificBiosciences/IsoSeq_SA3nUP). First, we first converted the bax.h5 files of each SMRT Cell run to the BAM format using bax2bam v0.0.8 (https://github.com/PacificBiosciences/bax2bam) with default parameters. Circular Consensus Sequence (CCS) calling was then done by ccs v3.4.1 (https://github.com/PacificBiosciences/ccs) under default settings. To remove cDNA primers and obtain full-length reads, we used lima v2.0.0 (https://github.com/pacificbiosciences/barcoding) with parameters “--isoseq --peek-guess”. Next we used “isoseq3 refine --require-polya” to trim the poly (A) tails. Finally, to generate polished, high-quality reads, the SMRT cells were merged and we used “isoseq3 cluster” with parameters “--use-qvs”.

Heterozygosity estimation

Heterozygosity was estimated using gce v1.0.0 (parameters: -b 1 -H 1 -m 1 -D 8), based on the pooled Illumina reads of the sib-mated females. Before running gce, we first used kmer_freq_hash in the gce package to compute the histogram of kmer frequencies (k-mer from 17 to 25).

Genome assembly

A hybrid assembly strategy was used to assemble the PacBio and Illumina reads. PacBio reads were first assembled using Flye v2.3.5b [96] with the parameters set to “--genome-size 830m” to obtain contigs. Then Redundans v0.14a [97], which uses an iterative scaffolding approach by calling SSPACE [98], was used to scaffold the PacBio contigs with the Illumina short reads, with parameters set to “--noreduction”. After that, SSPACE-Long v1.1 [99] was used to further scaffold with PacBio long reads using default settings (Table S14).

In addition, we used these PacBio RNA-seq data to evaluate the genome assembly. We mapped PacBio RNA-seq reads to the genome using minimap2 v2.12, with parameters “-ax splice -uf --secondary=no”. For those unmapped and poorly mapped transcripts (coverage <0.99 or identity <0.95), the Cogent pipeline (https://github.com/Magdoll/Cogent) was employed to explore their potential functions. Transcripts were partitioned into gene families based on k-mer similarity. For each gene family, contigs were reconstructed to represent the “union” of all coding bases. To check the potential functions of these transcripts, we used Blastn v2.7.1 [100] against the non-redundant sequence database (ftp://ftp.ncbi.nih.gov/blast/db/FASTA; accessed at November 2020).

Estimation of genome completeness

Benchmarking Universal Single-Copy Orthologs (BUSCO v3.0.1) [101] analyses against Insecta (n=1,658) and Endopterygota (n=2,442) were used to evaluate the integrity and quality of the genome assembly with parameters "-sp tribolium2012 -m geno".

Genome annotation

Prior to gene prediction using the assembled sequences, we identified repeat sequences in the genome of A. viridicyanea. LTR FINDER v1.05 [102], RepeatScout v1.0.5 [103] and PILER-DF v2.4 [104] were used to construct a library of repetitive sequences based on the A. viridicyanea genome. PASTEClassifier v1.0 [105] was used to classify these repeats, and the Repbase database [106] was used to merge them. Finally, RepeatMasker v4.0.5 [107] was used to identify and mask the genomic repeated sequences. All the parameters were set to “default” except RepeatMasker where the parameters were set to “-nolow -no_is -norna -engine wublast”. To check whether the PacBio reads could span most of the repeats (transposons here), we aligned the PacBio reads to the assembled genome using minimap2 v2.12 (parameters: -x map-pb). Then Bedtools v2.26.0 [108] was employed to get the read coverage of each repeat.

We combined the de novo, homology-based, and transcriptome-based predictions to identify protein-coding genes in the A. viridicyanea genome. The default settings for Genscan v1.1.0 [109], Augustus v2.4 [110], GlimmerHMM v3.0.4 [111], GeneID v1.4 [112] and SNAP (v2006-07-28) [113] were used to generate the de novo predictions. For the homology-based prediction, protein sequences of Drosophila melanogaster, Tribolium castaneum, Dendroctonus ponderosae and Anoplophora glabripennis were downloaded from NCBI and aligned to the assembled A. viridicyanea genome with GeMoMa v1.3.1 [114, 115] using default parameters. For transcriptome sequencing-based prediction, we firstly assembled the Illumina short reads into unigenes using Hisat v2.0.4 with parameters “--max-intronlen 20000，--min-intronlen 20” [116] and StringTie v1.3.4 set to default parameters [116]. Second, we combined the unigenes with the PacBio full-length transcripts, then we predicted genes based on these combined sequences with PASA v2.0.2 using the parameters “-align_tools gma，-maxIntronLen 20000” [117]. EVidenceModeler v1.1.1 was used to integrate all of the gene predictions [118]. Finally, to obtain a final gene dataset, PASA v2.0.2 was used for further modification including adding 5’-UTR and 3’-UTRs, obtaining alternative splicing, and extending or shortening, adding or subtracting exons for the EVM-predicted results.

We predicted the non-coding RNAs based on the Rfam v12.1 [119] and miRBase v21.0 [120] databases. Putative microRNAs (miRNAs) and ribosomal RNAs (rRNAs) were predicted using Infernal v1.1 [121], and transfer RNAs (tRNAs) were predicted with tRNAscan- SE v1.3.1 [122]. Based on homology to known protein-coding genes, putative pseudogenes were first searched in the intergenic regions of the A. viridicyanea genome using genBlastA v1.0.4 [123]. Then GeneWise v2.4.1 [124] was used to search for premature stop codons or frameshift mutations in those sequences.

The predicted genes were annotated by aligning them to the NCBI non-redundant protein (nr) [125], non-redundant nucleotide (nt) [125], Swissprot [126], TrEMBL [123], KOG [127], and KEGG [128] databases using BLAST v2.2.31 [129] with a maximal e-value of 1e^-5. We assigned gene ontology (GO) terms [130] to the genes using the BLAST2GO v2.5 pipeline [131].

Species phylogenetic analysis

We estimated the phylogenetic relationships of A. viridicyanea and an additional nine representative beetle species for which genomic data were available (Anoplophora glabripennis, Aethina tumida, Agrilus planipennis, Dendroctonus ponderosae, Diabrotica virgifera virgifera, Leptinotarsa decemlineata, Nicrophorus vespilloides, Onthophagus taurus and Tribolium castaneum). We used 1,751 conserved single copy orthologs for this analysis. The analysis was implemented in PhyML v3.0 [132] with parameters “-gapRatio 0.5 -badRatio 0.25 -bootstrap 1000”.

We also estimated divergence times among the species using MCMCtree (PAML v4.8 package) (http://abacus.gene.ucl.ac.uk/software/paml.html) [133, 134] with parameters set to a burn-in time of 10,000, the sample number was 100,000, and the sampling frequency was 2. We used the time divergence data in timetree (http://www.timetree.org/) for calibration (Nicrophorus vespilloides - Onthophagus taurus [245~296 MYA]; Dendroctonus ponderosae - Anoplophora glabripennis [149~240 MYA]). To examine gene family expansion and contraction among species, we used CAFE v4.1 [135] with “lambda -l 0.002” to automatically search for birth and death parameters (λ) of genes.

Gene family evolution

To provide insights into host plant specialization, we specifically focused on the evolution of chemosensory gene families, chemical detoxification supergene families, and plant cell wall-degrading enzymes. We first conducted multiple sequence alignment using Mafft (online version 7.305) [136] (http:// mafft.cbrc.jp/alignment/server/), applying the L-INS-I algorithm. The best-fit models for amino acid sequence evolution were selected using the Akaike Information Criterion (AIC) in Prottest v3.4.2. [137] (Table S15). Maximum likelihood trees were reconstructed using RaxML v8.2.9 [138] in the CIPRES Science Gateway v3.3 (https://www.phylo.org/) [139], with 1000 non-parametric bootstrap replicates. The resulting tree was viewed and edited with FigTree v1.3.1 [140] and Adobe Illustrator CS5.

ABC: ATP binding cassette transporter; bp: base pairs; BUSCO: Benchmarking Universal Single-Copy Orthologs; CCEs: carboxyl/cholinesterase; CE: carbohydrate esterase; CSP: chemosensory protein; Gb: gigabase pairs; GH: glycoside hydrolase; GR: gustatory receptor; GST: glutathione S-transferase; IR: ionotropic receptor; kb: kilobase pairs; LINE: long interspersed nuclear element; LTR: long terminal repeat; Mb: megabase pairs; ML: maximum likelihood; OBP: odorant binding protein; OR: odorant receptor; P450: cytochrome P450; PCW: plant cell wall; PL: polysaccharide lyase; RAxML: Randomized Axelerated Maximum Likelihood; RNAseq: RNA sequencing; SINE: short interspersed nuclear element; TE: transposable element; UGT: UDP-glucuronosyltransferase.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The raw sequence data reported in this paper are available at the Genome Sequence Archive [141] in National Genomics Data Center [142] (https://bigd.big.ac.cn/gsa; accession number: CRA002741). And the whole genome sequence data have been deposited in the Genome Warehouse in the National Genomics Data Center [142] (https://bigd.big.ac.cn/gwh; accession number: GWHAMMQ00000000).

Competing Interests

The authors declare that they have no competing interests.

Funding

This study was supported by National Natural Science Foundation of China (Grant Nos. 31472030, 31672334).

Author’s Contributions

H.J.X, X.K.Y. and R.S.C. designed and managed this project. H.J.X. and W.Z.L. maintained the beetle colony used for sequencing. Y.W.N.,Y.J.H., L.L.Z, X.C.C, X.W.Z and H.J.X performed genome assembly and genome annotation. H.J.X. and R.E.N. performed the phylogenetic analysis. H.J.X., K.A.S and Y.W.N. wrote the paper.

Acknowledgements

Not applicable.

Forbes AA, Devine SN, Hippee AC, Tvedte ES, Ward AKG, Widmayer HA, et al. Revisiting the particular role of host shifts in initiating insect speciation. Evolution. 2017;71:1126–1137.
Nosil P, Crespi BJ, Sandoval CP. Host-plant adaptation drives the parallel evolution of reproductive isolation. Nature. 2002;417:440–443.
Simon J-C, d’Alençon E, Guy E, Jacquin-Joly E, Jaquiéry J, Nouhaud P, et al. Genomics of adaptation to host-plants in herbivorous insects. Brief Funct Genomics. 2015;14:413–423.
Xue HJ, Li WZ, Yang XK. Genetic analysis of feeding preference in two related species of Altica (Coleoptera: Chrysomelidae: Alticinae). Ecol Entomol. 2009;34:74–80.
Xue HJ, Magalhães S, Li WZ, Yang XK. Reproductive barriers between two sympatric beetle species specialized on different host plants. J Evol Biol. 2009;22:2258–2266.
Xue HJ, Li WZ, Nie RE, Yang XK. Recent speciation in three closely related sympatric specialists: inferences using multi-locus sequence, post-mating isolation and endosymbiont data. PLoS One. 2011;6:e27834.
Xue HJ, Li WZ, Yang XK. Assortative mating between two sympatric closely-related specialists: inferred from molecular phylogenetic analysis and behavioral data. Sci Rep. 2014;4:5436.
Xue HJ, Wei JN, Magalhães S, Zhang B, Song KQ, Liu J, et al. Contact pheromones of 2 sympatric beetle species are modified by the host plant and affect mating. Behav Ecol. 2016;27:895–902.
Xue HJ, Zhang B, Segraves KA, Wei JN, Nie RE, Song KQ, et al. Contact cuticular hydrocarbons act as a mating cue to discriminate intraspecific variation in Altica flea beetles. Anim Behav. 2016;111:217–224.
Xue HJ, Segraves KA, Wei J, Zhang B, Nie RE, Li WZ, et al. Chemically mediated sexual signals restrict hybrid speciation in a flea beetle. Behav Ecol. 2018;29:1462–1471.
Laroche A, DeClerck-Floate RA, LeSage L, Floate KD, Demeke T, et al. Are Altica carduorum and Altica cirsicola (Coleoptera: Chrysomelidae) different species? Implications for the release of cirsicola for the biocontrol of Canada thistle in Canada. Biol Control. 1996;6:306–314.
Jenkins TM, Braman SK, Chen Z, Eaton TD, Pettis GV, Boyd DW, et al. Insights into flea beetle (Coleoptera: Chrysomelidae: Galerucinae) host specificity from conconrdant mitochondrial and nuclear DNA phylogenies. Ann Entomol Soc Am. 2009;102:386–395.
Reid CA, Beatson M. Disentangling a taxonomic nightmare: a revision of the Australian, Indomalayan and Pacific species of Altica Geoffroy, 1762 (Coleoptera: Chrysomelidae: Galerucinae). Zootaxa. 2009;3918:503–551.
Bernays EA, Chapman RF. Host-plant Selection by Phytophagous Insects. Chapman and Hall, New York; 1994.
Gassmann AJ, Levy A, Tran T, Futuyma DJ. Adaptations of an insect to a novel host plant: a phylogenetic approach. Funct Ecol. 2006;20:478–485.
Katakura H, Shioi M, Kira Y. Reproductive isolation by host specificity in a pair of phytophagous ladybird beetles. Evolution. 1989;43:1045–1053.
Via S. Reproductive isolation between sympatric races of pea aphids. I. Gene flow and habitat choice. Evolution. 1999;53:1446–1457.
Heidel-Fischer HM, Vogel H. 2015. Molecular mechanisms of insect adaptation to plant secondary compounds. Curr Opin Insect Sci. 2015;8:8–14.
Leschen RAB, Beutel RG. Handbook of Zoology, Band 4: Arthropoda: Insecta, Teilband / Part 40: Coleoptera, Beetles, Morphology and Systematics (Phytophaga), Vol. 3. Walter de Gruyter, Berlin; 2014.
Jolivet PH, Cox ML, Petitpierre E. Novel aspects of the biology of Chrysomelidae. Kluwer Academic Publishers, Dordrecht/Boston/London; 1994.
Schoville SD, Chen YH, Andersson MN, Benoit JB, Bhandari A, Bowsher JH, et al. A model species for agricultural pest genomics: the genome of the Colorado potato beetle, Leptinotarsa decemlineata (Coleoptera: Chrysomelidae). Sci Rep. 2018;8:1931.
Sayadi A, Barrio AM, Immonen E, Dainat J, Berger D, Tellgren-Roth C, et al. The genomic footprint of sexual conflict. Nat Ecol Evol. 2019;3:1725–
Sarah B, Laurent F, Heinz M-S. Genome assembly of the ragweed leaf beetle: a step forward to better predict rapid evolution of a weed biocontrol agent to environmental novelties. Genome Biol Evol. 2020; https://doi.org/10.1093/gbe/evaa102.
Fu XH, Li JJ, Tian Y, Quan WP, Zhang S, Liu Q, et al. Long-read sequence assembly of the firefly Pyrocoelia pectoralis GigaScience. 2017;6:1–7.
Ando T, Matsuda T, Goto K, Hara K, Ito A, Hirata J, et al. Repeated inversions within a pannier intron drive diversification of intraspecific colour patterns of ladybird beetles. Nat Commun. 2018;9:3843.
Evans JD, McKenna D, Scully E, Cook SC, Dainat B, Egekwu N, et al. Genome of the small hive beetle (Aethina tumida, Coleoptera: Nitidulidae), a worldwide parasite of social bee colonies, provides insights into detoxification and herbivory. GigaScience 2018;7:1–16.
Fallon TR, Lower SE, Chang CH, Bessho-Uehara M, Martin GJ, Bewick AJ, et al. Firefly genomes illuminate parallel origins of bioluminescence in beetles. eLife. 2018;7:e36495.
McKenna DD. 2018. Beetle genomes in the 21st century: prospects, progress and priorities. Curr Opin Insect Sci. 2018;25:76–82.
Wu YM, Li J, Chen XS. Draft genomes of two blister beetles Hycleus cichorii and Hycleus phaleratus. 2018;7:1–7.
Kraaijeveld K, Neleman P, Mariën J, de Meijer E, Ellers J. Genomic resources for Goniozus legneri, Aleochara bilineata and Paykullia maculata, representing three independent origins of the parasitoid lifestyle in insects. G3-Genes Genom Genet. 2019;9:987–991.
Wang K, Li PP, Gao YY, Liu CQ, Wang QL, Yin J, et al. De novo genome assembly of the white-spotted flower chafer (Protaetia brevitarsis). GigaScience. 2019;8:1–9.
Guan DL, Hao XQ, Mi D, Peng J, Li Y, Xie JY, et al. Draft genome of a blister beetle Mylabris aulica. Front Genet. 2020;10:1281.
Zhang LJ, Li S, Luo JY, Du P, Wu LK, Li YR, et al. Chromosome-level genome assembly of the predator Propylea japonicato understand its tolerance to insecticides and high temperatures. Mol Ecol Res. 2020;20:292–
Hammond PM. Species inventory. In: Groombridge B, editor. Global Biodiversity, Status of the Earth’s Living Resources. London: Chapman and Hall; 1992. p. 17–
Slipinski SA, Leschen RAB, Lawrence JF. Order Coleoptera Linnaeus, 1758. In: Zhang, Z.-Q. (Ed.) Animal biodiversity: An outline of higher-level classification and survey of taxonomic richness. Zootaxa. 2011;3148:203–208.
McKenna DD, Scully ED, Pauchet Y, Hoover K, Kirsch R, Geib SM, 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. 2016;17:227.
Wicker T, Sabot F, Hua-Van A, Bennetzen JL, Capy P, Chalhoub B, et al. A unified classification system for eukaryotic transposable elements. Nat Rev Genet. 2007;8:973–
Hunt T, Bergsten J, Levkanicova Z, Papadopoulou A, John OS, Wild R, et al. A comprehensive phylogeny of beetles reveals the evolutionary origins of a superradiation. Science. 2007;318:1913–1916.
Mckenna DD, Wild AL, Kanda K, Bellamy CL, Beutel RG, Caterino MS, et al. The beetle tree of life reveals that Coleoptera survived end-Permian mass extinction to diversify during the Cretaceous terrestrial revolution. Syst Entomol. 2015;40:835–880.
McKenna DD, Shin S, Ahrens D, Balke M, Beza-Beza C, Clarke DJ, et al. The evolution and genomic basis of beetle diversity. Proc Natl Acad Sci U S A. 2019;116:24729–24737.
Andersson MN, Grosse-Wilde E, Keeling CI, Bengtsson JM, Yuen MMS, Li M, et al. Antennal transcriptome analysis of the chemosensory gene families in the tree killing bark beetles, Ips typographus and Dendroctonus ponderosae (Coleoptera: Curculionidae: Scolytinae). BMC Genomics. 2013;14:198.
Smadja C, Butlin RK. On the scent of speciation: the chemosensory system and its role in premating isolation. Heredity. 2009;102:77–97.
Smadja CM, Canbäck B, Vitalis R, Gautier M, Ferrari J, Zhou JJ, et al. Large-scale candidate gene scan reveals the role of chemoreceptor genes in host plant specialization and speciation in the pea aphid. Evolution. 2012;66:2723–2738.
Zhang B, Zhang W, Nie RE, Li WZ, Segraves KA, Yang XK, et al. Comparative transcriptome analysis of chemosensory genes in two sister leaf beetles provides insights into chemosensory speciation. Insect Biochem Mol Biol. 2016;79:108–118.
Hallem EA, Carlson JR. Coding of odors by a receptor repertoire. Cell. 2006;125:143–160.
Carey AF, Wang G, Su CY, Zwiebel LJ, Carlson JR. Odorant reception in the malaria mosquito Anopheles gambiae. Nature. 2010;464:66–71.
Stensmyr MC, Dweck HKM, Farhan A, Ibba I, Strutz A, Mukunda L, et al. A conserved dedicated olfactory circuit for detecting harmful microbes in Drosophila. Cell. 2012;151:1345–1357.
Vosshall LB, Stocker RF. Molecular architecture of smell and taste in Drosophila. Ann Rev Neurosci. 2007;30:505–533.
Abuin L, Bargeton B, Ulbrich MH, Isacoff EY, Kellenberger S, Benton R, et al. Functional architecture of olfactory ionotropic glutamate receptors. Neuron. 2011;69:44–60.
Sánchez-Gracia A, Vieira FG, Rozas J. 2009. Molecular evolution of the major chemosensory gene families in insects. Heredity. 2009;103:208–216.
Zhou JJ. Odorant-binding proteins in insects. Vitam Horm. 2010;83:241–272.
Andersson MN, Keeling CI, Mitchell RF. Genomic content of chemosensory genes correlates with host range in wood-boring beetles (Dendroctonus ponderosae, Agrilus planipennis, and Anoplophora glabripennis). BMC Genomics. 2019;20:690.
Benton R. On the ORigin of smell: odorant receptors in insects. Cell Mol Life Sci. 2006;63:1579–1585.
Nichols AS, Chen S, Luetje CW. Subunit contributions to insect olfactory receptor function: channel block and odorant recognition. Chem Senses. 2011;36:781–790.
Mitchell RF, Schneider TM, Schwartz AM, Andersson MN, McKenna DD.. The diversity and evolution of odorant receptors in beetles (Coleoptera). Insect Mol Biol. 2020;29:77–
Hu P, Wang JZ, Cui MM, Tao J, Luo YQ. Antennal transcriptome analysis of the Asian longhorned beetle Anoplophora glabripennis. Sci Rep. 2016;6:26652.
Benton R, Vannice KS, Gomez-Diaz C, Vosshall LB. Variant ionotropic glutamate receptors as chemosensory receptors in Drosophila. Cell. 2009;136:149–162.
Rytz R, Croset V, Benton R. 2013. Ionotropic Receptors (IRs): Chemosensory ionotropic glutamate receptors in Drosophila and beyond. Insect Biochem Mol Biol. 2013;43:888–897.
Croset V, Rytz R, Cummins SF, Budd A, Brawand D, Kaessmann H, et al. Ancient protostome origin of chemosensory ionotropic glutamate receptors and the evolution of insect taste and olfaction. PLoS Genet. 2010;6:e1001064.
Nie RE, Andújar C, Gómez-Rodríguez C, Bai M, Xue HJ, Tang M, et al. The phylogeny of leaf beetles (Chrysomelidae) inferred from mitochondrial genomes. Syst Entomol. 2020;45:188–204.
Pelosi P, Iovinella I, Zhu J, Wang GR, Dani FR. Beyond chemoreception: diverse tasks of soluble olfactory proteins in insects. Biol Rev. 2018;93:184–200.
Leal WS. Odorant reception in insects: roles of receptors, binding proteins, and degrading enzymes. Ann Rev Entomol. 2013;58:373–391.
Pelosi P, Zhou JJ, Ban LP, Calvello M. Soluble proteins in insect chemical communication. Cell Mol Life Sci. 2006;63:1658–1676.
Koenig C, Bretschneider A, Heckel DG, Grosse-Wilde E, Hansson BS, Vogel H.. The plastic response of Manduca sexta to host and non-host plants. Insect Biochem Mol Biol. 2015;63:72–85.
Despres L, David J-P, Gallet C. The evolutionary ecology of insect resistance to plant chemicals. Trends Ecol Evol. 2007;22:298–307.
Rane RV, Walsh TK, Pearce SL, Jermiin LS, Gordon KHJ, Richards S, et al. Are feeding preferences and insecticide resistance associated with the size of detoxifying enzyme families in insect herbivores? Curr Opin Insect Sci. 2016;13:70–76.
Rane RV, Ghodke AB, Hoffmann AA, Edwards OR, Walsh TK, Oakeshott JG. Detoxifying enzyme complements and host use phenotypes in 160 insect species. Curr Opin Insect Sci. 2019;31:131–138.
Ahn SJ, Vogel H, Heckel DG. Comparative analysis of the UDP-glycosyltransferase multigene family in insects. Insect Biochem Mol Biol. 2012;42:133–
Dermauw W, Van Leeuwen T. The ABC gene family in arthropods: comparative genomics and role in insecticide transport and resistance. Insect Biochem Mol Biol. 2014;45:89–
Merzendorfer H. ABC transporters and their role in protecting insects from pesticides and their metabolites. Adv Insect Physiol. 2014;46:1–72.
Ramsey JS, Rider DS, Walsh TK, De Vos M, Gordon KHJ, Ponnala L, et al. Comparative analysis of detoxification enzymes in Acyrthosiphon pisum and Myzus persicae. Insect Mol Biol. 2010;19(Suppl. 2):155–164.
He WT, Jin ZX, Wang BQ. Research progress of Geranium nepalense J Aerosp Med. 2011;21:1200–1202.
Hahn MW, De Bie T, Stajich JE, Nguyen C, Cristianini N. Estimating the tempo and mode of gene family evolution from comparative genomic data. Genome Res. 2005;15:1153–1160.
Kanost MR, Arrese EL, Cao XL, Chen YR, Chellapilla S, Goldsmith MR, et al. Multifaceted biological insights from a draft genome sequence of the tobacco hornworm moth, Manduca sexta. Insect Biochem Mol Biol. 2016;76:118–147.
Shi HX, Pei LH, Gu SS, Zhu SC, Wang YY, Zhang Y, et al. Glutathione S-transferase (GST) genes in the red flour beetle, Tribolium castaneum, and comparative analysis with five additional insects. Genomics. 2012;100:327–335.
Scott JG. Cytochromes P450 and insecticide resistance. Insect Biochem Mol Biol. 1999;29:757–777.
Helvig C, Koener JF, Unnithan GC, Feyereisen R. CYP15A1, the cytochrome P450 that catalyzes epoxidation of methyl farnesoate to juvenile hormone III in cockroach Corpora allata. Proc Natl Acad Sci U S A. 2004;101:4024–4029.
Rewitz KF, O’Connor MB, Gilbert LI. Molecular evolution of the insect Halloween family of cytochrome P450s: phylogeny, gene organization and functional conservation. Insect Biochem Mol Biol. 2007;37:741–753.
Feyereisen R. Evolution of insect P450. Biochem Soc Trans. 2006;34:1252–1255.
Li XC, Berenbaum MR, Schuler MA. Plant allelochemicals differentially regulate Helicoverpa zea cytochrome P450 genes. Insect Biochem Mol Biol. 2002;11:343–351.
Scully ED, Hoover K, Carlson JE, Tien M, Geib SM. Midgut transcriptome profiling of Anoplophora glabripennis, a lignocellulose degrading cerambycid beetle. BMC Genomics. 2013;14:850.
Zhu F, Moural TW, Nelson DR, Palli SR. A specialist herbivore pest adaptation to xenobiotics through upregulation of multiple cytochrome P450s. Sci Rep. 2016;6:20421.
Lü FG, Fu K., Li Q, Guo WC, Ahmat T, Li GQ. Identification of carboxylesterase genes and their expression profiles in the Colorado potato beetle Leptinotarsa decemlineata treated with fipronil and cyhalothrin. Pestic Biochem Phys. 2015;122:86–95.
Enayati AA, Ranson H, Hemingway J. Insect glutathione transferases and insecticide resistance. Insect Biochem Mol Biol. 2005;14:3–8.
Friedman R. Genomic organization of the glutathione S-transferase family in insects. Mol Phylogenet Evol. 2011;61:924–932.
Meyer JM, Markov GV, Baskaran P, Herrmann M, Sommer RJ, Rödelsperger C. Draft genome of the scarab beetle Oryctes borbonicus on La Réunion Island. Genome Biol Evol. 2016;8:2093–2105.
Han JB, Li GQ, Wan PJ, Zhu TT, Meng QW. Identification of glutathione S-transferase genes in Leptinotarsa decemlineataand their expression patterns under stress of three insecticides. Pestic Biochem Phys. 2016;133:26–34.
Broehan G, Kroeger T, Lorenzen M, Merzendorfer H. Functional analysis of the ATP-binding cassette (ABC) transporter gene family of Tribolium castaneum. BMC Genomics. 2013;14:6.
Strauss AS, Peters S, Boland W, Gretscher RR, Groth M, Boland W, et al. ABC transporter functions as a pacemaker for sequestration of plant glucosides in leaf beetles. eLife. 2013;2:e01096.
Liu SM, Zhou S, Tian L, Guo EE, Luan YX, Zhang JZ, et al. Genome-wide identification and characterization of ATP-binding cassette transporters in the silkworm, Bombyx mori. BMC Genomics. 2011;12:491.
Watanabe H, Tokuda G. Cellulolytic systems in insects. Ann Rev Entomol. 2010;55:609–632.
Caldeŕon-Cort́es N, Quesada M, Watanabe H, Cano-Camacho H, Oyama K. Endogenous plant cell wall digestion: A key mechanism in insect evolution. Ann Rev Ecol Evol Syst. 2012;43:45–71.
Hare EE, Johnston JS. Genome size determination using flow cytometry of propidium iodide-stained nuclei. Method Mol Biol. 2011;772:3–12.
Nie RE, Wei J, Zhang SK, Vogler AP, Wu L, Konstantinov AS, et al. Diversification of mitogenomes in three sympatric Altica flea beetles (Insecta, Chrysomelidae). Zool Scr. 2019;48:657–666.
Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34:3094–3100.
Kolmogorov M, Yuan J, Lin Y, Pevzner PA. Assembly of long, error-prone reads using repeat graphs. Nature Biotechnol. 2019;37:540.
Pryszcz LP, Gabaldón T. Redundans: an assembly pipeline for highly heterozygous genomes. Nucleic Acids Res. 2016;44:e113.
Boetzer M, Henkel CV, Jansen HJ, Butler D, Pirovano W. Scaffolding pre-assembled contigs using. Bioinformatics. 2011;27:578–579.
Boetzer M, Pirovano W. SSPACE-LongRead: scaffolding bacterial draft genomes using long read sequence information. BMC Bioinformatics. 2014;15:211.
Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.
Waterhouse RM, Seppey M, Simao FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35:543–548.
Xu Z, Wang H. LTR_FINDER: an efficient tool for the prediction of full-length LTR retrotransposons. Nucleic Acids Res. 2007;35:W265–W268.
Price AL, Jones NC, Pevzner PA. De novo identification of repeat families in large genomes. Bioinformatics. 2005;21:i351–i358.
Edgar RC, Myers EW. PILER: identification and classification of genomic repeats. Bioinformatics. 2005;21(Suppl. 1):152–158.
Hoede C, Arnoux S, Moisset M, Chaumier T, Inizan O, Jamilloux V, et al. PASTEC: an automatic transposable element classification tool. PLoS One. 2014;9:e91929.
Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J. Repbase Update, a database of eukaryotic repetitive elements. Cytogenet. Genome Res. 2005;110:462–467.
Tarailo-Graovac M, Chen NS. Using RepeatMasker to identify repetitive elements in genomic sequences. Current Protocols in Bioinformatics. 2009;4.10.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–842.
Burge C, Karlin S. Prediction of complete gene structures in human genomic DNA. J Mol Biol. 1997;268:78–
Stanke M, Waack S. Gene prediction with a hidden Markov model and a new intron submodel. Bioinformatics. 2003;19 (Suppl. 2):215–225.
Majoros WH, Pertea M, Salzberg SL. TigrScan and GlimmerHMM: two open source ab initio eukaryotic gene-finders. Bioinformatics. 2004;20:2878–2879.
Blanco E, Parra G, Guigo R. Using geneid to identify genes. Current Protocols in Bioinformatics. 2007; 4.3.
Korf I. Gene finding in novel genomes. BMC Bioinform. 2004;5:59.
Keilwagen J, Wenk M, Erickson JL, Schattat MH, Grau J, Hartung F. Using intron position conservation for homology-based gene prediction. Nucleic Acids Res. 2016;44:e89.
Keilwagen J, Hartung F, Paulini M, Twardziok SO, Grau J. Combining RNA-seq data and homology-based gene prediction for plants, animals and fungi. BMC Bioinform. 2018;19:189.
Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11:1650.
Campbell MA, Haas BJ, Hamilton JP, Mount SM, Buell CR. Comprehensive analysis of alternative splicing in rice and comparative analyses with Arabidopsis. BMC Genomics. 2006;7:327.
Haas BJ, Salzberg SL, Zhu W, Pertea M, Allen JE, Orvis J, et Automated eukaryotic gene structure annotation using EVidenceModeler and the program to assemble spliced alignments. Genome Biol. 2008;9:R7.
Griffiths-Jones S, Moxon S, Marshall M, Khanna A, Eddy SR, Bateman A. Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res 2005;33(Database issue):121–124.
Griffiths-Jones S, Grocock RJ, van Dongen S, Bateman A, Enright AJ. miRBase: microRNA sequences, targets and gene nomenclature. Nucleic Acids Res. 2006;34(Database issue):140–144.
Nawrocki EP, Eddy SR. Infernal 1.1: 100-fold faster RNA homology searches. Bioinformatics. 2013;29:2933–2935.
Lowe TM, Eddy SR. tRNAscan-SE: a programfor improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res. 1997;25:955–964.
She R, Chu SC, Uyar B, Wang J, Wang K, Chen NS. genBlastG: using BLAST searches to build homologous gene models. Bioinformatics. 2011;27:2141–2143.
Birney E, Clamp M, Durbin R. GeneWise and Genomewise. Genome Res. 2004;14:988–995.
Marchler-Bauer A, Lu S, Anderson JB, Chitsaz F, Derbyshire MK, DeWeese-Scott C, et al. CDD: a Conserved Domain Database for the functional annotation of proteins. Nucleic Acids Res. 2011;39(Database issue):D225–229.
Boeckmann B, Bairoch A, Apweiler R, Blatter MC, Estreicher A, Gasteiger E, et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res. 2003;31:365–370.
Tatusov RL, Natale DA, Garkavtsev IV, Tatusova TA, Shankavaram UT, Rao BS, et al. The COG database: new developments in phylogenetic classification of proteins from complete genomes. Nucleic Acids Res. 2001;29:22–28.
Kanehisa M, Goto S.KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 2000;28:27–30.
Altschul SF, Gish W, Miller W, Myers EW, Lipman Basic local alignment search tool. J Mol Biol. 1990;215:403–410.
Dimmer EC, Huntley RP, Alam-Faruque Y, Sawford T, O'Donovan C, Martin MJ, et al. The UniProt-GO Annotation database in 2011. Nucleic Acids Res. 2012;40:D565–D570.
Conesa A, Gotz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–3676.
Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59:307–321.
Yang ZH, Rannala B. Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds. Mol Biol Evol. 2006;23:212–226.
Rannala B, Yang ZH. Inferring speciation times under an episodic molecular clock. Syst Biol. 2007;56:453–466.
Han MV, Thomas GWC, Lugo-Martinez J, Hahn MW. Estimating gene gain and loss rates in the presence of error in genome assembly and annotation using CAFE 3. Mol Biol Evol. 2013;30:1987–1997.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2017;bbx108.
Darriba D, Taboada GL, Doallo R, et al. ProtTest 3: fast selection of best-fit models of protein evolution. Bioinformatics. 2011;27:1164–1165.
Stamatakis A. RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;22:2688–2690.
Miller MA, Pfeiffer W, Schwartz T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In Gateway Computing Environments Workshop (GCE). New Orleans, LA: IEEE. 2010.
Rambaut A. FigTree v1.3.1: Tree figure drawing tool. 2009. Available:http://tree.bio.ed.ac.uk/software/figtree/.
Wang YQ, Song FH, Zhu JW, Zhang SS, Yang YD, Chen TT, et al. GSA: Genome sequence archive. Genom Proteom Bioinf. 2017;15:14–18.
National Genomics Data Center Members and Partners. Database Resources of the National Genomics Data Center in 2020. Nucleic Acids Res. 2020;48:D24–D33.

Table 1 BUSCO results showing completeness of the Altica viridicyanea genome assembly and annotation.

	Insect		Endopterygota
	Counts	Percentage	Counts	Percentage
Complete BUSCOs	1,635	98.6%	2340	95.8%
Complete and single-copy BUSCOs	1,540	92.9%	2211	90.5%
Complete and duplicated BUSCOs	95	5.7%	129	5.3%
Fragmented BUSCOs	7	0.4%	52	2.1%
Missing BUSCOs	16	1.0%	50	2.1%
Total BUSCO groups searched	1,658		2,442

Table 2 Composition of repetitive sequences in the Altica viridicyanea genome assembly.

Repeat type		Number of elements	Length (bp)	Rate (%)
Retrotransposons (transposable element class I)		1,033,903	356,902,348	41.27
	DIRS	12,458	7,468,113	0.86
	LINE	3,17,300	97,559,643	11.28
	LTR/uncertain	45,604	28,406,497	3.28
	LTR/Copia	15,362	7,752,670	0.9
	LTR/Gypsy	143,535	83,332,210	9.64
	LTR or DIRS	91	14,310	0
	PLE or LARD	485,362	155,872,484	18.02
	SINE	78	80,868	0.01
	TRIM	13,553	6,401,760	0.74
	Unknown	560	58,144	0.01
DNA transposons (transposable element class II)		803,681	226,943,053	26.24
	TIR	550,895	150,574,911	17.41
	MITE	10	635	0
	Crypton	145,294	43,476,772	5.03
	Helitron	29,008	8,403,268	0.97
	Maverick	51,701	31,213,827	3.61
	Unknown	26,773	5,423,361	0.63
SSR		583	165,863	0.02
Unknown		77,813	23,953,352	2.77
Total		1,848,679	544,002,544	62.91

Notes: DIRS, dictyostelium intermediate repeat sequence; LINE, long interspersed nuclear element; LTR, long terminal repeat; PLE, penelope-like elements; SINE, short interspersed nuclear element; LARD, large retrotransposon derivative elements; TIR, terminal inverted repeat.

Table 3 Statistics of gene prediction of Altica viridicyanea.

Method	Software	Gene number
ab initio	Genscan v1.1.0	15,170
	Augustus v2.4	31,813
	GlimmerHMM v3.0.4	58,872
	GeneID v1.4	14,970
	SNAP v2006-07-28	72,679
homology-based	GeMoMa v1.3.1 Drosophila melanogaster Tribolium castaneum Dendroctonus ponderosae Anoplophora glabripennis	9,310
		14,234
		13,066
		18,274
transcriptome-based	TransDecoder v2.0	73,432
	GeneMarkS-T v5.1	29,645
	PASA v2.0.2	23,200
Integration	EVidenceModeler v1.1.1	17,730

Download PDF

Journal Publication

published 07 Apr, 2021

Read the published version in BMC Genomics →

Editorial decision: Major revision
07 Jan, 2021
Review #2 received at journal
28 Dec, 2020
Review #1 received at journal
20 Dec, 2020
Reviewer #3 agreed at journal
09 Dec, 2020
Reviewer #2 agreed at journal
06 Dec, 2020
Reviewer #1 agreed at journal
06 Dec, 2020
Editor assigned by journal
03 Dec, 2020
Reviewers invited by journal
03 Dec, 2020
Submission checks completed at journal
03 Dec, 2020
Editor invited by journal
03 Dec, 2020

You are reading this latest preprint version

The draft genome of the specialist flea beetle Altica viridicyanea (Coleoptera: Chrysomelidae)

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Results And Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Tables

Table 2 Composition of repetitive sequences in the Altica viridicyanea genome assembly.

Supplementary Files

Status:

Journal Publication

Version 2