Identification and Expression Patterns of Chemosensory Genes in Male and Female Wax Moths, Galleria mellonella

The greater wax moth, Galleria mellonella (Linnaeus, 1758), is one of the most notorious pests of honey bee colonies (Kwadha et al., 2017; Svensson et al., 2014; Dweck et al., 2010). Its larva burrows into the edges of unsealed cells, and feeds on pollen, honey, beeswax and bee brood. The adult wax moth can cause galleriasis by constructing silk-lined tunnels through cell walls and on the comb surface, thus resulting in a significant loss to the beekeeping industry (Bombelli et al., 2017; Kwadha et al., 2017; Krams et al., 2015; Svensson et al., 2014). Besides, an infestation of honey bee colonies by G. mellonella larvae often leads to colony losses, absconding and small-scale migration (Ellis et al., 2013; Kwadha et al., 2017). As a result, G. mellonella-infested honey bee combs in hives exerts a serious side effect on bees (Svensson et al., 2014). Recently, the black queen cell virus and Israeli acute paralysis virus were found in G. mellonella larvae, and the spores of Paenibacillus spp. were also detected in fecal pellets of the larvae (Kwadha et al., 2017; Traiyasut et al., 2016; Hood et al., 2003; Charriere and Imdorf, 1999). This indicates that both G. mellonella adults and larvae are potential vectors of honey bee disease-causing pathogens. The damage caused by this moth is believed to contribute to the substantial decreases in both honey and native honey bee populations (Kwadha et al., 2017; Strauss et al., 2013). Given the growing concern about the health of honey bees and the economic downfall caused by G. mellonella infestation of honey bee colonies, it is necessary to develop alternative strategies for controlling this moth.

As a nocturnal species, freshly eclosed wax moths often fly away from the bee hives to mate at night. After mating, the gravid females re-enter the hives and lay eggs in small cracks or crevices (Kwadha et al., 2017). Upon hatching, G. mellonella larvae move into the bee combs where they begin to feed and, ultimately, damage the comb (Kwadha et al., 2017; Ellis et al., 2013). Therefore, controlling the mating process and egg-laying behaviors of G. mellonella may be a critical strategy to reduce the damage caused by this insect. Unlike most moths, G. mellonella has distinct reproductive behavior where the males secret sex pheromones to attract potential mates (Kwadha et al., 2017; Svensson et al., 2014; Han et al., 2003). Over the past decades, there has been a tremendous advancement in the understanding of insect chemical ecology following the development of pheromone-based trapping systems (Zhang et al., 2015a; Vogt and Riddiford, 1981), which provides an opportunity to develop convenience, cost-effective and sustainable techniques for pest control.

G. mellonella uses a fined-tuned olfactory system located in the sensory hairs (sensilla) on each antenna to detect pheromones and other odors (Grosjean et al., 2011). A diverse range of olfactory proteins, such as chemosensory protein (CSP), ionotropic receptor (IR), odorant-binding protein (OBP), olfactory receptor (OR) and sensory neuron membrane protein (SNMP), are highly expressed in the sensillum lymph (Zhao et al., 2016; He et al., 2015; Sanchez-Gracia et al., 2009). OBPs comprise 6 cysteine residues that form 3 disulfide bonds, thereby generating a hydrophilic pocket that binds to volatile compounds (Pelosi et al., 2014). CSPs are another class of soluble binding proteins enriched in the sensillum lymph, and are consisted of 4 conserved cysteines that form 2 disulfide bonds (Vieira and Rozas, 2011). CSPs are involved in the process of semiochemical detection, and can be found in the chemosensory/non-chemosensory organs (Zhao et al., 2016; Liu et al., 2014; Gu et al., 2012). ORs, seven-transmembrane proteins, can respond to odors and pheromones through coexpression with a conserved co-receptor (Orco) and subsequently trigger signal transduction pathways by converting the chemical signals of the active odor molecules into electrophysiological signals (Cao et al., 2014; Benton et al., 2009). IRs are associated with ionotropic glutamate receptors (iGluRs) and responsible for the recognition of ammonia and acids (Rogers et al., 2001). These reporters have been found in some insect species from different orders (Cao et al., 2014), but have only been comprehensively assessed in Drosophila melanogaster. SNMPs, which contain two orthologs of SNMP1 and SNMP2, can facilitate the ligand delivery to receptors (Rogers et al., 2001). SNMP1 is specifically expressed in the pheromone-responsive olfactory receptor neurons of D. melanogaster, and is responsible for mediating the sensitivity of these neurons to cis-vaccenyl acetate stimulation (He et al., 2019; Cao et al., 2014; Liu et al., 2012). Thus, identification of these chemosensory genes in G. mellonella can help to improve the current moth trapping systems.

In the past few years, antennal transcriptome sequencing has been successfully employed to detect a number of candidate olfactory genes in lepidopterans, including Bombyx mori (Fang et al., 2015), Chilo suppressalis (Cao et al., 2014), Sesamia inferens (Zhang et al., 2014), Loxostege sticticalis (Wei et al., 2017) and Mythimna separata (Chang et al., 2017). Besides, several other moth species, such as the hymenopteran Apis cerana, hemipteran Adelphocoris suturalis and dipteran Bactrocera dorsalis have been sequenced by Zhao et al. (2016), Cui et al. (2017) and Jin et al. (2017), respectively. More recently, Zhao et al. (2019) performed transcriptome sequencing on G. mellonella antennae and identified several chemosensory genes (e.g., 46 ORs, 25 IRs, 22 OBPs, 20 CSPs and 2 SNMPs). Similarly, Lizana et al. (2020) also identified 20 OBP genes from this moth species. However, the number of OBPs in G. mellonella is relatively lesser compared to other Pyralid moths, such as C. suppressalis, Ostrinia furnacalis, and Cnaphalocrocis medinalis, with 23, 26 and 26 OBP genes, respectively (Liu et al., 2017; Zhang et al., 2015b; Cao et al., 2014). Hence, the present research aimed to examine the differential expression patterns of chemosensory genes in the antennae of female and male G. mellonella, which will serve as a guide for future works, particularly those on the prevention and/or treatment of moth infestation using male pheromones as a bait.

G. mellonella pupae and larvae were sampled from the infected hives in an apiary at Shanxi Agricultural University (Shanxi, China). Newly emerged adult wax moths were used to preserve the stock culture. All G. mellonella larvae were reared on old honeycombs in an incubator at 34 ± 1°C with 65 ± 5% relative humidity in constant darkness. For transcriptomic sequencing, approximately 100 antennae per sex were collected from freshly emerged female and male moths (three replicates per sex). All samples were quickly frozen in liquid nitrogen and then kept at -80°C for later use.

Total RNA was isolated from the antennae using TRIzol™ Reagent (Invitrogen, CA, USA) in compliance with the manufacturer’s instructions. The residual DNA was removed using DNase I (Promega, WI, USA). RNA quality was assessed using a 1% agarose gel dissolved in electrophoresis buffer. RNA-Seq library preparation and Illumina sequencing were carried out by Novogene Bioinformatic Technology (Beijing, China).

Raw data in FASTQ format were processed by PERL scripts built in house. After eliminating the low-quality reads and those with adapter sequences or poly-N, clean sequences were retained. Next, the sequence duplication level as well as the Q20, Q30, and GC content of the high-quality reads were determined. Finally, clean data were assembled into unigenes using Trinity software (r2014-04-13p1; Trinity Software Solutions Inc., Waterford, VA, USA) (Grabherr et al., 2011), and each unigene was assigned a unique gene ID. All downstream analyses were then performed on the high-quality sequences. The Illumina sequencing data generated in this study were submitted to the sequence read archive (SRA) of the National Center for Biotechnology Information (NCBI) (accession No. SRR11446320).

The functions of new genes were analyzed by conducting a basic local alignment search tool (BLAST) search against seven databases, including Gene ontology (GO), NCBI nucleotide sequences (Nt), NCBI non-redundant protein (Nr), Swiss-Prot, protein families (Pfam), Kyoto encyclopedia of genes and genomes orthology (KO), and euKaryotic ortholog groups (KOGs). Proteins with the highest sequence identity for a given unigene were retrieved together with their respective functional annotation categories. The open reading frames (ORFs) of the candidate chemosensory genes were analyzed by the NCBI ORFfinder. The putative N-terminal signal peptides of OBP and CSP genes were estimated using the SignalP 5.0 Server. The conserved domains of OBP and CSP genes were predicted using the NCBI Conserved Domain Database. The transmembrane domains (TMDs) of IR, OR and SNMP genes were determined through the TMHMM version 2.0 Server. Subsequently, a phylogenetic tree was built according to the amino acid sequences of putative chemosensory genes in G. mellonella as well as the homologous sequences found in other lepidopteran species. After aligning the amino acid sequences through ClustalW Version 2.1, an unrooted neighbor-joining tree was established by MEGA5.2, and the branch support was evaluated with 1,000 bootstrap replications (Tamura et al., 2011).

HTSeq version 0.9.1 was applied to estimate the expression levels of candidate genes (Trapnell et al., 2010). After that, the Fragments Per Kilobase per Million mapped reads (FPKM) value of each chemosensory gene was calculated based on the following equation:

FPKM = (1,000,000 × C) / (N × L× 1000)

Where, C represents the number of reads uniquely aligned to the chemosensory gene, N denotes the total number of fragments mapped to all unigenes, and L indicates the number of bases in each chemosensory gene (Mortazavi et al., 2008).

RNA-Seq by Expectation-Maximization (RSEM) software was used to map the clean sequences to the transcriptomic unigenes based on the default settings (Li and Dewey, 2011). The differential expression between two transcriptomes (duplicate biological samples) was analyzed by DESeq Version 1.12.0 (Wang et al., 2010). Differentially expressed genes were those with an adjusted P-value of <0.05. GO enrichment and KEGG pathway analyses were conducted using the GOseq R package version 1.10.0 and KOBAS Version 3.0, respectively. Online tools and databases used in this study are listed in Table I.

Quantitative real-time polymerase chain reaction (qRT-PCR) was carried out to validate the levels of differentially expressed genes (DEGs). The specific primers were designed with Primer3.0Plus server,and elongation factor 1-alpha (Ef-1a) was employed as an internal standard for data normalization (Table II). RNA extraction and cDNA synthesis were performed with TRIzol™ Reagent (Invitrogen) and PrimeScript™ RT Reagent Kit with gDNA Eraser (Takara Bio Inc., Dalian, China), respectively, by following the manufacturer’s protocols. In PCR, 2 μL of 1:3 diluted cDNA was used as the template, with a reaction mixture of total volume of 20 μL. The qRT-PCR was conducted on an Applied Biosystems™ 7500 Fast Real-Time PCR System (Applied Biosystems, CA, USA), with the following reaction conditions: denaturation for 4 min at 95°C, followed by 40 cycles of 15 s at 95°C and 34 s at 60°C. The specificity of the qRT-PCR reaction was assessed by performing a

melting curve analysis of 20 s at 95°C, 30 s at 60°C, and 30 s at 95°C. Each candidate gene had 3 biological replicates, and 3 technical replicates were tested for each biological replicate. The relative quantification of PCR results was performed using the comparative CT method (2-ΔΔCt).

Data of the antennal transcriptome of G. mellonella were produced using the Illumina Hiseq™ 2500 platform (Illumina Inc.). In total, 135,622,632 and 134,292,198 raw reads were acquired from the libraries of female and male specimens, respectively. After eliminating adapters, ambiguous nucleotides and low-quality sequences, the female and male antennae yielded 132,262,518 and 130,841,026 clean reads, respectively. These clean reads were assembled into 372,571 transcripts, which accounted for 39.47 gigabases with a GC percentage of 41.78% (Table III). After merging and clustering, 188,278 unigenes were obtained (mean length = 781 bp and N50 length = 1161 bp), and 107,679 unigenes (~57%) were 200–500 bp in length. After annotation by tBLASTn, 108,047 (57.38%) unigenes were allocated to more than 1 database, and 7388 (3.92%) were allocated to all the 7 databases.

The 5 chemosensory gene groups (CSPs, IRs OBPs, ORs and SNMPs) were identified through keyword searching and manual analysis of annotated unigenes. In total, 121 chemosensory genes were annotated in the antennal transcriptome of G. mellonella, comprising 37 OBPs, 35 CSPs, 33 ORs, 14 IRs, and 2 SNMPs. The information on these unigenes, such as the gene ID, amino acid sequence length, ORF length and FPKM value, are provided in Table IV.

Among the 37 candidate OBP genes, 31 OBPs had a full-length ORF with signal peptide prediction sequences, whereas the remaining 6 OBPs (OBP6, OBP8, OBP12, OBP14, OBP15 and OBP22) corresponded to a partial sequence encoding 87–147 amino acids. Among the 31 full-length OBPs, OBP1, OBP2, OBP3, OBP4, OBP5, OBP7, OBP9, OBP11, OBP13, OBP16, OBP17, OBP23, OBP72-like, OBP69a-like, GOBP2, GOBP3, PBP1, PBP2, PBP3, PBP4, PBP5 and PBP6 were the typical OBPs with 6 cysteines; OBP18, OBP19 and OBP56d belonged to the Minus-C OBPs without the 2nd and 5th cysteines; OBP10, OBP21, OBP83a-like, OBP84a and GOBP1 were the Plus-C OBPs with additional cysteine residues. Figure 1 showed the phylogenetic tree constructed based on the amino acid sequences of OBPs in C. medinalis, C. suppressalis, O. furnacalis, Conogethes punctiferalisand G. mellonella. As expected, 37 Gmel OBPs were divided into 5 groups, namely, GOBP, PBP, typical OBPs, Plus-C OBPs and Minus-C OBPs.

C_1, C_2, C_3− biological replicate samples of female G. mellonella. X_1, X_2, X_3 − biological replicate samples of male G. mellonella.

A total of 35 unigenes were annotated to putative CSPs, all of which were predicted to have four cysteines, and they consisted of a full-length ORF encoding 97–147 amino acids. The neighbor-joining tree analysis revealed that all 35 sequences were clustered with at least 1 lepidopteran orthologous gene, and the CSPs were clearly observed (Fig. 2). The unigenes corresponding to CSPs were designated according to the obtained CSP data.

A total of 33 ORs were identified from the assembled unigenes, which belonged to the seven-transmembrane receptor superfamily. Among them, eight ORs (Orco, OR7, OR10, OR12, OR13a-like, OR29, OR67c-like and OR94a-like) contained full-length ORFs with 5-7 transmembrane domains. As could be seen from the phylogenetic tree in Figure 3, the OR2 sequence was 99% identical to CpunOR2 and CmedOrco, and thus we labeled it as GmelOrco, while other ORs were classified into distinct clades with known ORs.

We identified 14 transcripts encoding candidate IRs, including IR8a and IR25a (members of highly conserved IR co-receptors). Among these, 5 IRs (IR8a, IR21a, IR25a, IR64a and IR76b) contained full-length ORFs encoding a protein of 547 amino acids. As shown in Figure 4, most IRs were clustered together with their orthologs into a distinct clade.

SNMP transcripts with two TMDs were identified in the transcriptome of G. mellonella, including SNMP1 and SNMP2. Moreover, intact ORFs with 525 and 521 amino acids were observed for SNMP1 and SNMP2, respectively. Figure 5 shows the phylogenetic tree constructed according to 31 SNMP sequences from 19 species. Notably, the insect SNMPs were assigned to two highly conserved, distinguishable classes, namely, SNMP1 and SNMP2.

The expression levels of the 121 chemosensory unigenes in six cDNA libraries were determined with the FPKM index. The differential expression profiles revealed that several chemosensory genes (CSP1, CSP2, CSP3, CSP4, CSP6, CSP7, CSP10, CSP28, CSP33, OBP1, OBP11, OBP69a-like, OBP84a, GOBP1, GOBP2, PBP2 and PBP3) were overexpressed in female and male antennae (FPKM > 100), and CSP6 exhibited the highest expression (FPKM = 31189.81). Moreover, some genes were specifically detected in the sexual state, although their expression was extremely low. For example, CSP8, CSP25, CSP29, OR6, OR7, OR12, OR38, OR85c-like, OR85e, OR92a, OR94b-like and IR75p were expressed exclusively in the female antennae, while OR21 was detected only in the male antennae. Furthermore, OR, IR and SNMP genes exhibited relatively low expression in each sample (FPKM < 100), and the two SNMPs demonstrated higher FPKM values in the female antennae than in the male antennae.

Next, DEG analysis was carried out, and the Venn diagram in Figure 6A shows the overlap of 62,445 genes between the two groups. Notably, 57,407 genes were expressed exclusively in the female antennae, while 26,370 were found only in the male antennae. Moreover, 114 DEGs were identified, of which 5 were chemosensory-related genes such as 2 CSPs (CSP7 and CSP10), 2 OBPs (OBP72-like and OBP69a-like) and 1 OR (OR29). Among these DEGs, 66 were upregulated and 48 were downregulated in the female antennae (Fig. 6B). The cluster analysis revealed that the DEGs between females (C) and males (X) were assigned to two main groups: (i) gene upregulation and (ii) gene downregulation (Fig. 6C).

The GO enrichment analysis of DEGs was performed using the GOseq method with Wallenius non-central hypergeometric distribution (Young et al., 2012). The DEGs were mostly enriched in the “binding” category, followed by “metabolic process”, “cellular process” and “single-organism process” categories (Fig. 7). More olfactory-related GO terms were related to the upregulated DEGs compared with the downregulated DEGs. These olfactory-related GO terms were associated with the molecular function (binding, catalytic activity, protein binding and ion binding), biological process (metabolic process, cellular process and single-organism process), and cellular component (cell and cell part).

In addition, our results demonstrated that all DEGs were mapped to 28 reference KEGG pathways. These DEGs were remarkably enriched in “ubiquitin mediated proteolysis”, “signaling pathways regulating pluripotency of stem cells”, and “renal cell carcinoma and pathways in cancer” (Fig. 8). Besides, several enriched pathways were associated with cAMP, hypoxia-inducible factor- 1 (HIF-1), Hippo, mitogen-activated protein kinase (MAPK), T cell receptor, and glucagon signaling pathways.

To verify the results of the transcriptome analysis, we carried out qRT-PCR to detect the expression of 16 DEGs (2 CSPs, 2 OBPs, 1 OR, 1 gustatory receptor, 2 P450s, 1 juvenile hormone binding protein and 7 other genes). The results demonstrated that the expression patterns of 15 out of the 16 selected genes were consistent with those generated from the RNA-Seq analysis (Table V), thus implying the reliability of the sequencing data.

Pollinators provide an important ecosystem service by enhancing the yields of wild and crop plants globally (Li et al., 2015). Increasing evidence has shown the declines in both wild and domesticated insect pollinators (Kwadha et al., 2017). G. mellonella is considered a key factor for the decreases in both native and feral honey bee colonies, especially in tropical and sub-tropical regions (Kwadha et al., 2017). Olfaction is an important attribute of smell for the greater wax moths to reproduce and survive (Zhao et al., 2019). In the present work, the antennal transcriptomes of both sexes of the greater wax moths were sequenced and analyzed. We identified 121 chemosensory gene transcripts (37 OBPs, 35 CSPs, 33 ORs, 14 IRs and 2 SNMPs), and which might provide base data to elucidate the olfactory recognition mechanism of G. mellonella. Numerous candidate chemosensory genes identified in this study are comparable with those reported in a recent study (Zhao et al., 2019), with 22 OBPs, 20 CSPs, 46 ORs, 25 IRs and 2 SNMPs. The reduced number of ORs and IRs may be attributed to different sampling time and G. mellonella life-cycle, as the greater wax moth are more active during the first three days and their activities reduce with an increase in age (Li et al., 2019).

Courtship behavior is essential for animal reproduction (Grosjean et al., 2011; Jiang et al., 2016). To mate, animals have evolved a wide variety of pheromone release and detection patterns (Zhang et al., 2015c). G. mellonella is an interesting insect to study because males attract females (Kwadha et al., 2017). Pheromone-binding proteins (PBPs) are a subgroup of OBPs, which play a crucial role in regulating olfactory process (Liu et al., 2012). In this research, we identified six transcripts (PBP1, PBP2, PBP3, PBP4, PBP5 and PBP6) encoding candidate PBP genes based on their similarity with the PBPs from other lepidopterans and the physiological analysis. The number of PBP genes identified was slightly higher than that reported in other lepidopteran insects, and these six sequences were detected in the antennae of both female and male wax moths. Thus, we speculate that the PBPs in G. mellonella are different from those in other lepidopteran species, as the expression of PBPs is male-biased and antenna-predominant in most lepidopterans during the synthesis of male pheromones (Zhao et al., 2019). However, the function of PBPs in G. mellonella remains to be fully elucidated.

SNMPs were first discovered in pheromone-responsive neurons, which could influence the detection of pheromones (Rogers et al., 2001). At present, the molecular mechanism underlying the functions of insect SNMPs is poorly understood (Li et al., 2015). In this study, we successfully identified two SNMPs (SNMP1 and SNMP2) that were not differentially expressed between the female and male antennae, although a higher FPKM value was observed in the female antennae than in the male antennae. These findings are consistent with those reported in other known lepidopterans, suggesting that the SNMPs in G. mellonella may play a similar role as in D. melanogaster and other moths.

To assess the differential expression patterns of chemosensory-related genes in both female and male antennae, RNA-Seq was performed to compare the levels of DEGs. The results indicated that the number of downregulated DEGs (66) was slightly lower than that of upregulated DEGs (48) between the two samples. Among these DEGs, we identified three (OBP72-like, OBP69a-like and OR29) and two (CSP7 and CSP10) genes with remarkably higher expression levels in the female and male antennae, respectively. According to the functions of insect OBPs, CSPs and ORs (Li et al., 2015), the female-biased OBP72-like, OBP69a-like and OR29 genes are responsible for the detection of sex pheromones released by males or odors critical to female-specific behaviors (e.g., searching bee-comb hosts for oviposition), while male-biased CSP7 and CSP10 genes can detect odors critical to male-specific behaviors. Nevertheless, the sex-specific functions of these chemosensory-related DEGs should be further investigated.

GO analysis demonstrated that the 114 DEGs were most significantly enriched in molecular function (e.g., binding). These included protein binding, ion binding, heterocyclic compounding, and organic cyclic compound binding. Considering that the antennae are critical olfactory appendages in an olfactory system that interact with different types of chemical stimuli (Zhao et al., 2019), it is reasonable that the identified DEGs are involved in binding. The KEGG analysis revealed that signal transduction pathways (e.g., stem cells, T cell receptor, MAPK, glucagon, calcium, HIF-1, cAMP and Hippo signaling pathways) were found to be markedly enriched and closely associated with the antenna (as the primary organ of odor binding and signal transduction).

The analysis of antennal transcriptomic data revealed 121 putative chemosensory genes in G. mellonella and 114 DEGs between the male and female antennae. Furthermore, our method successfully detected chemosensory genes with very low expression, which could provide essential information to further investigate the underlying olfactory recognition mechanism in G. mellonella and serve as a platform for further functional analyses of the related genes. This may ultimately lead to the identification of a suitable male pheromone that can be used as a bait for trapping G. mellonella. Such pest management strategy can treat and/or prevent moth infestations of honey bees and, consequently, improve the health of honey bee colonies worldwide.

This work was supported by the National Natural Science Foundation of China (31640080) and China Agriculture Research System (CARS-44-SYZ 16). We acknowledge the help of Xinyu Li, Wenting Su, Yujia Feng, and Denglong Long with the experiments. We greatly appreciate the assistance of the scientists of Novogene Bioinformatics Technology Co., Ltd. in sequencing the DEG libraries. We also thank EditSprings for providing language editing services.

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.