Introduction

In recent years, the extensive use of pesticide-based pest management has led to drastic effects on our ecosystem and environment1, 2. Therefore, integrated pest management (IPM) has gained great attention as a strategy to protect crops from pest losses3, 4. Biological control is a key agro-system service and a pillar of IPM5, 6. Because of the key advantage of their capacity to both kill and reproduce at the expense of their hosts, parasitoids have been widely used in biological control against insect pests7, 8.

Aphids are key insect pests that are responsible for major agricultural losses, particularly as they are vectors of various plant viruses9. Aphids can be attacked by a wide variety of natural enemies, including several endoparasitoids. Different parasitic wasps usually have distinct hosts, and they can locate their target hosts accurately and efficiently6, 9. The success of parasitoids in locating their hosts in a complex environment depends mainly on the accurate recognition of a series of several chemical molecules. In most insects, the chemosensory system is involved in foraging, oviposition site selection, mate choice, and social communication (among social insects)10,11,12,13,14,15. It is likely that the semiochemicals from foods, hosts, mates, or partners are received by insect chemoreceptors at the membrane surface of chemosensory neurons such as olfactory receptor neurons (ORNs) and gustatory receptor neurons (GRNs)16,17,18. Chemoreceptors include three large, distinct families: odorant receptors (ORs), gustatory receptors (GRs), and ionotropic receptors (IRs)19.

Insect ORs were the first chemoreceptor family to be discovered in Drosophila melanogaster genome20. Until date, insect ORs have been identified in many species, including Apis mellifera 21, Macrocentrus cingulum 10, Conogethes punctiferalis 22, and Bombyx mori 23, with a high degree of divergence, both within and across species. These receptors are ligand-gated ion channels, composed of seven novel transmembrane domains with an inverted membrane topology, compared with mammalian ORs24,25,26. ORs are expressed in ORNs and can receive a variety of volatile chemicals, including pheromones and general odorants14, 19. The function of an insect OR depends on the presence of a non-ligand binding odorant receptor co-receptor (Orco), which functions as a ligand-gated ion channel27,28,29. In contrast to ORs, Orco is highly conserved across insect species.

After the annotation of ORs in D. melanogaster genome, GRs, a common ancestor to ORs and composed of seven transmembrane domains, were discovered in D. melanogaster 30. In insect, GRs are also conserved in their sequence and structure. They are highly expressed in the GRNs in taste organs31, 32. For instance, GRs have been shown to play a critical role in coordinating insect feeding behaviors. GRs located on the dendrites of taste sensilla recognize the taste stimuli from the environment, especially in foods32, 33. In addition, GRs are also involved in the detection of carbon dioxide34. Based on the functional research data, GRs have been classified into four clades: CO2, GR43a-like, sugar, and bitter19, 35.

A large number of ORNs express neither ORs nor GRs, but they express IRs, which are also ligand-gated ion channels, but with three transmembrane domains36,37,38. IRs have been identified across Protostomia and are regarded as an ancient family of chemosensory receptors. IRs in insects can be classified into two types: the “antennal IRs,” which are conserved across insect orders with chemosensory function, and the “divergent IRs,” which are species-specific and are assigned a tentative role in taste37, 38. Meanwhile, two IRs, IR8a and IR25a, appear to act as co-receptors with the function of turning IRs sensory cilia targeting and IR-based sensory channels37,38,39,40.

A. gifuensis has been selected as a potential biological-control agent for the green peach aphid Myzus persicae Sulzer, one of the most common pests of several crops in China and Japan; it has already been successfully used to control M. persicae on tobacco in Yunnan and many other regions of China41,42,43. During the predation and parasitism, natural enemies utilize herbivore-induced volatiles (HIPVs), green leaf volatiles (GLVs), or the body volatiles such as aphid alarm pheromone E-beta-farnesene (EBF) to locate its hosts44,45,46. For example, A. gifuensis is able to discriminate the healthy, mechanically damaged or infested by its original aphid44. Furthermore, both female and male of A. gifuensis represented a positive electroantennogram (EAG) response to EBF and several plant volatiles, such as linalool, cis-3-hexen-1-ol, (E)-2-hexenal45. And a lot of natural enemies such as Aphidius ervi, Aphidius uzbekistanicus, Adalia bipunctata show attractant behavior to EBF47. All of these results revealed that natural enemies including A. gifuensis have evolved a comprehensive chemosensory system to enhance their parasitism efficiently.

However, the previous research on A. gifuensis only focused on its ecological behavior and anatomy. The potential molecular mechanism involved in the ecological process is lacking. Until now, the chemosensory receptors of parasitoid wasps have been only characterized in Microplitis mediator 48, Nasonia vitripennis 49, Macrocentrus cingulum 10, and Chouioia cunea 15. Therefore, we analyzed the chemosensory receptors of A. gifuensis based on sex-specific antennal transcriptomes through next-generation sequencing technology. The comprehensive identification and expression profile of chemosensory receptor genes in A. gifuensis provide valuable information for understanding the intraspecific or interspecific chemical communications, which is crucial for potential genetic manipulation of their sensitivity to chemical cues from hosts, plants, and themselves in biological control systems.

Results

Transcriptome assembly summary

The male and female A. gifuensis antennal transcriptomes were generated using Illumina Hiseq2000. Collectively, there were 38,848 transcripts, and the longest transcript was 13,876 bp in length. We identified a total of 19,074 components, each of which contained at least one annotated gene. The N50 transcript length was 1,980 bp and the total length of the assembled transcriptome was about 45.75 Mbp (Table S1).

Functional annotation

Functional annotations for the assembled database of A. gifuensis transcriptome were generated through diverse protein datasets. A total of 29,302 unigenes were annotated: 26,969 (92.0%) in NR_Annotation, 21,411 (7.30%) in Nt_Annotation, 11,086 (37.8%) in COG_Annotation, 22,259 (78.0%) in Swiss-prot_Annotation, 12,552 (42.8%) in GO_Annotation, and 20,319 (69.3%) in KEGG_Annotation (Table S2). From the database in NR_Annotation, 16,830 (62.4%) had a strong match with an e-value less than 1e−45 (Fig. 1A). For the database in NR_Annotation, 1,278 (4.74%) showed a strong similarity (95–100%) to known proteins (Fig. 1B). Approximately 70% sequences matched to a hymenopteran sequence (Fig. 1C).

Figure 1
figure 1

Homology analyses results. The BLASTx annotations of Aphidius gifuensis antenna transcripts (A) E-value distribution, (B) Similarity distribution, and (C) Species distribution.

All the annotated unigenes were classified into three groups: biological process, cellular components, and molecular functions. In the biological process, the most represented biological processes were cellular process (8,294 antennae unigenes) and single-organism process (6,410 antennae unigenes). In the cellular components, the genes expressed in the antennae were mostly cell part- (6,000 antennae unigenes) and organelle-related (3,974 antennae unigenes). In the molecular functions, binding (6,226 antennae unigenes) and catalytic activity (6,014 antennae unigenes) were the highly expressed categories in antennae (Fig. 2). In total, 11,086 of the 29,302 unigenes with non-redundant database hits were grouped into 25 COG categories (Figure S1).

Figure 2
figure 2

Functional annotation of Aphidius gifuensis antenna transcripts based on gene ontology (GO) categorization. GO analysis was performed at the level of two or three main categories (cellular component, molecular function, and biological process).

Identification of chemosensory receptors

Odorant receptors

Sixty-two candidate ORs were identified (Table 1 and Fig. 3). Only one of the transcripts was an incomplete fragment, whereas all the other transcripts represented a full-length gene, containing complete open reading frames (ORF). The transcript name, length, best Blast P, e-value, and identity are presented in Table 1.

Table 1 Candidate odorant receptor transcripts identified in adult male and female A. gifuensis antennal transcriptomes.
Figure 3
figure 3

The number of chemosensory genes in different insect species. The digits by the histogram bars represent the numbers of chemosensory receptor genes in different subfamilies. A phylogenetic tree showing the phylogenetic relationships between these species is illustrated on the left. The data are obtained from the current study for Microplitis mediator, Acyrthosiphon pisum, Apis mellifera, Nasonia vitripennis, Solenopsis invicta, Bombyx mori, Tribolium castaneum, Drosophila melanogaster, Anopheles gambiae, Aedes aegypti, and Culex quinquefasciatus.

The odorant co-receptor in A. gifuensis was identified as having an intact open reading frame with seven transmembrane domains. With the exception of Orco, only 14 of the 62 ORs showed more than 50% identity with known ORs in the NCBI database (Table 1). The phylogenetic analysis of A. gifuensis ORs is presented in Fig. 4, which includes the identified ORs from D. melanogaster, A. mellifera, N. vitripennis, and M. mediator. The amino acid sequences for all used ORs are listed in Table S5. In the phylogenetic tree, 16 A. gifuensis ORs (AgifORs) (OR3, OR30, OR31, OR32, OR33, OR34, OR35, OR36, OR37, OR38, OR39, OR40, OR41, OR42, OR43, and OR44) clustered in a species-specific subgroup, while the other AgifORs grouped with the ORs of other species. The expression of the sex-specific AgifORs with the transcriptome data-based heat map is shown in Figure S2. Orco1, OR4, OR9, OR17, OR18, OR19, OR24, OR25, OR26, OR27, OR28, OR29, OR33, OR39, and OR49 were highly expressed in both female and male antennae.

Figure 4
figure 4

Maximum likelihood phylogenetic tree of odorant receptors (ORs). Included are ORs from Aphidius gifuensis (Agif), Microplitis mediator (Mmed), Apis mellifera (Amel), and Nasonia vitripennis (Nvit).

Gustatory receptors

We identified 15 candidate GRs in the A. gifuensis antennal transcriptomes (Table 2 and Fig. 3). All these candidate GRs were identified with an intact open reading frame. A phylogenetic tree was constructed with sequences from A. gifuensis, N. vitripennis, A. mellifera, and D. melanogaster (Fig. 5). Five GRs (AgifGR1, AgifGR3, AgifGR4, AgifGR5, and AgifGR6) were found in a clade with sugar receptors, which included GRs identified from N. vitripennis, A. mellifera, and D. melanogaster. The sex-specific expression of GRs can be seen in the phylogeny of all A. gifuensis GRs with the transcriptome data-based heat map (Figure S3). The expression profiles of these GRs were diverse.

Table 2 Candidate gustatory receptor transcripts identified in adult male and female A. gifuensis antennal transcriptomes.
Figure 5
figure 5

Maximum likelihood phylogenetic tree of gustatory receptors (GRs). Included are ORs from Aphidius gifuensis (Agif), Microplitis mediator (Mmed), Apis mellifera (Amel), Drosophila melanogaster (Dmel), and Nasonia vitripennis (Nvit).

Ionotropic receptors

The identified candidate IRs in the A. gifuensis antennal transcriptome are listed in Table 3, with the best blast results from the NCBI database. In the phylogenetic tree, all IRs were classified into five clades, including antennal IRs (IR1, IR4), IR8a (IR5, IR9), IR25a, IR75u (IR2, IR7), and divergent IRs (Fig. 6). The most highly expressed IR transcripts in both male and female antennae were IR8a.1, IR8a.2, IR25s, and Nmdar1 (Figure S4). IR5, IR7, IR6, and IR9 showed significant sex-specific expression patterns.

Table 3 Candidate ionotropic receptor transcripts identified in adult male and female A. gifuensis antennal transcriptomes.
Figure 6
figure 6

Maximum likelihood phylogenetic tree of ionotropic receptors (IRs). Included are ORs from Aphidius gifuensis (Agif), Microplitis mediator (Mmed), Apis mellifera (Amel), Drosophila melanogaster (Dmel), and Nasonia vitripennis (Nvit).

Tissue- and host-specific expression profile of candidate A. gifuensis chemosensory receptors

In order to evaluate the heat map results of the chemosensory receptors and define the expression pattern of the identified genes, the expression profile of 3 ORs, 6 GRs, and 8 IRs in different tissues and hosts were analyzed using qRT-PCR, considering their sex-specificity (Figs 7 and 8). Furthermore, the tissue- and host specific will help us to have an initial functional prediction of these chemosensory genes.

Figure 7
figure 7

The tissue-specific transcript abundances of Aphidius gifuensis chemosensory receptor genes. FA: female antennae, MA: male antennae, FB: female body, MB: male body. Each bar graph in Fig. 7 represent a relative expression patterns of a chemoreceptor gene individually without any comparisons across panel. The error bars represents standard errors and the small letters above each bar indicate significant differences in transcript abundances (p < 0.05).

Figure 8
figure 8

The host-specific transcript abundances of Aphidius gifuensis chemosensory receptor genes. M. persicae: A. gifuensis reared on the green peach aphid, Myzus persicae, for at least 1 year; S. avenae: A. gifuensis reared on the English grain aphid, Sitobion avenae, for at least 1 year; A. pisum: the A. gifuensis reared on the pea aphid Acyrthosiphon pisum, for at least 1 year. Each bar graph in Fig. 8 represent a relative expression patterns of a chemoreceptor gene individually without any comparisons across panel. The error bars represents standard errors and the small letters above each bar indicate significant differences in transcript abundances (p < 0.05).

All these selected target genes were successfully detected. Out of these, OR18, OR28, GR1, GR3, GR5, GR10, IR8a.1, IR8a.2, IR3, and IR6 showed a ubiquitous expression pattern in female and male tissues. Orco1, IR25a, and Nmdar1 were found to be significantly expressed in antennae, especially in females. On the contrary, GR4, GR5, and IR3 were highly expressed in the body. GR6 and GR10 showed a sex- and tissue-specific expression profile.

The English grain aphid, Sitobin avenae, pea aphid, Acyrthosiphon pisum, green peach aphid, Myzus persicae are the common natural pest aphids of A. gifuensis in China. The expression patterns of A. gifuensis reared on different aphid species were measured by RT-qPCR. Orco1, OR18, OR28, GR1, GR4, GR5, GR6, GR10, IR8a.1, IR8a.2, IR25a, and Nmdar1 were highly expressed in the female A. gifuensis reared on Sitobion avenae and A. pisum, compared with the A. gifuensis reared on M. persicae. Meanwhile, there were no significant differences in the expressions of OR28, GR3, GR4, GR6, GR10, IR8a.2, IR25a, and IR3 in the males. Compared with the A. gifuensis reared on M. persicae and A. pisum, GR3 and Nmdar1 were highly expressed in S. avenae-reared A. gifuensis.

Discussion

In this study, we identified 62 ORs, 15 GRs, and 23 IRs in the antennal transcriptomes of A. gifuensis. Since A. gifuensis is a key biological control agent, the identified chemosensory receptors represent a valuable genomic resource at the molecular level, for aphid-plant-parasitoid interactions.

In our transcriptome, 62 ORs were identified including one odorant co-receptor. The number of identified ORs in A. gifuensis is less than that in A. mellifera and N. vitripennis, which have a total of 170, 68, and 301 ORs, respectively21, 49. There could be several reasons for this difference. As OR expression is amenable to modulation by scent conditioning, and the laboratory-reared A. gifuensis have had no opportunity of exposure to the diverse variety of volatiles emitted from different plants and animals, some of the olfactory genes might not be well expressed. For example, we found that Orco1, OR18, and OR28 were highly expressed in the S. avenae- and A. pisum-reared female A. gifuensis, compared with the A. gifuensis reared on M. persicae. In addition, the physiological condition of the parasitoids can also affect the expression of their chemosensory receptor genes. A previous study reported that after blood feeding, the expression of OR1 in Anopheles gambiae was significantly decreased50. This revealed that maybe the expression of some chemosensory genes were too low to be detected by transcriptome under specific physiological conditions. Meanwhile, the sequenced tissue is another restricting factor. For example, in M. cingulum, McinGR2 and McinIR7e3 are specifically expressed in other tissues, such as legs, head with mouth parts and body tissues10. As shown in Table S4, the number of chemosensory receptors identified based on the transcriptome is lower than that identified based on the genome. Overall, the low number of ORs identified in A. gifuensis antennal transcriptome might result from the species difference, rearing conditions, sequenced tissue, sequencing depth, and other factors.

The OR27 was highly expressed in the whole body, indicating that it not only reacts with host odors but also plays other roles in the non-olfactory organs of A. gifuensis. For example, in the migratory locust, Locusta migratoria L., 11 conventional ORs, which are perceived as contacting pheromones, are highly expressed in non-olfactory tissues such as wings and legs17.

It has been found that Orco was responsible for adopting of the correct structure by OR, and worked as a selective ion channel during olfactory signal transduction19, 51. In Dendroctonus armandi, the silencing of Orco led to EAG declining to 11 major volatiles of its host52. In Aedes albopictus, the Orco gene was found to be crucial for transmitting olfactory signals and conventional ORs that contribute directly to odorant recognition29. RNA interference and behavioral assays in Locusta migratoria L. indicated that OR-based signaling pathways mediate their attraction to aggregation pheromones17. In ant, ORs are candidate CHCs receptors and orco co-receptor antagonist blocks CHC detection, which are the main social communication cues in ant colonies13, 53. Similar to this, CHCs have been demonstrated to be involved in the discrimination of aphid species by parasitoids and the regulation of parasitism strategy54. For example, A. gifuensis performed different on its original host aphids and the other aphid species55. It laid more eggs in the new introduced aphid species than its original host aphid to improve the success rate of parasitism. In our results, ORco1, OR18 and OR28 highly expressed in S. avenae and A. pisum reared parasitoids, whereas the expression of OR28 in A. pisum reared parasitoids is lower than that in S. avenae reared parasitoids. We hypothesize that the different expression profiles among different clones might be mainly resulted from the different information cues in different aphid species. Furthermore, in Nasonia, female CHC profiles can be perceived as sexual cues to attract males56. All of these results revealed that ORs in A. gifuensis might be not only involved in its intraspecific or interspecific chemical communications.

In natural ecology, parasitoids are obligate consumers of plant-derived foods, including carbohydrate - rich solutions such as nectar and homopteran honeydew, which has been demonstrated to be an information chemical of aphid parasitoids to locate their host aphids47, 57, 58. For example, the honeydew is mainly containing amino acids and several carbohydrates including sucrose, glucose, trehalose, erlose, fructose, maltose and maltotriose58, 59. In previous work, HarmGR4 expressed in female antennae was sensitive to D-fructose in Helicoverpa armigera 33. And the behavioral and electrophysiological experiments have also found the antennae to be involved in the perception of D-fructose. In addition, DmelGr64a expressed in GRNs was found to be required for the behavioral responses to glucose, sucrose, and maltose in D. melanogaster 60. More interestingly, taste receptors such as GR43a, GR64a, GR32a and GR28a expressed in Drosophila wing respond to sweet and bitter stimuli such as glucose and denatonium61. In the present work, the homology genes involved in sugars perception have been identified including GR1, GR3, GR4, GR5, GR6, GR13, and GR14 and one GR (GR10) was classified as fructose receptor. Meanwhile, GR1 and GR3 were expressed predominantly in the antenna, whereas GR4 and GR5 highly expressed in body. When reared on different aphid species, the diverse expression patterns were shown. GR4, GR5, GR6 and GR10 expressed highly in S.avenae and A. pisum reared female parasitoids whereas no difference was found in that of male except GR5. The different expression patterns between male and female might related with their different food types. The male parasitoid wasps mainly feed on pollen and nectar, whereas the female parasitoid wasps can also consume honeydew and the body fluid of the aphids57, 62, 63. All these results suggest that sugar related GRs in A. gifuensis might be involved in the discrimination of the honeydew and nutrition quality of pest aphids by antenna and ovipositor contact. And the further researches about the functions of GR1 and GR3 in the perception of sugars and behaviors regulations during the feeding and parasitism are needed done.

The number of identified IRs in this study is greater than those for other species. For example, there are 12 IRs in A. mellifera 21, 12 in N. vitripennis 49, 11 in M. mediator 48, and 13 in M. cingulum 10. Similar to N. vitripennis (two candidate IR25a orthologs), two candidate IR8a orthologs in A. gifuensis were identified, with 75% and 67% amino acid sequence identity. However, gene duplication for IR25a has not been detected. As the gene duplication of IR25a might be unique to some of the hymenopteran species including N. vitripennis and M. mediator, further research on the loss of IR25a duplication is needed. As with the relatively high antennal expression of the OR co-receptor 1 (Orco1), the most highly expressed IR transcripts in both male and female antennae were the putative IR co-receptors, IR8a, IR25a and IR76b, in addition to IR21a, which along with IR25a, seems to be involved in the detection of small changes of temperature38, 64, 65. In D. melanogaster, IR25a is expressed in different populations of sensory neurons, including those in the antenna and labellum and acts as a co-receptor with different odour-sensing IRs38. In this work, the highly expressed IR25a in antennae indicated that it might have a similar function in the antenna of A. gifuensis. Besides these IR co-receptors, another conserved IR has been identified is IR41, which along with IR64a and IR76b are considered to play vita roles in amine sensing66. Furthermore, IR25a, IR93a and IR40 of D. melanogaster have been demonstrated to participate the humidity preference behavior regulation mechanism67. All of these results revealed that IRs in D. melanogaster with a diverse roles in the interactions between D. melanogaster and environment. In this work, the homology genes of IR40a and IR76b were identified and named as AgifIR14, AgifIR16, AigIR6 and AgifIR8. However, due to the lacking investigation of IRs in the other insects, we only hypothesized that IRs in A. gifuensis might be involved into the similar functions with their homology genes in D. melanogaster.

In conclusion, the main purpose of this work was to identify the chemosensory receptors in A. gifuensis. And RT-qPCR of some selected genes were done to reveal an initial functional predication, which were supported by the functional investigation of their homology genes in other insects. Our results not only lay a solid foundation on the further investigation about the functions of these identified genes in A. gifuensis such as the CHCs discrimination, odor and sugar perceptions but also provide valuable information for understanding and investigating the intraspecific or interspecific chemical communications in the solitary parasitic wasps.

Materials and Methods

Insects rearing

A. gifuensis were collected from the pea aphid, Acyrthosiphon pisum Harris, which were reared on alfalfa. A laboratory colony was established and maintained at 21 °C with a 16 h light: 8 h dark photoperiod on A. pisum that were reared on broad bean (Vicia faba L., var. ‘Jingxuancandou’, Jinnong, Taigu, Shanxi, China).

For providing a host-specific experience, the A. gifuensis were reared on A. pisum, the green peach aphid Myzus persicae, and the English grain aphid Sitobion avenae, for at least one year.

RNA sequencing

Antennae of A. gifuensis were cut from newly emerged adult male or female wasps (1–2 days old) respectively, and were frozen in liquid nitrogen. This collection of antennae without any other tissues was immediately stored at −80 °C for further analysis. Total RNA was extracted from four hundred antennae of each sex for each replicate using TRIzol reagent (Takara Bio, Tokyo, Japan), as per manufacturer’s instructions. And there were three biological replicates for each sex. The RNA integrity was verified by 1% agarose gel electrophoresis and the quantity was assessed using a Nanodrop ND-2000 spectrophotometer. Synthesis of cDNA and Illumina library generation was completed at Beijing Genomics Institute (BGI) (Shenzhen, Shenzhen, Guangdong, China), using Illumina HiSeqTM2000 sequencing.

De novo Assembly and Gene Annotation

Transcriptome de novo assembly was carried out using a short reads assembling program—Trinity, which combines three independent software modules: Inchworm, Chrysalis, and Butterfly, to overcome the quality and polymorphism issues. In order to get comprehensive information about the genes, we aligned the unigenes larger than 150 bp to Nr, Nt, KEGG, Swiss-Prot, and COG databases, with e-value < 10−5. With Nr annotation, we used the Blast2GO program to get GO annotation of Unigenes. Next, the WEGO software was used to perform GO functional classification for all unigenes.

The unigene expression levels were calculated by fragments per kb per million reads (FPKM) method, using the formula, FPKM (A) = 103 (106 C)/NL. FPKM (A) was set as the expression level of Unigene A, and C was the number of fragments that uniquely aligned to Unigene A, N was the total number of fragments that uniquely aligned to all Unigenes, and L was the base number in the CDS of Unigene A. The FPKM method is able to eliminate the influence of different gene length and sequencing level on the calculation of gene expression. Therefore, the calculated gene expression can be directly used for comparing the differences in gene expression across samples.

Phylogenetic analysis of candidate chemosensory receptors

Amino acid sequences of the candidate ORs, GRs, or IRs were aligned using MAFFT, with FFT-NS-I iterative refinement method with JTT200 scoring matrix, unalignlevel 0.3, “leave gappy regions” set, and other default parameters. Bioedit Sequence Alignment Editor 7.1.3.0 (Ibis Pharmaceuticals, Inc., Carlsbad, CA, USA) was used for further manual editing. Phylogenetic trees were subsequently constructed by the Maximum likelihood (ML) method using PhyML3.1, based on the best-fit model LG + G estimated by ProtTest2.4. SH-like approximate likelihood ratio (aLRT-SH) supports were used to evaluate the reliability of internal branches. The trees were further edited using the ITOL tool. The identity scores of alignment were extracted using BioEdit software, and the heat map was constructed by ITOL based on a three-color scale. Phylogenetic trees were based on hymenopteran data sets. The OR data set contained 62 amino acid sequences from A. gifuensis, together with N. vitripennis (67), M. mediator (51), and A. mellifera (68). The GR dataset contained 6 amino acid sequences from A. gifuensis, together with sequences from N. vitripennis (67), A. mellifera (68), and D. melanogaster (95). The IR data set contained 9 A. gifuensis amino acid sequences, along with M. mediator (51), N. vitripennis (67), A. mellifera (68), and D. melanogaster (95) IR sequences. All amino acid sequences for the chemosensory receptors used in this study are shown in Table S5.

Chemosensory receptors in different insect species

The species phylogenetic tree was constructed based on the alignment results of cytochrome oxidase subunit I (COI) genes from different species, using Mega 6. The trees were further edited using the ITOL tool with the number of identified chemosensory receptors. The number of identified chemosensory genes in all insects is shown in Table S4. The number of chemosensory-related genes was collected from published papers. The GenBank numbers of COI are Microplitis mediator (GenBank ID: KJ459149.1), Acyrthosiphon pisum (GenBank ID: AB506720.1), Apis mellifera (GenBank ID: AY114465.1), Nasonia vitripennis (GenBank ID: EU746554.1), Solenopsis invicta (GenBank ID: JN808838.1), Bombyx mori (GenBank ID: EU141360.1), Tribolium castaneum (GenBank ID: KJ003352.1), Drosophila melanogaster (GenBank ID: KJ767244.1), Anopheles gambiae (GenBank ID: DQ465336.1), Aedes aegypti (GenBank ID: GQ165783.1), Aphidius gifuensis (GenBank ID: GU097658.1), and Culex quinquefasciatus (GenBank ID: GQ165766.1)

Expression analysis

Heatmap plots were generated for the binary logarithm of raw FPKM-plus 1 values. For each plot, the minimum value was set to the number type, with a value of zero, and displayed as yellow, the midpoint was set to percentile type, with a value of 100, and displayed as blue, and the maximum was set to the highest value type, and displayed as red. These plots were made and edited using ITOL tool.

Quantitative reverse transcription PCR was performed to validate the expression of candidate chemosensory receptors in A. gifuensis. The collection of antennae and body tissues without antennae of each sex were collected respectively (antennae: 400 of each sex; body tissues: 20 of each sex) and were frozen in liquid nitrogen. For the host aphid specific expression analysis, A. gifuensis reared on different aphids were collected (whole body and 20 of each sex) and frozen in liquid nitrogen. Total RNA of A. gifuensis was extracted using TRIzol reagent (Takara Bio, Tokyo, Japan), as per manufacturer’s instructions. The temple RNA was treated using Dnase I and incubated at 42 °C for 2 min to remove the genomic DNA. Next, the cDNA was synthesized from total RNA using Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics, Mannheim, Germany) according to the standard manufacturer’s protocol. Gene-specific primers were designed by Primer Premier 5 (PREMIER Biosoft International, Palo Alto, CA, USA), and are shown in Table S3. qPCR was conducted in 20 μl reactions containing 50 × SYBR Premix, Ex Taq (10 μL), primer (10 mM), sample cDNA (0.8 μL), and sterilized ultra-pure grade H2O (7.6 μL). Cycling conditions were 95 °C for 30 s, 40 cycles of 95 °C for 5 s, and 55 °C for 30 s. Each sample had three technical replicates and three biological replicates. Relative quantification was performed using the Comparative 2−ΔΔCT method. Transcription levels of these receptor genes were normalized by 18 S RNA, and the normalization of each gene was compared with the lowest expression level in different tissues68. The expression data among the different tissues and host aphids of each sex were subjected to one-way analysis of variance (ANOVA); means were separated using Duncan’s test at P  < 0.05.