Male- and Female-Biased Gene Expression of Olfactory-Related Genes in the Antennae of Asian Corn Borer, Ostrinia furnacalis (Guenée) (Lepidoptera: Crambidae)

The Asian corn borer (ACB), Ostrinia furnacalis (Guenée), is a destructive pest insect of cultivated corn crops, for which antennal-expressed receptors are important to detect olfactory cues for mate attraction and oviposition. Few olfactory related genes were reported in ACB, so we sequenced and characterized the transcriptome of male and female O. furnacalis antennae. Non-normalized male and female O. furnacalis antennal cDNA libraries were sequenced on the Illumina HiSeq 2000 and assembled into a reference transcriptome. Functional gene annotations identified putative olfactory-related genes; 56 odorant receptors (ORs), 23 odorant binding proteins (OBPs), and 10 CSPs. RNA-seq estimates of gene expression respectively showed up- and down-regulation of 79 and 30 genes in female compared to male antennae, which included up-regulation of 8 ORs and 1 PBP gene in male antennae as well as 3 ORs in female antennae. Quantitative real-time RT-PCR analyses validated strong male antennal-biased expression of OfurOR3, 4, 6, 7, 8, 11, 12, 13 and 14 transcripts, whereas OfurOR17 and 18 were specially expressed in female antennae. Sex-biases gene expression described here provides important insight in gene functionalization, and provides candidate genes putatively involved in environmental perception, host plant attraction, and mate recognition.


Introduction
The olfactory and chemosensory systems of Lepidoptera are important for several biologicallyimportant functions including adult mate attraction, oviposition site selection and host plant preference, and negative taxis [1]. Trichoid sensilla are located on moth antennae and composed of pore tubes through which volatile hydrophobic odorants from the environment can enter. The specific detection and subsequent behavioral responses to environmental volatiles are mediated by the initial binding and transport of hydrocarbons across the aqueous sensillar lymph by classes of odorant binding proteins (OBPs) and chemosensory proteins (CSPs). The OBPs are small hydrophilic proteins that have a conserved tertiary protein structure of 6 alpha-helices coordinated by 3 disulfide bridges [2], and CSPs have 4 alpha-helices that form 2 disulfide bridges [3]. Both OBPs and CSPs are localized in the sensillar lymph of trichoid sensilla [4]. CSP sequences are comparatively more highly conserved, whereas OBPs have diverged and are classified into subfamilies based on functional and phylogenetic evidence: Classic, Dimer, D7, pheromone binding proteins (PBPs)/general odorant binding proteins (GOBPs), chemical-sense-related lipophilic-ligand-binding proteins (CRLBPs), antennal binding protein group I (ABPI) and II (ABPII), and atypical Plus-C and Minus-C OBPs [5,6]. OBPs can retain little intraspecific homology outside of six conserved cysteine residues, but those that are missing the C2 cysteine are referred to as the minus-C OBP class and those with an excess of cysteines belong to the plus-C OBP class. Both OBPs and CSPs perform analogous functions of bind chemical cues encountered in the environment and transporting these cues within chemosensory tissues to receptors located at the neuron surface.
The ligands that are bound by a majority of OBPs, including the lepidopteran PBP subfamily, remain largely unknown nor are the subsequent behavioral responses fully understood, with the possible exception of GOBPs and other OBPs which may function in host plant volatile recognition, taste [7] or xenobiotic perception [8]. The first described OBP was an Antheraea polyphemus PBP that bound radioactively labeled sex pheromones from conspecific females and was hypothesized to function in the perception of conspecific pheromones [4]. Indeed, PBPs were shown to bind sex pheromones in vitro and subsequently hypothesized to play a potential role in the discrimination of pheromone cues by male Lepidoptera [9,10]. These PBPs may be involved in the pH-dependent binding and transport of female sex pheromones to odorant receptors (ORs) on olfactory receptor neurons (ORNs). Specifically a pH-dependent conformational change was shown to shift the position of the C-terminal tail in or out of the hydrophobic PBP binding pocket from Bombyx mori [11], Amyelois transitella [12], and Antheraea polyphemus for PBPs [13], which may be a feature common of lepidopteran PBPs [14]. At neutral pH pheromones are predicted to bind a PBP hydrophobic binding pocket, which allows diffusion of hydrocarbons across the sensillar lymph and prevention of breakdown by pheromone-degrading enzymes [14]. Similarly, Große-Wilde et al. suggested that PBPs of B. mori can mediate the bombykol-induced activation of BmOR1 [10], and subsequently showed that OBPs promote pheromone sensitivity in a ligand-specific manner [15]. Directional selection between orthologs of male Ostrinia nubilalis and O. furnacalis antennal expressed PBP3 was hypothesized to result from the evolution of selective binding between structurally distinct sex pheromones emitted by corresponding females of the same species [16]. Furthermore, neuron response to the Drosophila pheromone, cis-vaccenyl acetate (cVA), was shown to directly depend on the function of an OBP called LUSH [17,18] that acts to solubilize and transport cVA [19]. Analogous OBP polymorphisms have also been shown to elicit variant behavioral responses in Drospohila [20,21]. Recent studies have reported electrophysiological response to the odorant indole decrease when the Anopheles gambiae AgamOBP1 was knocked down using RNAi [22]. However, it remains unclear how broadly these dependencies apply to other insect systems since seemingly contradictory studies have reported that silk moths can respond to the pheromone bombykol in absence of the cognate PBP [23,24].
Insect OBPs and CSPs perform analogous roles in that both reversibly binding small ligands with dissociation constants in the micromolar range, despite differences in structure [25] and biological function. Similar to the PBP subfamily of OBPs, CSP conformational changes are predicted when in association with cognate ligands [26]. Thus is has been hypothesized that CSPs may be involved in insect chemical communication, although most specific functions have not yet been discovered. CSPs are expressed in a variety of tissues and may be have evolved a divergent cellular functions involved in environmental perception [27]. For example, CSPs expressed in the chemosensory sensilla [27] are believed to detect environmental carbon dioxide levels [28] and modulate behavioral phase changes in the migratory locust [29]. The transduction of chemical cues from the environment to neurons likely involves a pathway analogous to OBPs, but the function of this portion of the insect chemosensory system also remains largely unknown.
Antennal-mediated olfactory detection in Lepidoptera involves ORNs that project into the sensillar lymph of trichoid sensilla, where specifically-tuning of each neurons is achieved by the expression of a specific OR. Each OR forms a voltage-gated ion channel following heterodimerization with OR2, which is also referred to as the odorant receptor co-receptor (Orco) and is an ortholog of the Drosophila melanogaster OR, DmOr83b. Species in the genus Ostrinia are a model for the study of the olfactory system of Lepidoptera, and has been used investigate the selectivity and response of ORs to sex pheromones in order to understand the general mechanisms of chemosensory response. Female O. nubilalis and O. furnacalis respectively produce and emit a blend of E/Z-11-tetradecenyl acetate (E/Z11-14:OAc) and E/Z-12-tetradecenyl acetate (E/Z12-14:OAc), which evoke responses by males of the corresponding species. The ORNs of male Ostrinia are classified with respect to strength of impulses produced in laboratory electophysiological recordings when stimulated by female pheromones; large, medium, and small spiking neurons. Large spiking ORNs in O. furnacalis respond to both E-and Z-12-14: OAc components of conspecific female pheromone blends, whereas medium spiking neurons responded with equal intensities to Z9-, E11-and Z11-14:OAc. Analogously, O. nubilalis large spiking ORNs specifically responded to intraspecific female pheromone components [30].
The molecular basis of these specific ORN responses have been partially elucidated in relation to species-specific male behavioral responses. Pheromone stimulation of co-expressed Orco and other OR proteins (ORx) in the Xenopus oocyte system are used to assay for specific responses by measuring changes in membrane ion permeability. These two-electrode voltage clamp electrophysiology measures indicate that O. nubilalis OR6 and OR2 (OnOR6/2) respond specifically to Z11-14:OAc, but OnOR3/2 and 5/2 respond to the known antagonist Z9-14:OAc as well and female O. nubilalis E11-, Z11-14:OAc and O. furnacalis emitted E12-and Z12-14: OAc [31]. Due to independent isolation, it should be noted that the nomenclature of Ostrinia ORs has the OfurOR4 described by Miura et al. [32] being the direct ortholog of OnOR3 and OfurOR3 from other studies [33]. More importantly, Xenopus oocyte assays showed that OnOR3 and OfurOR3 in complex with Orco produce species-specific electrophysiological responses when stimulated by corresponding female sex pheromones, and thus hypothesized to be the ORs expressed in large spiking ORNs [33]. Site-directed mutagenesis of amino acids at position 148 of OfurOR3 from an alanine to the serine in OnOR3 resulted in the electrophysiological response of the mutant OfurOR3 to O. nubilalis female pheromones, and linked the specific change to functional variation in species-specific pheromone responses [33]. The aforementioned proteins may interact with sensory neuron membrane proteins (SNMPs) and ionotropic receptors (IRs) for olfactory signal transduction, and odorant degrading enzymes (ODEs) [1,34] to restore the sensitivity of the sensory neuron [35].
Despite these advances in the elucidation of sex pheromone perception and ORN response, little is known regarding the specific function(s) of many OBPs, their cognate ORs, nor the behavioral responses these sensory pathways elicit. OBPs have been characterized from fully assembled genomes, and wherein they comprise a diverse gene family with complex evolutionary histories that likely may have resulted from a high degree of functional diversification [36]. For example, 44 OBPs are encoded by 6 tandem duplicated gene clusters in the B. mori genome indicating that paralogs may have evolved diverged functions involved in chemosensory detection [37], and have undergone an enigmatic path of gene gain and loss compared to other arthropods. Differential expression of genes involved in the lepidopteran olfactory system may be important for understanding the evolution of this duplicated gene family, and well as potentially uncovering the molecular basis for variation in moth response to potential mates, host plants and selection of oviposition sites. Sex biased expression of gene family members is an example of subfuncationalization [38], wherein duplicated genes may be retained in the genome due to the derivation of novel expression patterns and reinforced by exclusionary sex-specific functions. The maintenance of sex-biased gene expression may be under strong positive selection in instances where a gene function in sexual attraction and mating, where loss of this fidelity may negatively impact reproduction. For example, male antennal-specific Manduca sexta MsexOR-1 and MsexOR-4 are suggested to function in ORN response to sex pheromone [39]. Of the 48 OR genes identified in the B. mori genome by Wanner et al. [40], BmOR3 was expressed 6 to 8-times higher in females and 12 OR transcripts expressed predominantly in female antennae, and only 3 ORs were shown to be male antennae-specific [40]. Analogous genomic studies are yet to report system-wide sex-biased gene expression analyses in Lepidoptera, neither have comparative orthologies to B. mori olfactory genes been previously identified in another moth species.
Expressed sequence tag (EST)/transcriptome approaches have also been used to identify chemosensory receptors in arthropod species and present an alternative approach when full genome sequence assemblies are unavailable [41][42][43][44][45], and have proven to be valuable for elucidating olfactory system function. The motivation for this study was to use a transcriptomic approach to identify gene components of the Ostrinia antennal olfactory system (CSPs, OBPs and ORs), and apply sex-biased expression data to formulate hypotheses for future functional genomic research. Our prediction of orthologous gene relationships between O. furnacalis and B. mori provide a valuable tool for comparative genomic and function analyses. This study provides important tissue-and sex-biased gene expression data of olfactory-related genes in the antennae of O. furnacalis and define putative one-to-one orthologous gene relationships in Lepidoptera for comparative functional analyses. Results of this study are discussed in a system-wide evolutionary context wherein expression bias has precluded or reinforced the functional diversification of lepidopteran olfactory response pathways, and gene regulatory changes in multiple system components may act in concert to modulate sex-specific behaviors.

Insects rearing and antennae collection
Pupal O. furnacalis were obtained from a laboratory colony at the Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China. Male and female pupae were placed into different gauze cages for eclosion respectively. After emergence adults were fed with cotton dipped in a 10% honey solution, and the solution was renewed daily. Antennae from 3-day old adults of both sexes were sampled separately, pooled by sex and flash frozen in liquid nitrogen. Antennae, maxillary palpus and legs of both sexes were also collected and put into liquid nitrogen for downstream real-time quantitative RT-PCR (RT-qPCR) validation of gene expression.

Antennal transcriptome assembly and functional gene annotation
Total RNA from male and female antennae was extracted separately using TRIzol reagent (Invitrogen) according to manufacturer instructions. Subsequent cDNA library construction and Illumina sequencing was performed at SinoGenoMax Co., Ltd, Beijing, China. Briefly, enrichment of mRNA from~20 μg of total RNA used oligo (dT)25 magnetic beads, and 1 st strand synthesis was conducted using AMV reverse transcriptase primed with an oligo (dT) primer followed by 2 nd -strand cDNA synthesis by using random hexamer-priming for reactions including RNase H and DNA polymerase I. The cDNA was fragmented, end-repaired, and ligated with library-specific barcoded adaptors. Library fragments were PCR amplified a minimum number of cycles in order to avoid normalization for downstream quantitative gene expression analyses. Amplified products were purified with QIAGEN MiniElute PCR Purification Kit (Qiagen, Venlo, Netherlands), and approximately equal molar proportions of male and female indexed libraries were sequenced on a single flow cell of an Illumina HiSeq 2000. Each library was sequenced in a second technical replicate on an independent Illumina HiSeq 2000 lane with an approximate equal molar ratio of each library loaded. Raw sequence data from all runs were obtained in fastq format.

Assembly and functional gene annotation
Illumina output (fastq formatted read data) were trimmed for quality scores (q) < 20. A single de novo assembly of the combined reads was performed from trimmed read data using the short read assembler ABySS [46]. To assess assembler performance, assemblies of all the reads were compared for k values ranging from 26 to 50 bp [47]. CD-HIT-EST was used to clusters similar sequences and removed redundant segments using the web interface at http:// weizhong-lab.ucsd.edu/cdhit_suite/cgi-bin/index.cgi?cmd = cd-hit-est [48], and the longest contig sequences were retained as the best assembly result. The program Getorf from the EM-BOSS package [49] was used to predict and extract the longest open reading frames (ORFs) in assembled contig sequences. Since splice variants and sequence heterogeneity in UTRs tended to uncouple contigs belonging to the same gene, an all-versus-all tblastn search was performed (E-value cutoff 10 −50 , and protein identities 95%), and putatively homologous sequence data were aligned using the CLUSTAL W algorithm and inspected manually.
Functional gene annotations were collected for all contig (unigene) sequences 150 bp using Blast2GO [50], where initial searches of the National Center for Biotechnology Information (NCBI) non-redundant (nr) protein database were conducted with the BLASTx algorithm, followed by collection of gene ontology (GO) terms from the GO database and retrieval of KEGG Pathway designations. The BLASTx output was then processed with wapRNA [51]. Contigs receiving putative annotation as OR gene family members were investigated in greater detail. Assembly of contigs from iterative ABySS assemblies were achieved using CAP3 (default parameters; [52]), and were re-annotated by manual blastx of the NCBI protein database and search of conserved domain database (CDD) and protein family (pfam) databases (http:// pfam.xfam.org/search).

Analysis of differential gene expression
Differential gene expression and transcript bias between male and female antennal read data were conducted by independent alignment of short reads from male-and female-specific libraries to the reference antennal transcriptome assembly. Specifically, the Burrows-Wheeler Alignment (BWA [53]) was used to align reads to the reference antennal transcriptome ( 5 allowed mismatches), and expression level for each gene was initially calculated as Reads Per Kilobase per Million mapped reads (RPKM) using an in-house Perl module [51]. The RPKM method eliminates the influence of different gene lengths and sequencing discrepancies when calculations of expression abundance are made [54], and corrected variance in RPKM values between male and female alignments were normalized based on the depth of reads aligned to the housekeeping gene beta-actin. Since RPKM estimates of differentially expressed genes (DEGs) show bias towards overestimation of read count and transcript length, a different method was used to detect significant variation between O. furnacalis male and female antennal transcript levels. The algorithm DESeq assumes most transcripts do not represent DEGs, and implements a scaling factor which is calculated as the median ratio of read counts for each gene as a ratio of the geometric mean across all replicated libraries. The DESeq Bioconductor package v.1.6.0 for the R statistical package [55][56][57] was used here to construct the MA-plot-based method with random sampling model to identify DEGs. Differences in transcript abundances between male-and female-specific libraries were plotted on a log 2 (foldchange) scale and significance thresholds set at > 5-fold [50].
The level of transcripts OfurOR1 to 20 were estimated in male and female antennae, maxillary palpus and legs using real-time RT-qPCR, and the subset inclusive of the pheromone receptor subfamily used to validate RNA-seq estimates of gene expression. Total RNA for each tissue was prepared as described above, and then treated with DNase I (Invitrogen, New York, USA) to remove trace amounts of genomic DNA prior to cDNA synthesis. First strand cDNAs were synthesized by use of the AMV Reverse Transcription Kit (Promega, Wisconsin, USA). Primers for the O. furnacalis β-actin gene (accession number: GU301782) and 20 putative annotated OR genes were designed using Primer Express (ABI, USA), and respectively used to amplify target genes and a reference gene in separate reactions (S1 Table). Quantitative RT-qPCR was performed using an ABI Stepone Plus instrument with 20.0 μl reaction mixtures containing 10.0 μl SYBR Green qPCR Master Mix, 1.0 μl each of primers (10μM), 1.0 μl first strand antennal, maxillary palp or leg cDNA template and 7.0 μl ddH 2 O. Reactions were setup in triplicate for each template across all primer pairs (technical replicates), and repeated for 3 independent samples (biological replicates). The PCR program was as follows: 95°C for 2 m, followed by 40 cycles of 95°C 10 s and 60°C for 40 s, and then melt curve analysis was performed to test locus-specificity of reaction products. The data were analyzed by the comparative 2 -ΔΔCT method [58], with transcript levels normalized by comparison C T estimates from the beta-actin reaction. Transcript expression levels between male and female were compared as fold-change using the male antennae levels arbitrarily set at one, and significance of any relative difference between male and female expression (two conditions) was determined using paired T-tests as described by Pabiner et al., [59]. The correlation analysis was executed by using cor function of the R Statistical Package, and regression analysis was executed with the lm function [60].

Phylogenetic analysis of olfactory-related proteins
The O. furnacalis contigs that received functional annotations as putative OR-, OBP-, and CSP-like genes were retrieved from our reference antennal transcriptome assembly, and derived O. furnacalis amino acid sequence and putative translations for transcripts were made using Getorf from the EMBOSS package [49]. OR protein sequences previously identified in lepidopteran insects were downloaded from GenBank (71 from Bombyx mori and 1 from Conogethes punctiferalis, respectively). Nomenclature for all B. mori OR orthologs were retained from previously published analyses; BmorOR1-48 [40] BmorOR19j, 22j, 23j [61]; BmorOR49-68 [62]). A multiple amino acid sequence alignment was generated from downloaded and O. furnacalis ORs using the MUSCLE algorithm implemented using default parameters of MEGA 5.2.2 [63]. Alignments were similarly constructed for OBP using 10, 5, 19, 17, 11 sequences respectively from GenBank accessions for Chilo suppressalis, Manduca sexta, Helicoverpa armigera, Spodoptera exigua and B. mori. The nomenclature among B. mori OBP orthologs was used as described by Gong et al. [37]. The find best model of sequence evolution option of MEGA 5.2.2 [40] was used to evaluate both aligned OR and CSP sequences.
Phylogenetic reconstructions for OR and CSP orthologs were performed independently using neighbor-Joining (NJ) methods to evaluate the amino acid sequence evolution using the Jones-Taylor-Thorton (JTT) model, with node support generated from 1,000 bootstrap pseudoreplcations of the data. OfurOR13, 40,45,46,47,50,53,54, and 55 were omitted due to short derived protein sequence. Among site rate variation was accounted for with a gamma distribution of 4.994 and 3.599 respectively for the OR and CSP trees. Lepidopteran CSP protein family phylogenies used the neighbor-joining (NJ) method with Poisson correction of genetic distances. Trees were not rooted. All phylogenetic analyses used MEGA 5.2.2 [63].

Antennal transcriptome assembly and functional gene annotation
A combined assembly of these two datasets resulted in 37,687 contigs of 300 bp, which had with a mean length of 818 bp and the N50 length of 1,022 bp (3.02 million bp). A total of 37,687 clusters were obtained from CD-HIT (S1 Text), and represented the longest ORFs for each transcript. Among the 37,687 assembled O. furnacalis transcripts, a total of 15,544 showed "hits" following BLASTx homology search to the NCBI non-redundant (nr) protein database (E-value cut-off 10 −5 ; S2 Table). These BLASTx search results identified a total of 89 transcript sequences with putative homology to olfactory-related genes; 56 ORs ( Table 1) and 23 OBPs (Table A in S3 Table), and 10 CSPs (Table B in S3 Table). Of these transcripts, 14 showed > 95% sequence identity to previously identified O. furnacalis gene products already represented in GenBank: 5 PBPs (Accession number: GU828024 to GU828028; 1 GOBP2 (DQ673101) and 8 ORs (AB467327, JX910526, JN169134, JN169136, JN169138, JN169142, JX910532, JX910533; [16,32]; remaining search results not shown). The additional 75 putative olfactory-related transcripts found in the current study were not previously described in O. furnacalis or other species of Ostrinia. Of the 56 OfurOR transcripts (File A in S2 Text), the putative complete CDS was obtained for 13 (File B in S2 Text), wherein C-terminal CDS was obtained with greater prevalence. GO annotations from the longest sequence in each CD-HIT cluster ("UniGenes") were used to obtain function gene annotations. These annotations showed a high percentage of transcripts in GO Level 1 Cellular Component, GO Level 2; cell part (86.78%), cell (86.78%), and organelle (80.17%). Additionally, GO Level 1 Molecular Function showed the highest proportion of annotations in GO Level 2 binding (64.86%), and GO Level 1 Biological Process showed the greatest number of annotations in cellular process (76.03%) and metabolic process (76.03%). Contigs that received putative OR gene family member annotation were investigated in greater detail, where CDD and pfam database search results indicated that all 56 putative OfurORs encode 7tm_6 motifs which comprise 7 transmembrane domains that are typical of receptor proteins (pfam02949; superfamily cl20237; remaining data not shown).
Transcript levels were estimated between male and female O. furnacalis by RT-qPCR for OfurOR1 to OfurOR20 using antennae-, maxillary palps-and leg tissue-derived cDNA as template. These results showed that transcripts for OfurOR2, 3,4,6,8,11,13,14, and 18 were only Constituent contigs are provided along with the number of derived amino acid residues encoded in the full or partial CDS. The Top database hit for each OfurOR gene corresponding the Bombyx mori ortholog are listed along with percent protein identity (% ID), and do not represent ortholog predictions as predicted in Fig 4). detected in antennal tissues, whereas transcript for all the remaining OfurOR genes were also present in maxillary palpus and leg tissues. Most OfurOR genes showed higher estimated expression levels in antennal tissues, with the exception of OfurOR16, 20 (Fig 3). Transcripts that showed highly biased expression in male antennae included sex pheromone receptor subfamily members, OfurOR1, and 3 to 8, as well as OfurOR2, 9 to 15, 19, and 20 (Fig 3). In contrast, OfurOR17 and 18 transcript levels were higher in female compared to male antennae ( Table 1).
Comparisons of relative level of each transcript showed statistically significant differences between male and female antennal tissue (T-test P-values 0.0023), and included an estimated 3-fold greater level of OfurOR2 in male antennae (P-value = 0.0001; remaining results not shown). Comparison of the relative difference in C T values estimated from RT-qPCR experiments to log 2 (fold-change) estimates from RNA-seq data demonstrated that the two methods are not in 100% agreement across 20 OfurOR genes that were tested (Table 1). Although inconsistencies were observed, the male-specific sex pheromone receptor class (OfurOR1, 3 to 8) and female biased class (OfurOR17 and 18) were highly analogous between the two methods. Correlation between M:F ratios between RNA-seq and RT-qPCR showed an r = 0.57 and Pvalue = 0.009 (Fig 4).

Phylogenetic analyses of olfactory-related proteins
A 102 amino acid long consensus alignment was generated for 71 ORs from B. mori [40,61,62], Orco from C. punctiferalis and 56 putative OfurORs identified in this study. This alignment showed a (Ser/Ala)-Tyr-(Ser/Thr) C-terminal motif among OfurOR10, 13, 16, 26, 27, 34 and 37, which is believed to be conserved among B. mori ORs with female-biased expression [40]. Phylogenetic reconstruction based on this alignment indicated that OfurOR2/Orco clustered with known odorant receptor co-receptors from B. mori and C. punctiferalis. OfurOR1-8 proteins clustered into a well-supported monophyletic clade along with sex pheromone receptors from B. mori. Strong node support could not be obtained in all cases for one-to-one relationships between O. furnacalis and B. mori OR orthologs (Fig 5), but clustering of orthologs was observed in several instances (e.g. between BmorOR16 and OfurOR16, and BmorOR44 and OfurOR44) which might indicate that these ORs are derived from a common ancestral gene. Results also indicated that lineage-specific amplification of ORs may have taken place within O. furnacalis, where, for example OfurOR14, 20 and 23 all appear to share BmorOR23 as a most-recent common ancestral gene. A total of 23 putative OBPs were identified in the combined O. furnacalis antennal transcriptome, including genes encoding 5 PBP-and 2 GOBP-like proteins ( Table A in S3 Table). Phylogenetic analysis of O. furnacalis OBPs indicated that derived O. furnacalis OBPs clustered into ABPI, ABPII, CRLBP, Minus-C, Plus-C, and PBP/GOBP subfamilies as defined previously for homologous gene products from B. mori [37] (Fig 6). Specifically, OfurPBPs and Ofur-GOBPs clustered with homologous genes from B. mori, but one-to-one ortholog relationships were not defined for O. furnacalis within the pheromone binding protein (PBP) clade. A single member of the Minus-C subfamily, and two CRLBP transcript homologs were identified from the combined O. furnacalis antennal transcriptome, which is fewer than the 8 and 5 orthologs annotated from the B. mori genome assembly. Other O. furnacalis OBP classes similarly lacked putative B. mori orthologs, which might be due to lack of expression in antennal tissues or under treatment conditions used in this study.
A total of 10 CSPs were identified in the O. furnacalis transcriptome ( Table B in S3 Table). Phylogenetic reconstruction predicted that nine of these 10 putative OfurCSPs clustered with orthologs from other lepidopteran species with high node support (bootstrap values 87) (S1 Fig). In contrast, OfurCSP9 failed to cluster with other insect CSPs.

Discussion
Tissue-wide analysis of differential gene expression using RNA-seq data provides a seemingly un-precedent opportunity to define tissue-, stage-, and treatment-dependent variation. Comparison of sexual variation in gene expression in antennal tissues conducted in this experiment identified 109 significantly DEGs (Fig 1), which included the up-regulated pheromone receptor subfamily of ORs in male antennae. Several novel significantly DEGs were also identified such as those encoding putative senescence and seminal fluid proteins, an ABC transporter, and egg yolk proteins (S4 Table). This set of DEGs likely represent valuable candidate genes for understanding cellular processes that vary between the sexes, but are beyond the scope of the present research and admittedly suffer from a disconnect between statistical and biological significance. Specifically, it is likely not advisable to assume that all fold-changes between transcripts have biological consequences or are relevant to differential function of the olfactory system in male compared to female O. furnacalis without further testing. Regardless, steps toward understanding the molecular function of olfactory system components in species of Lepidoptera is important for deciphering how individuals perceive the environment, and the role of chemosensory reception in adult mate attraction and reproduction, oviposition site selection and host plant preference [1], and was focused on in this study.
A full genome sequence as well as gene expression data has been obtained for the first model species for Lepidoptera, B. mori [64,65]. The OR gene family in B. mori contains an estimated 71 paralogs [40,61,62] that have a highly variable primary sequence outside of a semiconserved C-terminal domain [40] and a 7 transmembrane domains (pfam02949; [34]). Sequence discovery and functional analyses have been critical in determining OR function in B. mori, where these data have facilitated determination that BmorOR1 stimulation is sufficient to elicit male sexual response [66], that binding of BmorOR56 by cis-jasmone is capable of mediating moth attraction to mulberry host plants leaves [62], and that female-specific expressed BmorOR19, 30, and 45-50 proteins can be stimulated by plant volatiles [61]. Species from the genus Ostrinia, including O. furnacalis, have emerged as models for the study of male sex pheromone detection systems [67] for which the function of antennal-expressed ORs during male detection of female emitted pheromones have been partially elucidated [31,33]. Despite these advances in uncovering the molecular function of Ostrinia ORs in male perception of female sex pheromones, little is known regarding the extent (number and diversity), expression, or function of many OR gene family members in Ostrinia. Functional diversification among members of a gene family often coincides with changes in temporal and spatial gene regulation [68], and also involve sex-specific patterns of expression in Lepidoptera. Specifically, regulation of PBP3 is biased for high expression in male compared to female antennae [16], and analogous male-biased expression was also shown for O. nubilalis OR3, 4, 5 and 6 [31]. In conjunction with in vivo pheromone response data collected from the Xenipus oocyte system, gene expression patterns are suggestive of a functional role for OR3, 4, 5 and 6 in male pheromone response [31,33]. This male-biased transcription appears to be retained among the B. mori orthologs OR3, 4, 5 and 6 [69]. Comparative genomic analyses suggest that male-biased expression and female pheromone receptor function is retained in this lepidopteran OR subfamily despite species divergence that spans over 100 million years. Analogously, female-biased transcription of OR gene family members is predicted among transcripts in both B. mori [40,61] and O. furnacalis. In the present study, analysis of read depth among aligned RNA-seq reads suggested 16.7-, 50.0-and 3.2-fold greater transcript levels respectively for OfurOR18, 17, and 16 in O. furnacalis female compared to male antennae (Table 1), and patterns were validated by RT-qPCR (Fig 3). Significant correlation between RNA-seq and RT-qPCR estimates of transcript levels were observed, although this should be interpreted with caution given the variance between paired estimates (Fig 4). Also, OfurOR54 and 26 comparatively showed 25-and 11-fold up-regulation in female O. furnacalis antennae that was significant at our highest stringent cutoff (Log 2 (fold-change) > 5.0). Phylogenetic analysis of B. mori and O. furnacalis OR orthologs indicate potential orthology between female biased BmorOR19 and female biased OfurOR21 and 54, as well as the clade comprised of female biased BmorOR45, 46 and 47 with OfurOR18 ( Fig 5). In contrast, the female biased transcripts OfurOR17 and 26 were not predicted to show any close orthologous relationship with female biased transcripts from B. mori.
Given that several B. mori female-biased ORs are capable of binding host plant volatiles [61,62], it is conceivable that O. furnacalis orthologs may have retained similar functions, but further studies are required to investigate any potential evolutionary conservation of function. Instances in which clear orthology was not predicted between female biased transcripts from B. mori and O. furnacalis, it may be hypothesized that similar expression pattern could have evolved within a different set of genes (e.g. identity by state). Since Ostrinia females tend to oviposit on a greater diversity of host plants compared to the more mulberry-specific attraction exhibited by B. mori females it is conceivable that female O. furnacalis adaptations may have selected additional ORs (potentially OfurOR17 and 26) that allow response to volatiles emitted from a greater range host plants. Thus, selection for ORs that respond to different host plant volatiles may have resulted functional diversification of unrelated genes in Ostrinia. Alternatively, divergent selection upon the same ancestral genes for binding to different host plant volatiles may resulted in a high degree of change in the protein sequence such that convergence with unrelated ORs may obscure the phylogenetic and orthologous relationships. The response of these Ostrinia ORs to different host plant volatiles remains unknown until appropriate function assays are performed, such that putative function in female host plant recognition prior to ovipostion cannot yet be established.
Courtship involves the emission of low-intensity utrasonic waves from specialized wing scales of male Ostrinia which causes conspecific females to become motionless [70] and results in increased frequency of male mating success [70,71]. Hair pencils are located on the 8 th sternite of O. nubilalis males, which when presented during courtship with cognate females, were shown to enhance mating success [72] and later discovered to contain cells that produce a blend of hexadecnyl acetates (male pheromone) [73]. These male pheromones are likely detected by potential female mates prior to copulation when females are observed to rapidly move antennae, suggesting that pheromone perception by females may be important for mate selection [74,75]. The mechanisms of female chemoreception of male pheromones remains unknown, but potentially involves neuronal stimulation by female-specific ORs in a system analogous to that which has been elucidated in male antennae. Functional characterization OfurOR stimulation in response to male pheromones using those ORs validated as up-regulated in female antennal might provide insight into female sexual acceptance, and may likely be the focus of future studies.
In addition to the role of ORs in eliciting neuronal signals in response to specific volatiles, CSPs and OBPs also may form important components of the chemosensory system way of shuttling hydrophobic volatiles from the peripheral environment to the ORs. PBPs are a type of OBP which are proposed to bind and chaperone sex pheromones across the sensillar lymph [34] and may, although contentious and still unresolved [34], have co-evolved with the OR gene family to provide additional selectivity in the chemosensory system [16]. Our RNA-seq data agree with prior RT-qPCR results that showed O. furnacalis PBP3 is significantly up-regulated in male antennae (Log 2 (fold-change) = 5.31). Since lepidopteran PBPs can bind ligands that are structurally similar female pheromones females and are expressed in non-pheromone sensitive tissues [76,77], PBPs and OBPs might have a role in chemosensory reception that is independent of sexual response and may have roles in general environmental perceptions.
The Asian corn borer, is widely distributed in countries from China to Australia, including Japan, Korea, and the Philippines, and is highly destructive to cultivated corn plants due to larval feeding damage to leaf, stalk and seed tissue that in-turn causes reduced crop yields. Larval O. furnacalis are also found on alternative host plants including bell pepper, cotton, hops, millet, pearl millet, foxtail millet, sugarcane, sorghum, and ginger as well as many weedy native plant species [78]. The attraction of female O. furnacalis to host plants for oviposition in the landscape may be important for understanding host range and for understanding chemoreception in lepidopteran insects. The evolution of female chemoreception that elicits an attraction and oviposition on host plants that are suitable for development of larval progeny is complex (see review [79]), selection of females that oviposit on plant that best support larval development may likely have shaped this female chemosensory system [80,81]. The molecular mechanisms involved in host plant attraction have been shown to potentially involved specific binding by CSPs in hemipteran insects [82] and OBPs in Drosophila [8], but the expression of OBPs and CSPs in non-olfactory sensitive tissues shown in this and prior studies might suggest these protein have alternate functions [83]. The capacity of putative O. furnacalis OBP and CSPs to bind various ligands have yet to be performed, but defining the structure and expression of CSP and OBP, as well and ORs, in this study represents a significant initial step that will allow future investigation of biological function.