Identification and expression analysis of chemosensory receptor genes in an aphid endoparasitoid Aphidius gifuensis

Olfaction and gustation play critical roles during the host-location search process of insects. Several chemosensory receptor genes are thought to be involved in providing specificity to the olfactory sensory neuron responses. The aphid endoparasitoid, Aphidius gifuensis, has been used as a biological control agent against a variety of aphid species; this parasitoid is able to detect its target host(s) effectively during the parasitic process. To understand the mechanism of host detection in A. gifuensis, we assembled specific antennal transcriptomes of each sex through next generation sequencing technology to identify the major chemosensory receptor genes. Using a bioinformatics screen, we identified 100 olfactory receptors candidates (62 odorant receptors, 15 gustatory receptors, and 23 ionotropic receptors) from the sex-specific antennal transcriptome. In addition, combining with the demonstrated functions of chemosensory genes in other insects, the sex-, tissue-, and host-specific expression profile of chemosensory genes potentially revealed the candidate physiological functions. The identification and expression profile of chemosensory receptor genes in A. gifuensis provide valuable information for understanding and investigating the intraspecific or interspecific chemical communications in the solitary parasitic wasps.


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).
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).
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.
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,   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.
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.
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. 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.

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, respectively 21,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 decreased 50 . 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 tissues 10 . 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 Scientific RepoRts | 7: 3939 | DOI:10.1038/s41598-017-03988-z migratoria L., 11 conventional ORs, which are perceived as contacting pheromones, are highly expressed in non-olfactory tissues such as wings and legs 17 .
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 transduction 19,51 . In Dendroctonus armandi, the silencing of Orco led to EAG declining to 11 major volatiles of its host 52 . In Aedes albopictus, the Orco gene was found to be crucial for transmitting olfactory signals and conventional ORs that contribute directly to odorant recognition 29 . RNA interference and behavioral assays in Locusta migratoria L. indicated that OR-based signaling pathways mediate  their attraction to aggregation pheromones 17 . In ant, ORs are candidate CHCs receptors and orco co-receptor antagonist blocks CHC detection, which are the main social communication cues in ant colonies 13,53 . Similar to this, CHCs have been demonstrated to be involved in the discrimination of aphid species by parasitoids and the regulation of parasitism strategy 54 . For example, A. gifuensis performed different on its original host aphids and the other aphid species 55 . 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 males 56 . 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 aphids 47,57,58 . For example, the honeydew is mainly containing amino acids and several carbohydrates including sucrose, glucose, trehalose, erlose, fructose, maltose and maltotriose 58,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 denatonium 61 . 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 aphids 57, 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 temperature 38,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 IRs 38 . 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 sensing 66 . Furthermore, IR25a, IR93a and IR40 of D. melanogaster have been demonstrated to participate the humidity preference behavior regulation mechanism 67 . 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 HiSeq TM 2000 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) = 10 3 (10 6 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.
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 H 2 O (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 tissues 68 . 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.