Whole-genome identification of chemosensory genes in two sister blister beetles provides insights into chemosensory evolution

Background: Commonly known as blister beetles or Spanish fly, there are more than 1,500 species in the Meloidae family (Coleoptera:Tenebrionoidea) that produce the potent defensive blistering agent cantharidin, which has recently been exploited for cancer therapy. Hycleus cichorii and Hycleus phaleratus are exploited as traditional Chinese medicine over 2,000 years due to their ability to biosynthesize cantharidin. These blister beetles share highly similar environments and ecological niches. To understand the role of the chemosensory system in speciation and evolution in the beetles, we identified the chemosensory gene families in whole genome of both blister beetle comprehensively. Results: We identified 29 OBPs, 10 CSPs, 116 ORs, 80 GRs, and 35 IRs in H. phaleratus and 30 OBPs, 10 CSPs, 191 ORs, 92 GRs and 42 IRs in H. cichorii. The 7 groups of beetles’ ORs were recovered, but a new pattern of ORs in Coleoptera were surfaced, due to lack of whole OR repertoire for comparison to T. castaneum before. The OR number is varied, identified in H. cichorii 64% more than that in H. phaleratus (191 VS 116), which maybe an evolution event between these blister beetle’s speciation. A major Hycleus-specific expansion clade of bitter GRs is evident based on our phylogenetic tree, in which, clustered more than 50% Hycleus-bitter GRs. All the ten clade of beetles’ antennal IR identified in both Hycleus beetles. Interestingly, IR75 and IR41a were obviously expensed compare to other insects. Conclusions: The GR and IR expansion events may promote Hycleus genus evolution. Our data will provide a basis for future species protection and speciation focusing on the blister beetle.


Background
To adapt to varied environments and ecological niches, insects always rely upon their sensitive olfactory system to distinguish crucial chemical signals, such as pheromones and plant volatiles. Because olfaction is an important element of host and mate selection, the identification of changes in chemosensory genes between diverging populations or closely related species can provide valuable insights into the role of chemosensory adaptation in different environments and host speciation [1]. At the species level and above, gene amplification and amino acid mutation events may play essential roles during adaptive evolution. For example, some gustatory receptor (GR) and odorant receptor (OR) genes underwent rapid evolution during host specialization in Drosophila [2], and certain chemoreceptor genes with genetic divergence might play key roles in local adaptation and reproductive isolation in the pea aphid Acyrthosiphon pisum [3]. Using developing expression analysis and targeted mutagenesis experiments, Matsuo et al. [4] identified two genes encoding odorant binding proteins (OBPs) that were responsible for behavioral differences between Drosophila sechellia and its sibling species in their responses to taste perception and host plant preferences. Moreover, the evolution of differential gene expression may act as a key mechanism of chemosensory divergence. For instance, a comparison of chemosensory gene expression profiles of D. sechellia and two sibling species showed significant changes in gene expression that might be related to host specialization [5].
In insects, at least five gene families are involved in the detection of chemical signals.  [11][12][13][14], contact chemicals or carbon dioxide (GRs) [15,16], and nitrogen-containing compounds, acids, and aromatics (IRs) [17]. In contrast, the binding protein gene families are highly abundant in the sensillum lymph and usually function as binding of molecules that interact with receptors [18,19]. Together, these genes allow insects to seek and select hosts and mates and thus affect the adaptability changes in closely related insect species. Although several of studies have compared differences in chemosensory genes between sister species [20][21][22][23][24][25][26][27], comparatively few cases have aimed to understand the mechanisms that have shaped divergence or chemosensory speciation. Chemical cues with diverse ecological functions can be detected by focusing on the chemosensory genes expressed in the insects [28,29].
Moreover, key genetic factors can impact insect behavior and host use.
Commonly known as blister beetles or Spanish fly, there are more than 2,500 species in the Meloidae family, with more than 1,500 of these beetle species known to produce cantharidin. Cantharidin (C 10 H 12  These two beetles have largely overlapping sympatric ranges in China and a similar emergence phenology and appearance, except that H. phaleratus has a bigger body size.
The beetles are mainly distributed in southwest China and feed on the eggs of Locustidae during the larvae period and on legumes and cucurbitaceous crops in adulthood.
Therefore, understanding the changes that underlie the speciation between these two sister species and comparing the chemosensory proteins of extremely importance for the insects' survival and reproduction are interesting strategy. On other hand, in Coleoptera groups, the chemosensory genes have rarely been identified in whole genome except the flour beetle T. castaneum [37,38], lack of enough data for comprehensive analysis of gene evolution. Here, we identified the members of five gene families involved in chemosensory perception in Hycleus genomes and analysed its evolution in beetles.

Identification of chemosensory genes
The H. cichorii and H. phaleratus genome sequences were obtained from the GigaScience Database [39]. We identified the OBP, CSP, OR, GR, and IR genes using the following steps. First, we constructed a protein reference database of each family that included two parts: Dendroctonus ponderosae, Tribolium castaneum, Pyrrhalta aenescens and Pyrrhalta maculicollis, Chrysomela lapponica from GenBank and Acyrthosiphon pisum, Drosophila melanogaster, Bombyx mori, Apis mellifera, Dendroctonus ponderosae, and Tribolium castaneum from InsectBase. The second, a tBLASTn (Version 2.2.25) alignment was used to search the chemosensory genes onto both genome with an e-value cut-off of 1e-5. The third, we first linked the alignment hits into candidate gene loci exon by exon with house Perl script. We then extracted genomic sequences of candidate loci, together with 1 kb flanking sequences, using Gene-Wise [40] to determine gene models. Next, the conflict region was manually compared with other gene of the same family. At last, we checked the functional annotations of the public data and removed candidate genes that had inconsistent functions. Finally, repeat the whole processes used the identified putative proteins to search again, incase missing the new genes.

Phylogenetic tree and selection analysis
To analyze the characteristics of the Hycleus chemosensory genes and their relationships to other insects, ML trees for each family were constructed using the amino acid sequences derived from the putative CDSs and the published sequences of other species, including T. castaneum, D. ponderosae, C. lapponica and D. melanogaster. The published sequences were retrieved from NCBI and InsectBase (S1 Dataset). First, a multiple sequence alignment was generated through MAFFT (Version 7. 407) [41] using the chemosensory gene amino acid sequences of 6 species. Then, RAxML (v8.2.12) [42] was used to construct chemosensory phylogenetic trees with best models (LG for CSPs and OBPs, WAG for IPs, JTT for ORs and GRs), and node support was assessed using a bootstrap procedure based on 2,00 replicates. We viewed the phylogenetic tree with EvolView (http://evolgenius.info/evolview).
According to the phylogenetic tree, we selected the orthologous genes between these two Hycleus species to process in the selection analysis. The paired orthologs of the protein sequences were aligned using MAFFT (version 7.407) [41]. Then, the alignment information for the protein sequence was used as a guide for the nucleic acid sequence alignment. Finally, analyses of nonsynonymous (KA) and synonymous (KS) substitutions were performed by Kaks_Calculator (v2.0) [43] using the YN method.

Results
The annotated candidate chemosensory genes were identified by Blastp

OBPs.
We identified 29 and 30 OBPs in the two blister beetles, respectively; and 48 out of 59 OBPs with nearly full-length CDSs. We identified 28 classic, 4 plus-C and 27 minus-C OBPs based on their sequence properties [18,46,47]. The maximum likelihood (ML) phylogenetic tree suggested 28 pairs of orthologous OBPs between the two beetles ( Fig   S1). Analysis of selection performed on the orthologs between the two blister beetles suggested that 13 pairs ortholog (dN/dS ≤ 0.3436, Fisher's test P ≤ 0.05) was under purifying selection (Additional files 1).

Chemosensory proteins.
We identified ten CSPs in H. phaleratus and H. cichorii, respectively; all the CSPs represented full-length proteins including the conserved C-pattern, C1X 6 C2X 18 C3X 2 C4 [47]. Nine pairs of chemosensory protein orthologs were shared by the two blister beetles ( Fig S2). The selection pressure analysis showed that purifying selection acted on six CSP orthologs (Additional files 1).

Odorant receptors.
The coleopteran OR subfamilies were classified by Engsontia et al.
Group 1 and 6 contained ORs were mainly Tribolium castaneum-specific except for a few ORs. Group 2 contained a mixture of ORs from all five species. Group 3 and 5 only contain ORs mainly from Tribolium castaneum and Hycleus genus. Groups 4 were Tribolium castaneum-specific clade. Otherwise, group 7 was a specific expansion in bark beetles (DponOR) and blister beetles, and this clade exception ORs of Tribolium castaneum (Fig.  1). Specifically, HcicORco and HphaORco were identified, and ORco homolog group was clustered into group 7 (Fig. 1).
Notably, we identified 116 and 191 ORs in H. phaleratus and H. cichorii genome, and up to 64% ORs identified in H. cichorii more than that in its sister species. This expensed ORs may acts a key functional along to its species speciation.
We identified 86 pairs of orthologs between the two beetles based on the ML tree, in which, 49 pairs of orthologs were under purifying selection (Additional files 1).

Gustatory receptors.
We identified 80 and 92 GRs in H. phaleratus and H. cichorii, respectively. And 63 pairs of GRs ortholog were shared by the two blister beetles (Fig 2). Three pairs of CO 2 receptors were identified, and one by one orthologs presents among Hycleus genus and T. castaneum (Fig. 2). Eight and seven sugar receptors, which was in DmelGR64a-f cluster, were detected in this study (Fig. 2). Also, a pair of homologs (HcicGR56 and HphaGR56) were homolog to DmGR43a, which functions as a fructose receptor. Notably, a huge Hycleus-specific expansion clade is evident in our phylogenetic tree (Fig. 2). More than 1/2 identified Hycleus bitter GR were classified in this clade, that may act important functions in Hycleus genus speciation. Moreover, three species (or genus) expansion clade, Hycleus genus, T.castaneum, and D. melanogaster, were identified (Fig. 2). A total of 40 pairs of orthologs between blister beetles were under purifying selection (Additional files 1).

Ionotropic receptors.
IRs are more closely related to ionotropic glutamate receptors (iGluRs) [38,48]. A total of 35 and 42 IRs and 13 and 13 iGluRs were identified in H. phaleratus and H. cichorii, respectively ( Table 1). The DmelIR25a and DmelIR8a was considered the coreceptor IRs, and their orthologs of both species of Hycleus were identified in this study (Fig. 3). It is very conserved from beetle to fruit fly, a 1:1 orthologs appeared in our ML tree. This cluster are in the clade of iGluRs, corresponds to previous report [48,49]. All the ten of beetle's antenna IRs, including IR75, IR64a, IR40a, IR21a, IR68a, IR60a, IR41a, IR76b and IR93a, were identified (Fig 3 and Table 2). Eight of ten clades of these only detected one copy in both blister beetles, and five of that were 1:1 ortholog in all the six insects (Fig 3 and Table 2). Notably, a total of 11 and 15 copies of IR41a were annotated in H.
phaleratus and H. cichorii, respectively. However, that corresponding IR in other species no more than 2 copies (Fig 3 and Table 2). It was expensed obviously in blister beetles compare to other insects. IR41a is not a unique instance. Also, the IR75 were expensed in both blister beetles relatively, 10 and 16 copies in these two beetles but up to 4 copies in other insects (Fig 3 and Table 2).
A total of 35 pairs of orthologous iGluRs/IRs were found between both species, and 30 pairs of orthologs were under purifying selection (Additional files 1).

Discussion
Chemosensory traits are key components of speciation in many insects, and changes in these traits can influence insect behavior and subsequently mediate survivability.
Although evaluating changes in chemosensory gene is a necessary beginning for understanding how chemosensory traits contribute to adaptability, surprisingly few studies have compared how sister taxa differ in this field. Moreover, chemosensory gene identification within Coleoptera groups has rarely been reported. These genes have been characterized in whole genome level only the flour beetle T. castaneum [37]. Several species have been identified partly by RNA-seq, such as two bark beetles (I. typographus and D. ponderosae) [20], two leaf beetles (P. aenescens and P. maculicollis) [27], scarab beetle Hylamorpha elegans [50], and a leaf beetle C. lapponica [49]. In the present study, [27] and the leaf beetle C. lapponica (32 OBPs and 12 CSPs) [49]. An expansion of the minus-C OBPs was found in this Hycleus genus (27 of minus-C OBPs) comparison to D. melanogaster, the similar pattern reported in Tribolium and bark beetles [20], two beetles in the Pyrrhalta genus [27] and the leaf beetle C. lapponica [49]. The expanded Minus-C OBPs may one specific evolution event in Coleoptera chemosensory, and it suggests that these genes might play an important role in chemosensation.

Receptor gene families
The OR family, which was studied mostly in these three receptors. Based on previous studies [20,37], seven groups of beetle ORs were recovered in this study. Notably, some new patterns have been found in this study (Fig. 1). Such as, the group 7 was a specific expansion only in bark beetles (DponOR and ItypOR) was reported before [20]; this clade was also expansion in blister beetles was be found in our result. The group 3 and 5 almost the flour beetle specific ORs was reported before [20,25,27], but we found lots of Hycleus ORs also in these clades. This may due to the most reports were identified ORs via RNA-seq except T. castaneum (341) [37]. These species only tens of ORs could been found in RNA-seq, such as the leaf beetles P. maculicollis (22) and P. aenescens26) [27], the bark beetles I. typographus (43) and D. ponderosae (49) [20], the long horned beetle Megacyllene caryae (57) [25], and the leaf beetle C. lapponica (42) [49]. This indicted more beetle's chemosensory family should be annotated comprehensively via whole genome that is better for comprehensive study its evolution in insects.
Moreover, the number of ORs in each blister beetles was obvious differ, ORs identified in H. cichorii 64% more than that in H. phaleratus (191 VS 116). However, the number of others only slightly vary (Table 1). These data hint ORs specific expansion in H. cichorii along both specie differentiation, which occurred around 23 million year ago [56]. ORs is one of the key bridges between animals and its surrounding, and act functions with food, predators and potential mates. These two beetles have largely overlapping sympatric ranges in west south of China. To our field survey experience, in the past years, the H. phaleratus population has declined in the field due to destruction of the living environment by human activity and human catching. In contrast, the H. cichorii population has not declined obviously as the same situation. These ORs maybe acts the keys roles for a stronger adaption ability than H. phaleratus. These ORs needs verified in the future.
For GRs, always divides into three groups by its recognizing materials, CO 2 , sugar and better receptors. we identified three pairs of CO 2 receptors, and nine and eight sugar receptors in two blister beetles (Fig 2). These two groups always more conservative in insects. Notably, a Hycleus-specific expansion clade was found in our phylogenetic tree, and more than 1/2 identified Hycleus-bitter GR were in this clade. That may act important functions in Hycleus genus. On the other hand, there are lack whole GRs identify from genome of other species for comparative analysis. Most published GR family by RNA-seq, only few GRs was been detected due to the GRs expressed very low. such as in P. aenescens (16), P. maculicollis (10) [27], I. typographus (2) and D. ponderosae (6) [20], the leaf beetle C. lapponica (8) [49].
IRs is also membrane proteins, are more closely related to ionotropic glutamate receptors (iGluRs) [38,48].We identified 35 and 42 IRs in H. phaleratus and H. cichorii. IRs could be divided into two groups: 'antennal IRs' was always conserved among species, which likely define the olfactory receptor family of insects, and the 'divergent IRs' was more divergent, which are expressed in peripheral and internal gustatory neurons, implicating this family in taste and food assessment [38]. All the ten clades of beetle antennal IRs were identified in both blister beetles, and five of these sub-family were conserved with 1:1 ortholog in all the six insects ( Fig. 3 and Table 2). However, the IR41a (11 and 15 in blister beetles VS 1-2 in others) and IR75 (10 and 16 in blister beetles VS 3-4 in others) were evidently expensed compare to others (Fig 3 and Table 2). Moreover, there have 5 out 11 copies of IR41a and 7 out of 10 copies of IR75 were detected expression (FPKM ≥1). This data hints that IR41a and IR75 would be specific expansion in the Hycleus genus, Candidate chemosensory genes for speciation Divergent chemosensory genes linked to differences in the sensory tuning of sister species represent a response to divergent selection of chemosensory traits [1]; thus, these genes are candidates for chemical-mediated speciation, which directly or indirectly reflects the ability to adapt to the environment. For instance, a comparison of the gene amplification of chemosensory gene families between the sister beetle species suggested some candidate genes that might have been involved in this ability. In present study, an evident expansion of ORs identified in H. cichorii compare to that in H. phaleratus (191 VS 116). These data hint ORs specific expansion in H. cichorii may promoted both specie differentiation. Which specific ORs and how acts the keys role should be study in the next.
On the other hands, we identified 221 pairs of gene between both blister beetles, of which 138 pairs were undergone purifying selection and no gene under positive selection (Additional files 1). This data may hint no positive clue of gene mutation linking to the speciation, the other ways should be attended, such as gene expression change. At least, expression changes between orthologs should be attended in the future study.

Conclusions
In the present study, we firstly performed a comprehensive whole genome identification of a pair of blister beetles (Meloidae: Hycleus) with a focus on chemosensory gene families. A New pattern of ORs in Coleoptera were surfaced, due to lack of whole OR repertoire for comparison to T. castaneum before. Several fascinating evolution events in blister beetles was disclosed, such as a clade of Hycleus-specific expansion bitter-GR and IR75 and IR41a evident expansion in Hycleus, which will help to unfold this genus evolution.      Red: 70%, yellow: 50-70% and grey: 50% support.