Antennal transcriptome analysis of the piercing moth Oraesia emarginata (Lepidoptera: Noctuidae)

The piercing fruit moth Oraesia emarginata is an economically significant pest; however, our understanding of its olfactory mechanisms in infestation is limited. The present study conducted antennal transcriptome analysis of olfactory genes using real-time quantitative reverse transcription PCR analysis (RT-qPCR). We identified a total of 104 candidate chemosensory genes from several gene families, including 35 olfactory receptors (ORs), 41 odorant-binding proteins, 20 chemosensory proteins, 6 ionotropic receptors, and 2 sensory neuron membrane proteins. Seven candidate pheromone receptors (PRs) and 3 candidate pheromone-binding proteins (PBPs) for sex pheromone recognition were found. OemaOR29 and OemaPBP1 had the highest fragments per kb per million fragments (FPKM) values in all ORs and OBPs, respectively. Eighteen olfactory genes were upregulated in females, including 5 candidate PRs, and 20 olfactory genes were upregulated in males, including 2 candidate PRs (OemaOR29 and 4) and 2 PBPs (OemaPBP1 and 3). These genes may have roles in mediating sex-specific behaviors. Most candidate olfactory genes of sex pheromone recognition (except OemaOR29 and OemaPBP3) in O. emarginata were not clustered with those of studied noctuid species (type I pheromone). In addition, OemaOR29 was belonged to cluster PRIII, which comprise proteins that recognize type II pheromones instead of type I pheromones. The structure and function of olfactory genes that encode sex pheromones in O. emarginata might thus differ from those of other studied noctuids. The findings of the present study may help explain the molecular mechanism underlying olfaction and the evolution of olfactory genes encoding sex pheromones in O. emarginata.

Lepidoptera sex pheromones are divided into two main types based on their chemistry [16]. Type I pheromone components have 10-to 18-carbon, even numbered straight chain acetates, aldehydes, and alcohols. Type II pheromones consist of polyunsaturated C 17 -C 23 straight chains, skipped conjugated polyenic hydrocarbons and the corresponding epoxide derivatives [17]. Type I pheromones occur in about 75% of all studied moth species, whereas type II pheromones occur in about 15% of identified Lepidopteran pheromones [17]. These two major types of sex pheromones are produced through distinct pathways that involve different biosynthetic sites, substrates, and enzymes, as well as respectively employ specific endocrine regulatory mechanisms. However, both types of pheromones have the same function in mate recognition and attraction in moths [16,18].
The piercing fruit moth Oraesia emarginata Fabricius (Lepidoptera: Noctuidae) is an important pest of fruits such as citrus, pear, peach, and plum. The larvae feed on plants belonging to the Menispermaceae. Adult moths obtain nutrition from ripe fruits. Mated females lay eggs on Menispermaceae plants (Fig 1) [27]. The electroantennographic and behavioral responses of O. emarginata to volatiles from ripe fruits [28] and the repellency of a volatile compound, sec-butyl β-styryl ketone have been studied [29]. However, little is known about the olfactory mechanism of O. emarginata. Type II pheromones were identified as female sex pheromones in Oraesia species. The major and minor sex pheromone components of the related O. excavate were identified as cis-9,10-epoxy-(Z)-6 -heneicosene and cis-9,10-epoxy-(Z,Z)-3,6-heneicosadiene [30]. Although the sex pheromone of female O. emarginata was not published, it was similar to epoxide components from a preliminary identification (Du et al., unpublished data). In the present study, we achieved significant coverage of olfactory genes with de novo transcriptome and measured gene expression using real-time quantitative reverse transcription PCR analysis (RT-qPCR) for comparison between the sexes. We also discuss the diversification of olfactory genes for the recognition of type I and type II pheromones.

Insects
O. emarginata larvae were collected from fields in Gannan City of Jiangxi Province, China and reared in the laboratory at 25 ± 1˚C and 75 ± 5% relative humidity with a 14-h light/10-h dark photoperiod. Our field collection activities did not impact endangered or protected species. Larvae were fed fresh leaves of Cocculus orbiculatus until pupation. Emergence of males and females was checked every morning, and adults were separately maintained in ventilated wooden cages (35 cm × 35 cm × 50 cm). Emerging adult moths were fed with 10% glucose water soaked into cotton. 90-bp read length for the paired-end reads by BGI (Shenzhen, Guangdong, China). Dirty reads containing adapters and unknown or low-quality bases were discarded from the raw reads to obtain clean reads for analysis. De novo transcriptome assembly was conducted with the short reads assembly program, Trinity (r20140413p1, min_kmer_cov:2) [31]. BLASTx (v2.2.28+) alignment (E value < 0.00001) between unigenes and protein databases (NCBI non-redundant protein database, Swiss-Prot, Kyoto Encyclopedia of Genes and Genome (KEGG), and Clusters of Orthologous Groups (COG)) was successively performed. Gene ontology (GO) annotations of the unigenes were determined using Blast2go (http://www. blast2go.org/) [32].

Olfactory gene analysis
The candidate olfactory gene was manually obtained from gene annotation. In addition, a 50% ORF length cutoff was used in identifying putative genes to prevent a gene from being counted twice. The candidate OBPs and CSPs were searched for the presence of N-terminal signal peptides using SignalP4.0 (http://www.cbs.dtu.dk/services/SignalP/) using default parameters [33]. The signal peptides likely contained significant phylogenetic information and were included in the phylogenetic analyses of OBPs and CSPs [34]. Amino acid sequence alignment was performed using CLUSTALX2.1 using default parameters [35]. For phylogenetic analysis, known amino acid sequences of olfactory genes from other insects were downloaded (S1 File). Phylogenetic analyses were conducted using the maximum likelihood method of MEGA 6.0, which was based on the Jones-Taylor-Thornton (JTT) substitution model, partial deletion gaps with 95% site coverage cutoff, a nearest neighbor interchanges (NNI) heuristic search, and other default parameters [36]. Node support for the phylogenetic tree was assessed using the bootstrap method with 1,000 bootstrap replicates.

Profiling analysis of gene expression based on the antennal transcriptome
Gene expression levels were calculated using the fragments per kb per million fragments (FPKM) method based on the results of antennal transcriptome analysis. The number of fragments that uniquely aligned to a gene was divided by the total number of fragments that uniquely aligned to all genes and by the base number in the CDS of that gene [37]. The FPKM method can eliminate the influence of different gene lengths and sequencing levels on the calculation of gene expression.

RT-qPCR analysis of olfactory gene expression in the antennae
Single-stranded cDNAs were synthesized from 1 μg of total RNA using the ReverTra Ace qPCR RT Kit (Toyobo, Kita-ku, Osaka, Japan) following the manufacturer's recommendations. RT-qPCR was performed with SsoFast™ EvaGreen 1 Supermix (Bio-Rad, Hercules, CA, USA), following the manufacturer's protocols, in a CFX-96™ PCR Detection System (Bio-Rad). The cycling conditions were an initial cycle at 95˚C for 30 s, followed by 39 cycles of 95˚C for 5 s and 60˚C for 5 s. Dissociation curves with 0.3˚C/s melt rates were used to check for the presence of non-specific dsDNA SYBR Green hybrids. Only primers with a single PCR amplification product were used in the subsequent analyses. The amplification efficiency of each primer was calculated from the slope of the standard curve [38]. The PCR primers used are listed in S1 Table. Ubiquinol-cytochrome c reductase (UCCR) and arginine kinase (AK) were used as reference genes. The difference in gene expression was measured by using the 2 -ΔΔCt algorithm [39]. Differential gene expression between females and males was measured, with the female antennae used as reference. Expression levels of target genes were normalized independent of each reference gene with the algorithm, and then averaged. When the gene expression of the female antennae was very low, the gene expression of the male antennae was used as control. RNA extraction was repeated three times for each sample, and two or more RT-qPCR replicates were prepared for each sample.

Data analysis
Data analysis was conducted using SPSS 17.0 (SPSS Inc., Chicago, IL, USA). The significance of the difference between means was determined using the student's t-test. The critical P value for each test was set at 0.05.

De novo antennal transcriptome assembly
Using the Illumina HiSeq™ 2000 sequencing system, 117,410,034 raw reads were obtained from the antennal samples. After removing low-quality (< Q20) adaptor and contaminating sequence reads, 103,301,292 (a total of 9,297,116,280 bp) clean reads were generated from antennae, and 42,992 unigenes were assembled (N50 = 1,098), with a mean length of 713 bp. More than 58% (24,954) of the unigenes were aligned to sequences in various protein databases. GO annotation was performed to obtain information on their molecular function, biological process, and cellular location (S1 Fig). The raw sequence of the transcriptome has been deposited to the National Center for Biotechnology Information (NCBI) (GenBank Accession Number PRJNA358570; https://www.ncbi.nlm.nih.gov/bioproject/PRJNA358570).
A total of 20 candidate chemosensory protein (CSP) genes were identified in O. emarginata, with a mean length of 128 aa. The full ORF of the 16 CSP genes were obtained (Table 3, Fig 5). In the phylogenetic tree, OemaCSP9 and OemaCSP16 were clustered the homologous genes of other insect species into two conserved groups (Fig 5). The bootstrap values of 5 CSPs (OemaCSP1, 2, 7, 8, and 10) were < smaller than 50%, although these were clustered with studied CSPs of the Lepidopteran species. Four conserved cysteines were found in all CSP genes, but OemaCSP16 differed from the other CSPs in terms of the number of amino acids (Fig 6). Six candidate ionotropic receptor (IR) genes and 2 sensory neuron membrane protein (SNMP) genes were identified in O. emarginata, and their mean lengths were 535 aa and 522 aa, respectively (Tables 4 and 5). All O. emarginata IRs and SNMPs were clustered with Lepidopteran IRs and SNMPs, respectively, with the bootstrap values > 80% (Figs 7 and 8). The full ORF of 2 SNMP genes was obtained.

Expression of olfactory genes with RNA sequences
The FPKM values of the chemosensory receptors were < 60, and OemaORco showed the highest FPKM value (Tables 1 and 4). The FPKM value of OemaOR29 was higher, but those of the other candidate PRs were lower than the general ORs, including OemaOR14, 25, 27, and 32 ( Table 1). The FPKM values of OemaIR75p and OemaIR21a were larger than those of the coreceptors OemaIR25a and OemaIR8a (Table 4). In contrast to chemosensory receptors, 39.0% of the OBP and 52.4% of the CSP genes showed FPKM values > 300, including 3 candidate PBPs (Tables 2 and 3). OemaPBP1 showed the highest FPKM value among all OBPs, and OemaCSP19 had the highest FPKM value among all chemosensory genes. The FPKM value of OemaSNMP1 was < 20, but that of OemaSNMP2 was > 500 ( Table 5).

Phylogeny of pheromone recognition genes of types I and II pheromones
In the phylogenetic tree, 4 orthologous PRs clusters for type I pheromones were obtained (Cluster PRI-PRIV), and candidate PRs of the noctuid species (excluding O. emarginata) formed subclusters of these 4 clusters, with high bootstrap support (! 89, Fig 10). OemaOR29   and ObruOR1 (the only identified pheromone receptor for type II sex pheromones from the geometrid O. brumata) belonged to cluster PRIII (Fig 10). Other candidate PRs of O. emarginata were not grouped with any of these 4 clusters, but 5 (OemaOR3, 4, 21, 26, and 28) were clustered, with a bootstrap support of 78 (Fig 10). The PBPs and GOBPs of all test species were clustered into 3 (Cluster PBPI-PBPIII) and 2 (Cluster GOBPI-II) apparent clusters, with good bootstrap support (! 52), respectively ( Fig  11). OemaPBP3 and OemaGOBP1 were clustered with orthologous PBPs and GOBP1s of the other noctuids for type I pheromones, respectively (bootstrap support ! 56) (Fig 11). However, OemaPBP1, OemaPBP2, and OemaGOBP2 were not clustered within PBPs and GOBP2s from other noctuid species for type I pheromones. OemaPBP2 was clustered with MsexPBP2, with a bootstrap value of 74 (Fig 11).

Discussion
The unique life history of O. emarginata might have driven the increase in the number of chemosensory genes O. emarginata has a unique life history. The larvae feed on Menispermaceae plants, but adults suck on the juices of ripe fruits. Mating behavior is mediated by female sex pheromones. Mated females oviposit on Menispermaceae plants. Odorant classes from different species might thus be different [52]. Moths of O. emarginata must recognize a range of different odors with diverse chemical structures emitted from conspecifics, fruits, or orchard background and larval host plants. The olfactory acuity and discriminatory power in O. emarginata may have evolved to fulfill its ecological needs. We found 104 candidate olfactory genes in the antennae of O. emarginata, including 35 ORs, 41 OBPs, 20 CSPs, 6 IRs, and 2 SNMPs. In these 104 olfactory genes, 2 ORs (OemaOR24 and 35) and 5 CSPs (OemaCSP1, 2, 7, 8, and 10) were not effectively clustered with those of other Lepidopterans (bootstrap values < 50) in the phylogenetic analysis. In addition, 8 OemaORs (OemaOR11, 14, 17, 19, 20, 25, 27, and 32) were clustered into the clade of OfurOR34, MsexOR42, and AdisOR9 (bootstrap value = 87) (Fig 2), and 7 OemaOBPs (OemaOBP4, 11, 13, 18, 23, 27, and 35) were clustered with AipsOBP4, SlitABP1, SlitOBP12, SexiABP1, HvirABP2, HarmOBP7, and HarmOBP7.2 (bootstrap value = 61) in the phylogenetic trees (Fig 3). Some of those genes might be species-specific to O. emarginata and used to recognize the odors produced by the Menispermaceae and fruits. The number of chemosensory binding proteins (including OBPs and CSPs) was slightly smaller than in B. mori, which included the whole genome, but larger than in other moth species studied using the same protocol (antennal transcriptome). These other species included polyphagous insects such as S. litura ( Table 6). The larger number of chemosensory binding proteins might be due to the life history of O. emarginata and the larger database in our study. We found a total of 103,301,292 reads that were assembled into 2,202,660 contigs, and compared to 55,288,304 reads assembled into 105,971 contigs in S. litura [51]. However, the number of chemosensory receptors was lower than in most other moths ( Table 6). The low expression level of chemosensory receptor genes (FPKM < 60) and short read length (250 bp) of the transcriptome analysis might have resulted in short sequences for many chemosensory receptor genes. However, the long sequence of the chemosensory receptor genes (about 400 aa and 800 aa for OR and IR, respectively) [53,54] and the criterion of 50% ORF length cutoff might have excluded numerous chemosensory receptors with short sequences. No gustatory receptor gene was identified in the antennae, which suggests that the antennae of O. emarginata are not major taste organs. The proboscis, which harbors considerably fewer sensilla than antennae, are believed to specialize in taste reception in some moths [37,55]. In addition, the long sequence of gustatory receptor genes (about 400 aa) and the criterion of 50% ORF length cutoff might have excluded some gustatory receptors with short sequences.

Olfactory genes with sex-specific expression
We identified 2 candidate PRs (OemaOR29 and 4) and 2 candidate PBPs (OemaPBP1 and 3) that showed male-biased expression and might be involved with female sex pheromone recognition in O. emarginata. Our results were consistent with the study on the sex pheromone recognition in a sibling speciesm O. excavate, which produces two sex pheromone compounds at the ratio of 86: 14[30]. OemaOR29 was clustered with ObruOR1 and AsegOR3 in the phylogenetic tree, which recognized the pheromonal tetraene of O. brumata, 3Z,6Z,9Z-19:H and the  triene 3Z,6Z,9Z-21:H separately [56]. OemaPBP1 and OemaPBP3 were ranked in the clusters PBPI and PBPIII in the phylogenetic analysis, respectively, which showed an equally consistent association with male-specific pheromone sensitive sensilla [57]. Orthologous genes in the clusters PBPI and PBPIII play critical and minor roles in female sex pheromone perception, respectively [58][59][60][61]. OemaOR29 and OemaPBP1 showed the highest FPKM values in all ORs and OBPs, respectively, and might be used to recognize the main sex pheromone component.
OemaOR4 and OemaPBP3 might be involved in the recognition of the minor sex pheromone component. Further studies are needed to verify the function of these genes. Five candidate pheromone receptor genes (OemaOR3, 21, 26, 28, and 30) showed femalebiased expression, and OemaOR26, and OemaOR28 were specifically expressed in females. The function of these genes is unknown, but these might be used by females to recognize male pheromones. Production of short-range pheromones has been reported in male butterflies [62]; these function in female mate selection, act as an aphrodisiac, and arrest female departure [63,64].
Besides the candidate PR genes, some genes with sex-specific expression were detected; for example, OemaOR13 was female-specific. These genes might also be correlated with sex specific behaviors such as the recognition of oviposition cues by females [65][66][67].

Diversification of olfactory recognition to sex pheromones
Type II pheromones have mainly been found in the moth superfamilies Geometroidea and Noctuoidea [17], but olfactory genes for type II pheromones were only identified in the geometrids A. selenaria cretacea [68,69] and O. brumata [56] and the erebids L. dispar [70][71][72] and Hyphantria cunea [73]. The sex pheromone of female O. emarginata was not published, but it   was similar to the epoxide components of a preliminary identification (Du et al., unpublished data). In addition, cis-9,10-epoxy-(Z)-6 -heneicosene and cis-9,10-epoxy-(Z, Z)-3,6-heneicosadiene were identified as the major and minor sex pheromone components from a sibling species, O. excavate [30]. In the present study, 7 candidate PRs and 3 candidate PBPs were obtained from the noctuid O. emarginata using antennal transcriptome analysis.
The diversification of olfactory recognition to sex pheromones has been verified for type I pheromones in noctuids such as A. segetum, H. armigera, and S. litura, and the phylogeny of moth PRs and PBPs for type I pheromone identified several apparent orthologous clusters (cluster PRI-PRIV for PRs and cluster PBPI-PBPIII for PBPs). PRs and PBPs from different clusters specifically respond to different type I sex pheromone components [59,74]. Although the functions of PRs for type II pheromone recognition were not identified, phylogenetic analysis clustered 3 candidate PRs of H. cunea [73] and 7 candidate PRs of O. emarginata into three groups. These findings are indicative of the diversification in olfactory recognition to type II pheromones.
Phylogenetic analysis did not separate the PRs and PBPs for types I and II pheromones, thereby suggesting that PRs and PBPs for types I and II pheromones evolved from a common ancestor. However, type I pheromones differed from type II pheromones in its chemical characteristics. OemaOR29 and ObruOR1 belonged to cluster PRIII of type I pheromone recognition, which is under strong purifying selection (a very small dN/dS values), and did not respond to any type I sex pheromone components [75]. On the contrary, ObruOR1 was verified to specifically recognize the pheromonal tetraene of O. brumata, 3Z,6Z,9Z-19:H, and the orthologous receptor AsegOR3 responded strongly to the triene 3Z,6Z,9Z-21:H instead of any female sex pheromone of A. segetum [56]. Cluster III might be specialized in the recognition type II sex pheromone components. In addition, 6 other candidate PRs of O. emarginata were not grouped within any of the four PR clusters of type I sex pheromones, but 5 of these were grouped into a specific cluster, with a bootstrap support value of 78. The candidate main sex pheromone-binding protein OemaPBP1 was not clustered into the subgroup of PBP1 genes from other noctuid species in the phylogenetic tree. These results indicate that the olfactory genes for sex pheromones in O. emarginata might differ from those of other noctuid species,

Conclusions
A total of 104 candidate olfactory genes, including 7 candidate PRs and 3 candidate PBPs were identified from the noctuid O. emarginata. Seven olfactory genes of O. emarginata were not effectively clustered with those of other Lepidoptera, and OemaORs and OemaOBPs in 2 clusters were strongly expanded. These changes in olfactory genes in O. emarginata might correlate with its unique life history. Most candidate PRs and PBPs (except for OemaOR29 and OemaPBP3) of O. emarginata were not clustered with other noctuid species. OemaOR29 was grouped into cluster PRIII of type I pheromones, which recognized type II pheromones instead of type I pheromones. Noctuid species might thus have undergone diversification of the pheromone recognition gene for types I and II sex pheromones. Our results increase our understanding of the molecular mechanism of O. emarginata olfaction and the evolution of olfactory genes associated with sex pheromones.
(DOC) S1 File. Amino acid sequences of the olfactory genes used in the phylogenetic analysis. (TXT)