The genome and stage‐specific transcriptomes of the carrot weevil, Listronotus oregonensis, reveal adaptive mechanisms for host specialisation and symbiotic interactions

TL;DR

This study provides the genome and transcriptomes of the carrot weevil, revealing genetic adaptations for its narrow host range and symbiotic relationships.

Contribution

The study identifies unique genomic features, including horizontal gene transfer and gene expansions, linked to host specialization and symbiosis in L. oregonensis.

Findings

L. oregonensis has 14,637 annotated genes with distinct expression patterns across developmental stages.

The species shows reduced chemoreceptor genes and expanded detoxification gene families compared to related species.

Evidence of horizontal gene transfer includes cell-wall degrading enzymes and a toxin gene likely from prokaryotes.

Abstract

Throughout their evolution, insects have become specialised to occupy diverse ecological niches. The carrot weevil, Listronotus oregonensis, is an important agricultural pest that exhibits a very specific host range. In this study, we characterised the genome and transcriptomes of each developmental stage of L. oregonensis and its Wolbachia endosymbiont to gain deeper knowledge of the genetic determinants controlling its biology. We annotated 14,637 genes and showed expression profiles across the developmental stages. We also compared orthologous genes between L. oregonensis and nine other species, with particular focus on chemoreceptors and detoxification genes. We identified 24 distinct odorant‐binding protein genes and 41 genes for receptors involved in stimulus perception, relatively low numbers compared with other species, which would be consistent with a narrow host range. In contrast, we found a high number of detoxification genes, with significant expansion of certain gene families. Among the annotated genes, 46 were putatively acquired through horizontal gene transfer, with 17 showing strong evidence for this, including several cell‐wall degrading enzymes. The phylogeny of a cytolethal distending toxin gene also suggests an initial transfer from a prokaryotic source and vertical dissemination in members of Curculionidae through recent evolution. The presence of the endosymbiotic bacterium Wolbachia (supergroup A) was confirmed in all tested L. oregonensis individuals from several regions in northeastern North America and showed very little diversity. This study enhances our understanding of the genomic, functional, and evolutionary aspects of a significant agricultural pest and makes important and useful databases available to the scientific community.

Comprehensive genome and stage‐specific transcriptomes reveal 14,637 genes in Listronotus oregonensis, advancing genetic insights into a key agricultural pest.Compared with related species, L. oregonensis exhibits reduced chemoreceptor gene families but expanded detoxification gene clusters, reflecting its narrow host range and adaptability.Evidence of horizontal gene transfer includes cell wall‐degrading enzymes and a cytolethal distending toxin.

Comprehensive genome and stage‐specific transcriptomes reveal 14,637 genes in Listronotus oregonensis, advancing genetic insights into a key agricultural pest.

Compared with related species, L. oregonensis exhibits reduced chemoreceptor gene families but expanded detoxification gene clusters, reflecting its narrow host range and adaptability.

Evidence of horizontal gene transfer includes cell wall‐degrading enzymes and a cytolethal distending toxin.

Figures

Life stage‐specific gene expression in Listronotus oregonensis and its endosymbiont Wolbachia. (a) Number of expressed genes in different stages and, in parentheses, the number of genes unique to these stages. (b) Cell localisation of expressed proteins.

Principal component analysis of the gene expression profile in each developmental stage of Listronotus oregonensis. L1 to L4 are the larval stages, male adults are denoted by ♂ and females by ♀.

Phylogenomic tree obtained from the OrthoFinder comparison of all high‐confidence orthologous genes. Insect families are indicated on the right side. Numerical values on tree branches show OrthoFinder support values based on the proportion of species trees derived from single‐locus gene trees supporting each bipartition.

Annotation of the Wolbachia genome and its prophage region. (a) Gene positions and organisation of prophage insertions (middle circle). Colours identify prophage regions that are complete (green) or incomplete (orange). The predicted tail and head regions are indicated in the inner circle. (b) Structural organisation of genes in the complete prophage region 6, with the different components shown in different colours.

Tables

TABLE 1. Genome completeness assessment using single‐copy orthologue genes (BUSCO) for Listronotus oregonensis and Wolbachia assemblies.

BUSCOs	L. oregonensis	Wolbachia
Complete	839 (82.8%)	330 (90.7%)
Complete and single‐copy	826 (81.5%)	321 (88.2%)
Complete and duplicated	13 (1.3%)	9 (2.5%)
Fragmented	75 (7.4%)	24 (6.6%)
Missing	99 (9.8%)	10 (2.7%)
Total searched	1013	364
odb10	Arthropoda	Rickettsiales

TABLE 2. Annotation of transcripts in Listronotus oregonensis and Wolbachia.

Category	L. oregonensis	Wolbachia
Total	21,010 | 100%	1340 | 100%
No hits	2781 | 13.2%	30 | 2.2%
Blast hits	18,229 | 86.7%	1310 | 97.8%
GO	8065 | 38.4%	43 | 3.2%
KEGG	7277 | 34.6%	702 | 52.4%
COG	12,533 | 59.6%	981 | 73.2%
pFAM	13,528 | 64.4%	963 | 71.9%
SignalP	1528 | 7.3%	45 | 3.4%
DeepTMHMM	1220 | 5.8%	240 | 17.9%

TABLE 3. Identity and annotation of Listronotus oregonensis genes putatively acquired by horizontal gene transfer.

Gene ID	Annotation	Origin	E‐value
LiOr723.t1	Cytolethal distending toxin	Bacteria	9.58e‐36
LiOr723.t2	Cytolethal distending toxin	Bacteria	7.09e‐42
LiOr724.t1	Cytolethal distending toxin	Bacteria	3.48e‐15
LiOr1549.t1	ORF2	Viruses	2.07e‐170
LiOr2322.t1	Pectinesterase	Bacteria	1.59e‐82
LiOr4266.t1	Rhamnogalacturonate lyase	Eukaryota@Metazoa	1.64e‐145
LiOr4267.t1	Rhamnogalacturonate lyase	Eukaryota@Metazoa	6.68e‐150
LiOr4528.t1	O‐acetyl‐ADP‐ribose deacetylase	Bacteria	1.93e‐22
LiOr4543.t1	Pectinesterase	Bacteria	2.51e‐86
LiOr5391.t1	DNA primase	Eukaryota@Fungi	1.96e‐47
LiOr5633.t1	Rhamnogalacturonate lyase	Eukaryota@Metazoa	2.37e‐145
LiOr5633.t2	Rhamnogalacturonate lyase	Eukaryota@Metazoa	3.91e‐145
LiOr8069.t1	Cytolethal distending toxin	Bacteria	3.13e‐48
LiOr12562.t1	Pectinesterase	Bacteria	2.64e‐28
LiOr13508.t1	Acyl‐CoA‐binding protein	Bacteria	4.55e‐18
LiOr16754.t1	Cytolethal distending toxin	Bacteria	9.58e‐13
LiOr17006.t1	Pectinesterase	Bacteria	2.04e‐87
LiOr17336.t1	Pectinesterase	Bacteria	7.07e‐86
LiOr18761.t1	Rhamnogalacturonate lyase	Eukaryota@Metazoa	5.67e‐98

TABLE 4. Comparison of the number of odorant‐binding proteins (OBPs), odorant receptors (ORs), gustatory receptors (GRs), ionotropic receptors (IRs), cytochrome P450s (CYPs), carboxyl/cholinesterases (CCEs), and glutathione S‐transferases (GSTs) identified in different insect species.

	OBP	OR	GR	IR	CYP	CCE	GST
Anthonomus grandis	29	12	13	21	37	11	4
Callosobruchus chinensis	19	10	18	22	41	14	4
Callosobruchus maculatus	18	7	15	24	37	13	4
Ceutorhynchus obstrictus	23	9	10	18	40	11	4
Dendroctonus ponderosae	28	16	15	23	33	12	4
Listronotus bonariensis	25	6	15	21	49	18	3
Listronotus oregonensis	24	13	13	15	36	14	4
Rhynchophorus ferrugineus	33	7	15	20	41	13	4
Sitophilus oryzae	36	19	20	21	47	15	4
Tribolium castaneum	51	21	21	20	80	13	4

Full text

INTRODUCTION

Insects are among the most diverse taxonomic groups on Earth, showcasing remarkable adaptability and evolutionary plasticity (Grimaldi & Engel, 2005). This adaptability is driven not only by genetic diversity and natural selection but also by interactions with other organisms, such as pathogens and symbionts, which can introduce new genetic material through processes like lateral gene transfer or influence host biology via mutualistic relationships (Werren et al., 2008). These processes have allowed insects to acquire extraordinary capabilities, from surviving extreme environmental conditions to exploiting novel ecological niches (Storey & Storey, 2012).

Horizontal gene transfer (HGT) refers to the non‐sexual transmission of genetic material between different species (Kidwell, 1993). This process has played a significant role in the evolution of all taxonomic groups, including eukaryotes (Keeling & Palmer, 2008). For example, nematodes have gained the ability to parasitise and feed on plants by acquiring genes from fungi and bacteria that degrade plant cell walls (Danchin et al., 2010). HGT is also believed to facilitate insects' adaptation to herbivory by enabling them to use genes from prokaryotes to detoxify plant metabolites (Wybouw et al., 2016). A recent study demonstrated this phenomenon in the whitefly Bemisia tabaci, which utilises an HGT‐acquired phenolic glucoside malonyltransferase gene to degrade phenolic glucoside compounds in host plants, giving it a significant feeding advantage over species of the related genus Trialeurodes on the same hosts (Xia et al., 2021). Therefore, HGT provides an evolutionary advantage to insects by aiding in their adaptation to various hosts or conditions, thereby enhancing insect fitness and overall evolutionary flexibility.

Insects interact with a diversity of organisms, sometimes establishing an intimate symbiosis, with positive (mutualism) or negative (parasitism) outcomes. The intracellular bacteria Wolbachia is probably the best known endosymbiont, as it can drastically alter its host's biology through reproductive manipulations (Werren et al., 2008). These kinds of interactions can provide many benefits to the host, such as allowing the acquisition of essential nutrients or protection against parasitic infections. The association with Wolbachia has evolved into an obligate one not only in various insect species, but also in other non‐related groups such as filarial nematodes, where the host cannot complete its life cycle without the presence of the symbiont (Slatko et al., 2010). The comparison of different Wolbachia strains by mining genome sequences from different groups did not reveal the presence of strict cophylogeny, which indicates that horizontal transfers occur frequently between groups (host switching) and that Wolbachia infection could be a dynamic driver of evolution in the host (Vancaester & Blaxter, 2023). These bacterial endosymbionts have been shown to be a strong driver of adaptation in weevils and a major factor in the success of agricultural pests, therefore representing interesting targets for novel control strategies (Morera‐Margarit et al., 2021). The presence of the bacteriophage WO has also been confirmed in several strains of Wolbachia, and has been shown to influence the effect of the bacterium on its host, adding another layer of complexity (Masui et al., 2000).

The carrot weevil, Listronotus oregonensis (Leconte) (Coleoptera: Curculionidae), a native species of eastern North America, poses a significant threat to crops in the Apiaceae family, including carrot, parsley and celery (Boivin, 2013). The larvae of this pest burrow into plant roots, causing direct damage that can potentially lead to substantial yield losses when no control measures are implemented (Justus & Long, 2019). Current management strategies rely heavily on insecticide applications. However, the emergence of a second generation has led to increased treatment frequency, potentially fostering the development of insecticide resistance (Gagnon & Bourgeois, 2024; Telfer et al., 2019). This shift in pest phenology and local adaptations is already resulting in higher levels of damage, necessitating the exploration of alternative control methods (Gagnon & Blatt, 2024).

Several biological control methods targeting the carrot weevil have been attempted so far with mixed success, involving different species of parasitoids (Huber et al., 1997) or entomopathogenic nematodes (Bélair & Boivin, 1995). The parasitic nematode Bradynema listronoti (Zeng et al., 2007) has demonstrated good control potential by sterilising female weevils (Gagnon et al., 2019). This castration appears to be the result of hormonal changes and resource reallocation, but it remains unclear which genetic processes are involved. Trials with various exotic parasitoids have also been attempted without successful establishment. Interestingly, the importance of odour and pheromone production and recognition for host location and infection in these trials has been identified as a key factor (Cournoyer & Boivin, 2004). Similarly, the carrot weevil must inevitably recognise the specific odours of host plants to locate them (Gagnon et al., 2021). Yet, genetic comparison has not shown any structure between carrot weevil populations developing on different host plants (Bessette et al., 2022). In all these cases, understanding the genetic basis behind host range and biotic interactions will be crucial to continue developing new sustainable management practices.

To gain deeper insights into the carrot weevil's biology, we combined genome annotation with comparative genomics and transcriptome analysis. By comparing the carrot weevil genome with those of related species, we inferred phylogenetic relationships and identified shared genes that may underpin its ecological and evolutionary traits. Stage‐specific transcriptomes were analysed to determine gene expression patterns across its life stages, shedding light on developmental processes and potential targets for pest management. Furthermore, we characterised the strain diversity of Wolbachia in L. oregonensis and investigated possible HGT events. These findings provide a foundation for understanding the genetic and symbiotic factors driving the carrot weevil's pest status and support the development of alternative pest management strategies. These genomic and transcriptomic resources will be of great importance to support future studies on the genetic determinants driving weevil host range and biotic interactions.

MATERIALS AND METHODS

Biological material

The population of carrot weevils (L. oregonensis) used in this study comes from a rearing colony at the Saint‐Jean‐sur‐Richelieu Research and Development Centre of Agriculture and Agri‐Food Canada. The initial specimens were collected from a carrot field in Sainte‐Clotilde‐de‐Chateauguay (45.162230, −73.673545) and multiplied on carrots for several generations. The different stages (eggs, four larval stages [L1 to L4], pupae, adult) were obtained throughout its life cycle. Field specimens were also collected from three fields of each of the Canadian provinces of Quebec, Nova Scotia, and Ontario and the U.S. state of Ohio. These field samples were used to determine the prevalence of Wolbachia and genetic diversity.

The endosymbiotic bacterium Wolbachia was found to be naturally present in this weevil population. To assemble the genome of this endosymbiont, a pool of 25 weevils from the rearing colony was used. Weevils were placed in a killing jar with ethanol and surface‐sterilised by dipping them in 70% ethanol. They were transferred to sterile water, and the head and the digestive systems were then removed under a stereomicroscope to limit the risk of bacterial contamination. The weevils were then bulked into a 10‐mL tube containing 5 mL of cold Eagle's Medium (B9638 – Sigma‐Aldrich, St. Louis, Missouri, United States), and vortexed for 10 min with 50% v/v of 3 mm sterile borosilicate glass beads. Following vortexing, the resulting cell lysate was transferred to a clean 10‐mL tube and centrifuged at 300g for 5 min to precipitate the intact cells and weevil debris. The supernatant was further centrifuged at 3000g for 10 min to remove cell debris. This cleared cell lysate was filtered through a 0.45‐micron syringe filter, and the bacterial cells were then enriched by centrifuging the sample at 16,500g for 10 min. The resulting Wolbachia‐enriched pellet was transferred to a sterile microtube and stored at −80°C to await downstream DNA extraction.

DNA/RNA extraction, library preparation and sequencing

For L. oregonensis, total RNA was extracted from 50 eggs; 35 L1, 10 L2, one L3, and one L4 individuals; one pupa; one adult male; and one adult female using the RNeasy Plus Mini Kit (Qiagen). The quality of the RNA extracts was verified with a Bioanalyzer (Agilent), and only those with a RIN (RNA Integrity Number) value higher than 8 were sent for sequencing. Three replicates of each life stage were prepared. All libraries and sequencing were carried out at the McGill University and Génome Québec Innovation Centre using the NovaSeq6000 platform (2 × 150 bp) (Illumina).

For the endoparasitic Wolbachia, gDNA was amplified from the bacterial cell pellet isolated from mature weevils using the REPLI‐g Advanced DNA Single Cell Kit following the manufacturer's recommended protocol (Qiagen). A library was prepared from 1.5 μg of amplified DNA using the nanopore long‐read compatible Ligation Sequencing Kit V14 (SQK‐LSK114) (Oxford Nanopore Technologies, Oxford, United Kingdom). The sequencing was performed using MinION R10.4 flow cells (Oxford Nanopore Technologies) at the Saint‐Jean‐sur‐Richelieu Research and Development Centre.

Genome assembly and gene prediction

The draft genome of L. oregonensis was generated using the hybrid genome assembler HASLR combining Oxford Nanopore and Illumina reads (details in Bessette et al., 2022). For the present study, this 1.3‐Gb assembly was masked using RepeatMasker (Smith et al., 2015) and annotated using BRAKER3 (Gabriel et al., 2024), guided using both RNASeq data, collected from all life stages, as well as protein evidence recovered from the OrthoBD 11 Arthropoda database (Kuznetsov et al., 2023). The predicted gene features were then filtered to isolate those having RNASeq support (>50 supporting paired‐end sequences).

The Wolbachia genome was assembled from long reads using Flye. Briefly, R10.4‐generated long reads were basecalled and trimmed with high‐accuracy basecalling models. These reads were then inputted into Flye and assembled using the meta‐genome parameters. The circular sequence corresponding to the Wolbachia genome was then extracted and submitted to the NCBI PGAP for gene annotation (Haft et al., 2018; Li et al., 2021; Tatusova et al., 2016). Completeness was assessed using BUSCO (Seppey et al., 2019) and the rickettsiales_odb10 list of reference genes. The phylogenetic position of the isolate was tested by comparing the genome sequence with those of representatives of supergroups A and B (Ramírez‐Puebla et al., 2015), using the autoMLST pipeline (Alanjary et al., 2019). After uploading our Wolbachia genome sequence, the pipeline identified the nearest reference organisms from NCBI RefSeq before selecting the best single‐copy genes to run the phylogenetic analysis.

Functional annotation

The functional annotation of the predicted protein sequences was carried out using the eggNOG genome‐wide functional annotation web‐server under default parameters (http://eggnog-mapper.embl.de/; Cantalapiedra et al., 2021; Huerta‐Cepas et al., 2019). The presence of signal peptide signatures, an indicator of extracellular protein secretion, was detected using signalP 6.0 “Eukarya” settings for L. oregonensis and signalP 6.0 “Other” settings for Wolbachia (Teufel et al., 2022). In addition, the subcellular localisation of proteins in L. oregonensis was predicted using the eukaryote‐specific tool DeepLoc2 under default settings (Thumuluri et al., 2022). Likewise, the gram‐negative prokaryotic‐specific predictor PSORTb was used to predict the subcellular localisation of Wolbachia proteins (v3.0.2 – https://www.psort.org/psortb/; Yu et al., 2010). The online PHAge Search Tool Enhanced Release (PHASTER) web‐server was also utilised to identify any phage segments integrated in the Wolbachia genome (Arndt et al., 2016; Zhou et al., 2011).

Life stage‐specific gene expression

The life stage‐specific transcriptomes of L. oregonensis were analysed using RNA‐seq data generated from each individual life stage. Briefly, RNA reads, generated as described above, were quality‐trimmed to a minimum phred‐score of 30 using Fastp (Chen et al., 2018). Trimmed paired‐end read sets were then mapped against the genome sequence using the splice‐aware aligner Hisat2 (Kim et al., 2019). Mapping was performed separately for each replicate of a given life stage and for each sampled life stage (Supplementary Table S1). Read counts were extracted, per species, from the mapping files using the featureCount tool (Liao et al., 2014), with the parameters –countReadPairs, −t exon, and ‐g gene_id set. Gene transfer format (GTF) annotations, generated from the BRAKER3 pipeline, were supplied and reads overlapping these features were extracted. A principal component analysis was used to assess the expressional landscape across the different life stages. Here, we define a stage‐specific gene as a gene whose expression was higher than the minimum 50 reads threshold in only one of the life stages.

Horizontal gene transfer

To evaluate the prevalence of HGT events in the L. oregonensis genome, the Alienness versus Predictor (AvP) automated pipeline was used (v1.0.10– Koutsovoulos et al., 2022). This pipeline, run under default settings, uses both alienness scores (AI and AHS) as well as phylogenetic tree inference to generate a list of robust HGT candidates. As input, the predicted amino acid sequences were blasted against the non‐redundant (nr) database using the Diamond‐BLASTp algorithm, using an E‐value cutoff score of 1E‐05 (Buchfink et al., 2021). To validate that these candidates were indeed genes acquired through HGT and not contamination, the downstream AvP hgt_local_score analysis was performed using the BRAKER3 generated gff3 as input. Each putative HGT gene was evaluated for the presence of introns. Genomic and transcriptomic PCR amplicon validation were then performed for two high‐confidence HGT candidates using the following primers: cdtB_F (5′‐GGGAGGAAACAGAGTCAA‐3′) and cdtB_R (5′‐GGGGTAAATCGTAGCTAT‐3′) to validate the cytolethal distending toxin B gene, and PME_F (5′‐CAAGGATGGAATGTGGAA‐3′) and PME_R (5′‐TGGAACCAGAACCGTAAA‐3′) to target the pectinesterase gene.

Comparative genomics

To determine the phylogenetic relationships of L. oregonensis, we retrieved protein sequences from five species of weevils of agricultural importance, including the cabbage seed weevil, Ceutorhynchus obstrictus (subfamily: Ceutorhynchinae) (PRJEB47903), the boll weevil, Anthonomus grandis (subfamily: Curculioninae) (PRJNA767408), the Argentine stem weevil, Listronotus bonariensis (subfamily: Cyclominae) (PRJNA640511), the palm weevil Rhynchophorus ferrugineus (subfamily: Dryophthorinae) (PRJNA599933), and the rice weevil Sitophilus oryzae (subfamily: Dryophthorinae) (PRJNA566109); two related agricultural pests that are not true weevils, comprising the cowpea weevil, Callosobruchus maculatus (Coleoptera: Chrysomelidae) (PRJEB30475) and the adzuki bean weevil, Callosobruchus chinensis (Coleoptera: Chrysomelidae) (PRJEB52276); and two genetic models, the mountain pine beetle, Dendroctonus ponderosae (Coleoptera: Curculionidae: Scolytinae) (PRJNA846874) and the red flour beetle, Tribolium castaneum (Coleoptera: Tenebrionidae) (PRJNA15718) (Supplementary Table S2). We then used OrthoFinder v2.5.4 (Emms & Kelly, 2015; Emms & Kelly, 2019) with the default parameters to define orthogroups using Diamond v0.9.14 and keeping only the longest protein isoform of each gene. Phylogenetic trees were computed using the tree root inference module in OrthoFinder with multiple sequence alignments (Emms & Kelly, 2017).

Specific gene families associated with odour perception, chemoreceptors and detoxification were analysed in greater detail. Individual protein family databases were created by downloading all the sequences from UniProtKB with Insecta‐linked taxonomy (taxonomy_id: 50557) for the following terms: “Odorant‐binding protein,” “Odorant receptor,” “Gustatory receptors,” “Ionotropic receptor,” “Cytochrome P450,” “Carboxyl/Cholinesterase,” “Acetylcholinesterase” and “Glutathione S‐transferase.” For each species, the proteomes were screened against these curated databases using blastp with a reporting E‐value threshold of 1E‐25.

Wolbachia diversity and cytoplasmic incompatibility

The prevalence and genetic diversity of Wolbachia in L. oregonensis was assessed using weevils collected from three fields each in the provinces of Quebec, Nova Scotia and Ontario, and the state of Ohio as well as four individuals from the lab rearing (as above, Quebec origin) and four from a distinct lab rearing initiated with a population from Ontario, for a total of 20 samples. DNA was extracted from whole weevils using the DNeasy Blood & Tissue Kits (Qiagen). Amplification of the variable Wolbachia surface protein (WSP) gene was conducted using the primers wsp81F (5′‐TGGTCCAATAAGTGATGAAGAAAC‐3′) and wsp691R (5′‐AAAAATTAAACGCTACTCCA‐3′) to confirm if an individual was a carrier of Wolbachia or not (Zhou et al., 1998). The 16S rRNA gene segment was amplified using the *Wolbachia‐*specific 16S primers 16S_281F (5′‐CTATAGCTGATCTGAGAGGAT‐3′) and 16S_1372R (5′‐YGCTTCGAGTGAAACCAATTC‐3′) (Wang et al., 2014). Likewise, sequence variations in the cytoplasmic incompatibility gene CifB were investigated in all above‐mentioned individuals by amplifying a 1913‐bp segment using the primers cifB_Fwd (5′‐TCATACCACATCAGAGCA‐3′) and cifB_Rev (5′‐CATAAACCACAATTCCAACC‐3′). Both 16S and CifB amplicons were sequenced using the Genome Quebec Sanger sequencing platforms.

RESULTS

Genome assemblies

In this experiment, the draft genome of the carrot weevil (PRJNA726155) was annotated using RNASeq and protein data (NCBI BioProject PRJNA1019676 (Ste‐Croix et al., 2023b), individual accession numbers in Supplementary Table S1). Gene expression was monitored in all major developmental stages using a total of 363 million sequencing reads. This resulted in the annotation of 19,111 genes, 14,637 having expression (>50 supporting reads), encoding some 21,010 transcripts (Ste‐Croix et al., 2025). Genome completeness was assessed using the Arthropoda (n = 1013) BUSCO gene reference sets. The results indicated 82.8% genome completeness, with few duplications of BUSCO genes (1.3%) (Table 1). Likewise, the genome of the endosymbiotic bacteria Wolbachia was assembled in a single, 1,567,375‐nucleotide circular sequence using a combination of 2,517,677 enriched long reads of Wolbachia (Ste‐Croix et al., 2023a). Completeness analysis using Rickettsiales BUSCO genes yielded a score of 90.7%, with low levels of BUSCO gene duplications (2.5%) (Table 1). Automated gene annotation using the Prokaryotic Genome Annotation Pipeline (PGAP) of the National Center for Biotechnology Information (NCBI) highlighted 1340 protein‐coding genes (Supplementary Table S3). A whole‐genome phylogenetic comparison with other Wolbachia strains indicated that this isolate belongs to the Wolbachia supergroup A (Supplementary Figure S1).

Functional annotation

Among the predicted genes, 86.7% and 97.8% generated significant BLAST hits when compared with the nr database for L. oregonensis and Wolbachia, respectively (Table 2). These genes were further annotated according to their biological process, molecular function, and/or cellular component (Table 2). Finally, the cellular location of the expression of each gene was predicted for both organisms (Figure 1b).

Life stage‐specific gene expression in Listronotus oregonensis and its endosymbiont Wolbachia. (a) Number of expressed genes in different stages and, in parentheses, the number of genes unique to these stages. (b) Cell localisation of expressed proteins.

Life‐stage‐specific gene expression

An analysis of the different developmental stages of L. oregonensis showed consistent patterns of gene expression and stage‐specific expression across all major stages. The adult stages had the highest number of expressed genes (>8300), while expression in all the other stages was consistently less, ranging from 6680 (L3) to 7727 (L2) (Figure 1a; Supplementary Table S4). Only a small number of genes were uniquely associated with a single developmental stage. For instance, several cuticle proteins were exclusively and highly expressed during the pupal stage, whereas the vitellogenin protein LiOr.4819.t1 was detected solely in adult females, whilst testis‐specific genes were observed exclusively within adult males.

Although the total number of expressed genes appeared to be similar in the different developmental stages, a closer examination of gene expression patterns revealed significant expression clusters corresponding to critical developmental periods. For instance, three distinct clusters were found in L. oregonensis, representing the larval stage, the adult stage, and the egg/pupal stages (Figure 2).

Principal component analysis of the gene expression profile in each developmental stage of Listronotus oregonensis. L1 to L4 are the larval stages, male adults are denoted by ♂ and females by ♀.

Horizontal gene transfer

The AvP pipeline was employed to examine HGT occurrences. A total of 46 genes were identified as having strong evidence for HGT events (Supplementary Table S5). To avoid calling contaminant sequences as HGT, all the genes without introns were removed; 17 genes remained. The identified HGT events are attributed to multiple sources, including 11 genes of bacterial origin, one of viral origin, one of fungal origin, and three putative ancient HGT events with strong phylogenetic support among other insects. Notably, the high‐confidence candidates include a cytolethal distending toxin of bacterial origin and cell wall‐degrading enzymes (CWDEs), such as a bacterial pectinesterase (Table 3). When compared with two databases (NCBI‐nr and UniProtKB), the amino acid sequence identified as cytolethal distending toxin showed homology with only one other eukaryote, the cabbage seed weevil (Ceutorhynchus assimilis, 62.4% identity), a closely related species, but with numerous bacterial genera often associated with gastrointestinal pathologies in humans and animals such as Escherichia, Salmonella, Helicobacter and Campylobacter (40.7–50% identity) (Supplementary Figure S2). We also retrieved a homologous sequence in the sister species L. bonariensis with 79.3% identity, confirming a lateral transfer from prokaryotes, followed by vertical transmission in this group of weevils (Supplementary Figure S2). PCR validation of the cdtB HGT candidate LiOr8069 yielded the expected fragment sizes of 5496 bp from genomic DNA (gDNA) isolated from legs and 404 bp from cDNA. Similarly, LiOr17006, encoding a pectinesterase predicted to have been acquired from a Pseudomonadota via HGT, produced the expected 3619 bp fragment from gDNA and a 214 bp fragment from cDNA.

Comparative genomics

The proteome of L. oregonensis was compared with those of various insects in order to identify shared and unique genes, highlighting similarities and differences in their proteome. The phylogenetic tree of insect species confirmed the close relationship between L. oregonensis and the species from the same genus, L. bonariensis, and its position inside a monophyletic group containing all the members of the Curculionidae family (Figure 3). The two “false weevils” from the leaf beetle family (Chrysomelidae), C. maculatus and C. chinensis, were clustered together in a distinct clade, while the model organism T. castaneum (Tenebrionidae family) was positioned as an outgroup. An OrthoFinder analysis of the proteomes of the 10 insect species was able to classify 217,386 of the 239,901 genes tested (90.6%) in 21,810 orthogroups (Supplementary Table S6). Listronotus oregonensis had 94.7% of its genes assigned to one of the orthogroups, with 398 orthogroups (covering 1851 genes) specific to the species (Supplementary Table S7). Six of these orthogroups were particularly rich in genes (316 genes), with their annotations linked to transposition (e.g., transposases and reverse transcriptases).

Phylogenomic tree obtained from the OrthoFinder comparison of all high‐confidence orthologous genes. Insect families are indicated on the right side. Numerical values on tree branches show OrthoFinder support values based on the proportion of species trees derived from single‐locus gene trees supporting each bipartition.

In order to analyse in more detail the functional capacities of L. oregonensis, its chemoreceptor and detoxification genes were compared more closely with those of the other insects under study (Table 4). A vast repertoire of odorant‐binding proteins was identified (Supplementary Table S8), along with odorant receptors (Supplementary Table S9), gustatory receptors (Supplementary Table S10), ionotropic receptors (Supplementary Table S11), cytochrome P450 monooxygenase (Supplementary Table S12), carboxyl/cholinesterase (Supplementary Table S13), and glutathione S‐transferase (Supplementary Table S14). Of particular interest is one gene (CYP6A14) for cytochrome P450 monooxygenase, which was found to be particularly well represented in L. oregonensis (Figure 4).

Annotation of the Wolbachia genome and its prophage region. (a) Gene positions and organisation of prophage insertions (middle circle). Colours identify prophage regions that are complete (green) or incomplete (orange). The predicted tail and head regions are indicated in the inner circle. (b) Structural organisation of genes in the complete prophage region 6, with the different components shown in different colours.

Wolbachia genomic organisation, diversity and cytoplasmic incompatibility

The analysis of the Wolbachia genome revealed nine potential prophage regions ranging in length from 5.2 Kb to 44.2 Kb, which represent 201 of the 1340 predicted Wolbachia genes (Figure 4). Only one of these regions was found to be complete, containing both attachment (att) sites and essential prophage components (head, tail, and baseplate components) (Figure 4b). All the L. oregonensis individuals tested were positive for the presence of the WSP gene by PCR, confirming the high prevalence of Wolbachia in carrot weevil populations. To assess the potential effect of Wolbachia on cytoplasmic incompatibility in L. oregonensis, we further characterised the CifA and CifB genes. This revealed a solitary CifA‐CifB pair positioned between coordinates 284,032 and 288,973, with CifA positioned upstream of CifB (Figure 4a). CifB was determined to be a peptidase‐type cytoplasmic incompatibility factor, with a papain‐like cysteine peptidase domain identified between positions 854 and 1010 of the protein. These characteristics identify it as a type I CidB (Beckmann et al., 2019; LePage et al., 2017). The sequencing of CifB segments in several individuals from across northeastern North America did not reveal any significant differences in sequences, suggesting a single Wolbachia biotype (Supplementary Figure S3). This finding was also validated through 16S amplicon sequencing (100% homology across all samples, GenBank accession PV500978.1, Mimee et al., 2025).

DISCUSSION

The remarkable ability of insects to adapt to and utilise resources across diverse ecological niches, occasionally specialising to exploit specific resources, is impressive. One such instance of insect resource specialisation can be seen in the carrot weevil, a significant pest affecting several Apiaceae crops. The carrot weevil is found in all production regions, despite the fact that these regions are limited and often geographically scattered. While isolation appears to be a factor driving recent differentiation in carrot weevil populations, the host range in North America has been conserved, with phenotypes remaining very similar (Bessette et al., 2022). In this study, we sought to unravel key aspects of the carrot weevil's biology by focusing on its gene content and stage‐specific expression, as well as symbiotic relationships.

The annotation of L. oregonensis's genome revealed 14,637 genes, with strong evidence of their expression in at least one of the developmental stages. One of the main challenges in accurately predicting genes in L. oregonensis is the high degree of genome fragmentation. This issue largely stems from the difficulties current hybrid assemblers face in resolving the abundant and large “TA” repetitive elements distributed throughout the genome, compounded by its relatively large size of ~1.3 Gb; as highlighted by the BUSCO completeness estimates. To address these difficulties, we used the BRAKER3 gene predictor, which adopts a more conservative approach, reducing the risk of reporting false positives, though at the cost of potentially missing some true bona fide genes. By contrast, the annotation using BRAKER2 (Hoff et al., 2019) on the same dataset predicted ~34,000 genes (available for comparison at https://github.com/Mimee-Lab/Bradynema_Listronotus). While the present annotation represents a considerable improvement over BRAKER2, generating a robust set of genes, there is still room for refinement as assembly and prediction methods continue to advance. A comparison of the amino acid sequences in L. oregonensis with those of other insects revealed that over 91.2% of genes (96% of orthologues) were shared with the other species. On the basis of this set of orthologous genes, the phylogenetic position of L. oregonensis was validated and its gene set compared with those of other weevil species. The phylogenetic analysis conducted in this study situates Listronotus species (subfamily Cyclominae) on a separate branch from Curculioninae and Scolytinae, which is aligned with the findings of studies using alternative phylogenomic methods (Shin et al., 2018).

Chemoreceptors

By relying on sensory neurons situated on their antennae or other external organs, insects are able to detect, taste and locate suitable food sources. These sensory capacities are governed by numerous chemoreceptor genes belonging to a diverse multigene family, encompassing gustatory and odorant receptors (Benton, 2015). In this study, we identified at least 24 genes that encode odorant‐binding proteins (OBPs) in L. oregonensis, along with 41 genes encoding diverse chemosensory receptors, including odorant, gustatory and ionotropic receptors. Certain genes displayed expression patterns that varied across developmental stages. Specifically, the Obp1 gene LiOr16393 and its homologue Obp33 LiOr16394 were expressed exclusively in adults, suggesting a potential role in mate localisation. While pheromones constitute the primary ligands of this OBP in Delia antiqua and Drosophila melanogaster (Rihani et al., 2021), the closest homologue in T. castaneum appears to function as a general odorant‐binding protein, responsive to a broad range of organic volatiles, including 2‐hexanone, β‐ionone, 4,8‐dimethyl‐decanal, 2‐heptenal and 6‐methyl‐5‐hepten‐2‐one (Montino et al., 2021). Many of these compounds are prominent volatiles of Apiaceae host plants, such as carrot, the principal food source of L. oregonensis. This suggests that the corresponding OBP in L. oregonensis may contribute to mate localisation by mediating the detection of food‐associated volatiles from oviposition sites. This interpretation is further supported by its localised expression within the antennae and legs of all three species (Montino et al., 2021; Rihani et al., 2021). In contrast, the Obp14 gene LiOr3915 was almost exclusively expressed within the pupal stages, suggesting a role other than chemosensing.

Compared with OBPs, the chemosensory receptor genes, including IRs, ORs and GRs, exhibited low expression levels across all life stages. Interestingly, L. oregonensis was also predicted to possess substantially fewer receptor genes than the several hundred typically observed in other insects (Robertson, 2019). This pattern aligns with the broader trend across insects, in which receptor gene family size is generally correlated with dietary range: specialists, such as obligate parasites or symbionts restricted to a single host, typically encode relatively few receptors, whereas generalists often exhibit extensive expansions of both GR and IR (Robertson, 2019). Consistent with this, T. castaneum, a generalist feeder on stored grains, possesses a substantially larger receptor repertoire, including nearly 200 more GR genes and 103 IR genes, whereas the comparatively limited receptor set in L. oregonensis likely reflects its narrow host preference for fresh carrots (Gagnon et al., 2021).

Detoxification

The success of an insect hinges not only on its ability to identify a suitable host but also on its capacity to neutralise its host's chemical defences, enabling it to feed successfully. The size and composition of detoxification gene families, such as cytochrome P450 monooxygenase (CYP), carboxyl/cholinesterase (CCE) and glutathione S‐transferase, have been shown to be correlated with food preferences (Pearce et al., 2017). Because these same gene families are also frequently involved in xenobiotics catabolism, it has been hypothesised that their expansion could also favour the emergence of insecticide resistance. In this study, the number of CYP genes was found to be similar in most of the species compared (mean = 44), except for the two model organisms, D. ponderosae, which had fewer (n = 33), and T. castaneum, which had significantly more (n = 80). A similar pattern was observed in L. oregonensis, where the total number of CYP copies identified (n = 179) closely matched that of T. castaneum (n = 221). However, despite overall similarity in gene counts, the CYP6A14 gene was exceptionally abundant (n = 42), with nearly 10 times more copies than in most other insects examined and pointing to a marked, lineage‐specific expansion. The involvement of this family (CYP6) in insecticide resistance has been confirmed in numerous systems and is associated with P450 blooms (lineage‐specific family expansions) (Nauen et al., 2022). In this analysis, 31 copies of CYP6A14 showed consistent background levels of expression across all sampled life stages, with one copy (LiOr13520) exhibiting nearly 10‐fold higher expression. With some populations of L. oregonensis already showing resistance to the insecticide Phosmet, understanding the functional relevance of this expansion and whether copy‐specific increases in expression correlate with elevated resistance to this compound could not only clarify its molecular basis but also inform future management strategies.

CCE gene diversity was also highest in the Listronotus species, highlighting these species' potential for detoxifying plant secondary metabolites. Both species have a developmental stage that resides within plant tissue, making them more exposed to plant defence compounds. Surprisingly, only two regions coding for putative acetylcholinesterase (AChE) genes were identified in L. oregonensis, the lowest number among all species tested (which had four to eight regions). Phosmet, an acetylcholinesterase inhibitor, was until recently among the most widely used insecticides in carrot cultivation and has been associated with suspected cases of resistance (Telfer et al., 2019). In this case, resistance does not appear to require copy number amplification of AChE, but could still arise from a point mutation altering the target site, as observed in other species (Hemingway et al., 2004).

Horizontal gene transfer

In L. oregonensis, a total of 46 genes were identified as probable products of HGT, with 17 of these being robustly supported. Several of these genes, which were actively transcribed in the weevils, were found to code for a plant CWDE pectinesterase of bacterial origin, in line with the already well‐supported hypothesis that the evolution of insect herbivory is linked to HGT events (Wybouw et al., 2016). At least three distinct versions of this enzyme were identified, each displaying differential expression aligned with the developmental requirements of L. oregonensis. For example, the pectinesterase gene LiOr17006 showed markedly elevated expression during actively feeding larval stages, with levels more than 38‐fold higher than those observed across all other stages. In contrast, the cytolethal distending toxin genes (LiOr723, LiOr724, LiOr8069 and LiOr16754) were relatively less expressed, and almost exclusively in the egg and pupal stages. The phylogenetic analysis of these genes supports the hypothesis of a horizontal transfer from a prokaryotic source to an ancestor of Curculionidae and its subsequent vertical transmission. Surprisingly, a gene coding for this toxin would also have been acquired in Drosophila from a bacteriophage (Verster et al., 2019), but does not share significant homology with the ones presented here. In Drosophila, it has been demonstrated that these genes were overexpressed during parasitoid infection and conferred strong resistance to the host insect compared with individuals in which the gene had been altered by mutagenesis (Verster et al., 2023). Given the increased expression in the more vulnerable stages of the carrot weevil, the hypothesis of a role in parasite defence cannot be excluded here either.

Wolbachia

Based on phylogenomic analyses, the Wolbachia strain isolated from L. oregonensis belongs to supergroup A. Several Wolbachia strains cause cytoplasmic incompatibility, which can have a significant genetic impact by hindering reproductive success in pairings of infected males and uninfected females (Hoffmann et al., 1986). The cytoplasmic incompatibility mechanism is believed to involve the coordinated expression of a toxin‐antitoxin pair, namely the cytoplasmic incompatibility factors CifA and CifB, located in the prophage WO (LePage et al., 2017). In this study, a single CifA‐CifB pair was identified and was found to be similar across all the subpopulations tested. The prevalence was high, with all individuals tested being positive for the presence of Wolbachia. Consequently, no limitations are foreseen in mating compatibility across North American populations of L. oregonensis.

AUTHOR CONTRIBUTIONS

Dave T. Ste‐Croix: Writing – original draft; writing – review and editing; investigation; conceptualization; methodology; validation; visualization; formal analysis; data curation. Annie‐Ève Gagnon: Funding acquisition; conceptualization; writing – review and editing; project administration; supervision; resources. Benjamin Mimee: Conceptualization; investigation; funding acquisition; writing – original draft; methodology; validation; visualization; writing – review and editing; formal analysis; project administration; resources; supervision; data curation.

FUNDING INFORMATION

This work was supported by the Alternative Pest Management Solutions initiative (APMS project J‐002846) of Agriculture and Agri‐Food Canada.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

Supporting information

Table S1. List of samples, use and accession numbers. Table S2. Accession numbers of sequences used for comparative genomics. Table S3. List of predicted genes in Wolbachia with annotations. Table S4. List of predicted genes in Listronotus oregonensis with annotations and their expression across all developmental stages. Table S5. Horizontal gene transfer (HGT) in Listronotus oregonensis. Table S6. Overall statistics from OrthoFinder. Table S7. Per‐species statistics from OrthoFinder. Table S8. Diversity of odorant‐binding protein genes identified in insect species. Table S9. Diversity of odorant‐receptor genes identified in insect species. Table S10. Diversity of gustatory‐receptor genes identified in insect species. Table S11. Diversity of ionotropic‐receptor genes identified in insect species. Table S12. Diversity of cytochrome P450 monooxygenase genes identified in insect species. Table S13. Diversity of carboxyl/cholinesterase genes identified in insect species. Table S14. Diversity of glutathione S‐transferase genes identified in insect species.

Figure S1. Phylogenetic relationship of different Wolbachia isolates based on 14 auto‐selected genes by the autoMLST pipeline. Figure S2. Comparison of the cytolethal distending toxin subunit B (cdtB) genes found in Listronotus oregnensis and closely related species or bacterial species. (a) Amino acid sequence alignment. (b) Phylogenetic tree inferred from a maximum‐likelihood (RAxML) analysis. (c) Pairwise amino acid identity matrix. Figure S3. Comparison of sequence homology for the Wolbachia cifB gene among Listronotus oregonensis individuals collected from fields (F) in Quebec (QC), Ontario (ON), Ohio (OH) and Nova Scotia (NS) or from rearing (R). Nucleotides highlighted in blue are identical across all samples.