Identification of chemosensory genes from the antennal transcriptome of Indian meal moth Plodia interpunctella

Olfaction plays an indispensable role in mediating insect behavior, such as locating host plants, mating partners, and avoidance of toxins and predators. Olfactory-related proteins are required for olfactory perception of insects. However, very few olfactory-related genes have been reported in Plodia interpunctella up to now. In the present study, we sequenced the antennae transcriptome of P. interpunctella using the next-generation sequencing technology, and identified 117 candidate olfactory-related genes, including 29 odorant-binding proteins (OBPs), 15 chemosensory proteins (CSPs), three sensory neuron membrane proteins (SNMPs), 47 odorant receptors (ORs), 14 ionotropic receptors (IRs) and nine gustatory receptors (GRs). Further analysis of qRT-PCR revealed that nine OBPs, three CSPs, two SNMPs, nine ORs and two GRs were specifically expressed in the male antennae, whereas eight OBPs, six CSPs, one SNMP, 16 ORs, two GRs and seven IRs significantly expressed in the female antennae. Taken together, our results provided useful information for further functional studies on insect genes related to recognition of pheromone and odorant, which might be meaningful targets for pest management.


Introduction
Indian meal moth, Plodia interpunctella (Hübener) (Lepidoptera: Pyraloidea, Pyralidae), is a notorious stored-product pest worldwide [1]. The larvae infest a variety of processed foods, including fruits, nuts, cereals, powdered milk, chocolate, birdseed, and pet food [2], causing extensive damage by impairing dry weight, germination, nutritional value, and quality grade. It is difficult to control P. interpunctella by conventional insecticides, because it often inhabits our kitchen, closet and warehouse, and its larvae are mixed with our processed foods. Accordingly, several novel strategies have been developed to monitor and control P. interpunctella. Among these novel methods, sex pheromone is widely acceptable due to its safety and efficiency. Meanwhile, host volatiles have been thought to affect the oviposition behavior of P. PLOS

Insects material and RNA extraction
Plodia interpunctella was the laboratorial population which was reared for more than 20 generations in our laboratory. The larvae were reared on crushed grains of wheat under constant conditions (28±1˚C, 60±5% RH and 14:10 L:D photoperiod). Mature larvae were sorted by sex according to the black spot in the middle of male back. Antennae were excised from 3-day-old unmated moths, immediately frozen in liquid nitrogen and ground with a pestle. Total RNA was extracted from 100 antennae for each sex. The evaluation of RNA purity, RNA concentration and RNA quality were conducted following our previous method [13].
cDNA library preparation for transcriptome sequencing cDNA library were constructed following previous method [17]. Briefly, 3 μg RNA per sample was used as input material for the RNA sample preparation. Sequencing libraries were generated using NEBNext1Ultra™ RNA Library Prep Kit for Illumina1 (NEB, USA) following manufacturer's instructions. Newly isolated mRNA was further purified using with Oligo (dT) magnetic beads and sheared into 200-700 nucleotides sections using fragmentation buffer.
The fragmented mRNA was used as templates for first-strand cDNA synthesis using random hexamer primers. Subsequently, second-strand cDNA was synthesized using DNA polymerase I and RNaseH. All remaining overhangs were passivated via polymerase. After adenylation of 3 0 ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated for hybridization. In order to select cDNA fragments of preferentially 150~200 bp, the library fragments were purified using an AMPure XP system. Then 3 μL USER Enzyme (NEB, USA) was incubated with size-selected, adaptor-ligated cDNA at 37˚C for 15 min followed by incubation at 95˚C for 5 min before PCR reaction. Subsequently, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Amplicons were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The cDNA library of P. interpunctella was sequenced on Illumina Hiseq™ 2500 using paired-end technology in a single run by Beijing Biomake Company (Beijing, China).

Clustering and sequencing
Following a previous report [17], clustering and sequencing were performed on a cBot Cluster Generation System and an Illumina Hiseq 2500 platform, respectively.

Sequence analysis and assembly
Raw reads of fastq format were firstly processed through in-house perl scripts. In this step, clean reads were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads. At the same time, Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. Cleaned reads shorter than 60 bases were removed because the short reads might represent sequencing artifacts [26]. The qualified reads were assembled into unigenes using short reads assembling program-Trinity [10]. The obtained contigs were annotated against the NCBI non-redundant protein (NR) database using BLASTn (E-value<10 −5 ) and BLASTx (E-value<10 −5 ) programs [11]. To annotate the assembled sequences with Gene Ontology (GO) terms, the Swiss-Prot BLAST results were imported into BLAST2GO, a software package that retrieves GO terms, allowing determination and comparison of gene functions [27]. The unigene sequences were also aligned to the Clusters of Orthologous Groups of proteins (COG) database to predict and classify the unigene sequences [28]. Pathway annotations for unigenes were determined using Kyoto Encyclopedia of Genes and Genomes (KEGG) ontology [29]. Finally, the best matches were used to identify coding regions and determine the sequence direction [30].

Olfactory gene identification and phylogenetic analysis
The annotations of OBP, CSP, SNMP, OR, IR and GR genes in P. interpunctella were verified by BLASTx and BLASTn programs NCBI. The complete coding region was predicted using the open reading frame (ORF) finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) based on the results given by BLASTx. After completing the alignments of the candidate chemosensory genes using ClustalX (2.1), phylogenetic reconstruction for the analysis of OBPs, CSPs, ORs, IRs and GRs was performed by MEGA5.0 software using the neighbor-joining method with 1000 Bootstrap iterations [31]. In addition, the evolutionary distances were assumed by using the Poisson correction method [11].

Analysis of differentially expressed genes and qRT-PCR verification
To compare the differential expression of chemosensory genes between the male and female antennal transcriptomes of P. interpunctella, the read number of each olfactory-related gene was converted to FPKM (fragments per kilobase of exon model per million mapped reads) [32].
qRT-PCR was performed to quantify the expression levels of olfactory-related genes in male and female antennae. Total RNA was extracted from 100 antennae as above description. cDNA from antennae of both sexes was synthesized using the SMART TM PCR cDNA synthesis kit (Clontech, Mountain View, CA, USA). The β-actin gene (SRP05157) was used as an internal control in each sample, and it was selected as a housekeeping gene in our qRT-PCR test. Real-time PCR was performed on an ABI 7500 using SYBR green dye binding to doublestranded DNA at the end of each elongation cycle. Primer sequences were designed using the Primer Premier 5.0 program (S1 Table). Real-time PCR was conducted with our previous method [13]. Briefly, 10.0 μL of 2×SYBR Green PCR Master Mix, 0.4 μL of primer, 2.0 μL of sample cDNA (100 ng μL -1 ) and 7.2 μL of sterilized ultrapure water were mixed to form a 20 μL reaction system.After an initial denaturation step at 95˚C for 3 min, amplifications were carried out with 40 cycles at a melting temperature of 95˚C for 10 s and an annealing temperature of 60˚C for 30 s. To check reproducibility, qRT-PCR test for each sample was performed with three technical replicates and three biological replicates.

qRT-PCR analysis
Relative quantification was determined using the comparative 2 -ΔΔCt method [33]. All data were normalized to endogenous β-actin levels from the same individual samples. The relative fold change was assessed by comparing the expression level in male moths to that in females [34]. The results were presented as the means of the fold change in three biological duplicates. The comparative analyses of chemosensory genes between sexes were determined by one-way analysis of variance (ANOVA) using SPSS 19.0, with p-value of 0.05 considered significant.

Results
Sequence analysis and assembly cDNA library of Plodia interpunctella was constructed using the TRINITY de novo assembly program, and short-read sequences were assembled into 150,633 transcripts with a mean length of 1,491 bp and an N50 of 3,567 bp. A total of 20,261 scaffolds (13.45%) were longer than 1,000 bp, and 36,148 scaffolds (24.00%) were longer than 2,000 bp. The scaffolds were subjected to cluster and assembly analyses. Subsequently, 87,300 unigenes were obtained with a mean length of 699 bp and an N50 of 1,282 bp (Fig 1, Table 1). The length distribution of unigenes revealed that 26,054 unigenes (29.84%) were longer than 500 bp and 12,485 unigenes (14.30%) were longer than 1,000 bp ( Table 1). The raw reads of P. interpunctella transcriptome have been deposited into the NCBI SRA database (accession number: SRR6002827 and SRR6002828), and the Transcriptome Shotgun Assembly (TSA) project has been deposited at DDBJ/ENA/GenBank under the accession GFWQ00000000. The version described in this paper is the first version, GFWQ01000000. The detailed TSA sequences could be obtained from Genbank (https://www.ncbi.nlm.nih.gov/Traces/wgs/?val=GFWQ01&display= contigs&page=1).

Sequence annotation
The unigene annotation showed that 27,920 unigenes (31.98%) significantly matched in the NR database and 15,815 unigenes (18.12%) had significant matches in the Swiss-Prot database. A total of 31,921 unigenes (36.56%) were successfully annotated in the NR, Swiss-Prot, KEGG, GO and COG databases ( Table 2), whereas 55,379 unigenes (63.44%) were unmapped in those databases.
For GO analysis, 15,893 unigenes (18.21%) could be assigned to three GO terms as follows: cellular components, molecular functions and biological process (Fig 3). The "cellular components" and "molecular functions" were most represented by 18.79% and 21.04% transcripts, respectively. In the "cellular components" ontology, the terms were mainly distributed in cell  (20.71%) and cell part (20.71%). In the "molecular functions" ontology, the terms of binding function and catalytic activity were the most represented (39.91% and 39.80%, respectively) (Fig 3).
To predict and classify the functional genes, all unigenes were searched against the COG database. A total of 10,106 unigenes could be assigned to 25 specific categories according to the COG classification results. "General function prediction" (2,494, 24.68%) was the largest group, and the categories of "cell motility" (20, 0.20%) and "nuclear structure" (11, 0.11%) were the smallest groups (Fig 4). In addition, 290 pathways were predicted in the KEGG database, representing 15,016 unigenes.

Identification of olfactory-related genes
In the present study, we identified 117 olfactory-related genes from antennal transcriptome of P. interpunctella, including 29 OBPs, 15 CSPs, three SNMPs, 47 ORs, nine GRs and 14 IRs. All genes were named according to a four-letter code (first letter of the genus name followed by the first three letters of the species name) + OR + number according to the ORF lengths. Analysis of differential expression of unigenes indicated that 1,031 genes showed differences between the antennal transcriptomes of male and female P. interpunctella, including 93 upregulated and 938 down-regulated genes using female result as the reference standard.

Candidate ORs
We identified 47 OR genes in the antennal transcriptomes of P. interpunctella, in which 36 Pin-tORs had intact ORFs with lengths ranging from 219 bp to 1,422 bp with four to seven transmembrane domains (Table 5).
In the neighbor-joining tree of ORs (Fig 11), the PintOR1 was clustered into the ORco family and four PintORs (PintOR5, PintOR7, PintOR22 and PintOR30) were clustered into the pheromone receptor (PR) family. Two groups of ORs (PintOR14 and PintOR35, PintOR29 and Pin-tOR26) were clustered into the same branch with bootstrapping values of 98 and 87, respectively. All of the other PintORs were distributed on various branches throughout the phylogenetic tree.

Candidate GRs
In the present study, we identified nine candidate PintGR encoding transcripts from antennal transcriptome of P. interpunctella. Five PintGR genes had intact ORFs with lengths ranging from 198 bp to 1,461 bp. The BLASTx results indicated that seven identified PintGRs shared relatively higher amino acid identities (>50%) with Lepidoptera GRs in NCBI (Table 6).
In the neighbor-joining tree of GRs (Fig 13), PintGRs were present on various branches throughout the cladogram. PintGR1 and PintGR3 were clustered into the same branch, with a bootstrapping value of 65.

Candidate IRs
In the present study, we identified 14 candidate PintIR genes encoding transcripts from antennal transcriptome of P. interpunctella (Table 6). Nine PintIRs had intact ORFs with lengths ranging from 384 bp to 2,706 bp. In the neighbor-joining tree of IRs (Fig 15), PintIR1 and Pin-tIR2 were phylogenetically clustered into the highly conserved IR8a and IR21a sub-families, respectively. The FPKM analysis revealed that all PintIRs showed a low expression level (FPKM value ranged from 0.36 to 113.52). The qRT-PCR results indicated that PintIR1, Pin-tIR3-5, PintIR10, and PintIR13-14 were highly expressed in the female antennae (1.2 to 5.3 times compared with males) (Fig 14). Identification of chemosensory genes of Plodia interpunctella

Discussion
In recent years, RNA-Seq transcriptome sequencing technology has been widely used due to the development of high-throughput sequencing technology, resulting in great progress in non-model organisms [11,[37][38][39]. In the present study, we used NGS technology to analyze  Table. https://doi.org/10.1371/journal.pone.0189889.g011 Identification of chemosensory genes of Plodia interpunctella the antennal transcriptome of P. interpunctella. Sequence analysis and assembly results demonstrated that Illumina sequencing technology could effectively and rapidly captured a large portion of the transcriptome, providing molecular foundations for rapid characterization of functional genes and better reference of target genes [40]. The unigene annotation showed that 55,379 unigenes (63.44%) were unmapped in those databases, which could be attributed to the short sequence reads generated by the sequencing technology. It also suggested that the unmapped sequences could represent unannotated or new genes. In fact, fewer than 5% of unmapped unigenes are likely to represent new genes. Generally, the 5' ends of sequences show less conservation than the body. Therefore, partial transcripts (unigenes representing the 5' CDS, but not the body) may not be found matches in the databases. For GO analysis, the antennal unigenes were annotated into different functional groups [16], which were similar to those in the antennal transcriptomes of Conogethes punctiferalis [13], Spodoptera littoralis [41] and Helicoverpa armigera [11]. Therefore, we inferred that the success rates of functional annotation of genes highly depended on the sequence length of the splicing unigene: the shorter the length of the sequence, the less possibility of the annotation. Others reasons might also result in partial information failure, such as the Identification of chemosensory genes of Plodia interpunctella incompleteness of P. interpunctella gene transcription group information, and/or the insufficiency of the sequence of partial RNA-Seq sequencing data in public database.
Olfactory-related genes might be used as potential targets for management programs of P. interpunctella. As the first step of odor detection [6], OBPs have attracted wide interests of researchers [13,17,42]. In the present study, we identified 29 PintOBP genes from antennal transcriptome of P. interpunctella. The number of identified PintOBPs was equivalent to that from H. armigera (26) [11], Dendrolimus kikuchii (27) [17] and Agrotis ipsilon (33) [21], and it was significantly greater than that from Cnaphalocrocis medinalis (12) [12], C. punctiferalis (14) [13], Manduca sexta (18) [43] and S. exigua (11) [44]. The small number of OBPs in above species could be attributed to that the actual number of OBPs was less in P. interpunctella, or there should be more OBPs that were not caught by the sequencing. Therefore, we speculated that the transcriptomic sequencing might not be strong enough to detect all the OBPs, especially for some OBPs with low expression levels in the antennae [45].
The OBP trees from five Lepidopteran species indicated that after a long history evolution, the Lepidopteran OBPs differentiated into several branches (Fig 5), which was consistent with previous reports [46]. In the evolutionary tree for GOBPs and PBPs, these two sub-families were clustered respectively, indicating that these genes might have the same ancestor gene and differentiate along sex isolation and speciation. The qRT-PCR results indicated that nine PintOBP genes (PintOBP4, PintOBP6, PintOBP9, PintOBP13, PintOBP17 PintOBP20, Pin-tOBP22 and PintPBP2-3) were significantly expressed in the male antennae, suggesting that these OBPs played essential roles in the detection of sex pheromones. Eight PintOBPs (Pin-tOBP5, PintOBP7, PintOBP12, PintOBP15-16, PintOBP18, PintPBP1 and PintGOBP1) were significantly expressed in the female antennae, revealing that these OBPs played important roles in the detection of general odorants, such as host plant volatiles [21].
CSPs represent a newly-discovered class of soluble carrier proteins with similar functions to OBPs in insect chemoreception [47]. CSPs have been found in insect chemosensory tissues and non-chemosensory organs, such as antennae [11], legs [48], labial palps [49], tarsi [50], brain [51], proboscis [52], pheromone gland [53][54] and wings [55]. We identified 15 putative CSP encoding transcripts, and found that six PintCSP genes were significantly expressed in the female antennae. These PintCSPs might play important roles in the detection of general odorants, such as host plant volatiles. Identification of chemosensory genes of Plodia interpunctella OR proteins are key players in insect olfaction [56]. We identified 47 PintOR genes in antennal transcriptome of P. interpunctella. The number of PintORs identified in this study was less than that identified from the antennal transcriptomes of Bombyx mori (72) [57], C. punctiferalis (62) [58] and Ostrinia furnacalis (56) [59]. However, the difference in identified OR gene numbers might be caused by sequencing methods and depth, or sample preparation. In the neighbor-joining tree of ORs, four PintORs (PintOR5, PintOR7, PintOR22 and Pin-tOR30) were clustered into the PR family, indicating that parts or all of them contributed to sex pheromone detection. The qRT-PCR results indicated that PintOR5 and PintOR22 were highly expressed in the male antennae, suggesting they are highly related to sex pheromone. PintOR7 and PintOR30 specifically expressed in the female antennae. The expression profiles of these sequences showed that not all of them were male-specific [60]. Recent studies also showed that some PR genes are expressed in both sexes [54]. The OR tree showed that the Pin-tORco (PintOR1) was highly conserved.  Table. In recent years, 12 HarmIRs in H. armigera [11], 17 SlitIRs in S. littoralis [61] and 15 Cpo-mIRs in C. pomonella [22] have been identified. In this study, we identified 14 PintIRs, including highly conserved IR co-receptors PintIR1 and PintIR2 (IR8a and IR21a) from antennal transcriptome of P. interpunctella. Therefore, we speculated that IRs were relatively highly conserved sequences, implying that IRs had conservative features.
Several recent reports simultaneously tested the qRT-PCR expression of olfactory-related genes in various tissues of insect, including bodies, heads, legs or abdomens [13-15, 20-21, 54]. In present study, we only focused on the qRT-PCR analysis of P. interpunctella antennae. To the best of our knowledge, P. interpunctella moths do not eat anything, suggesting they have no food demands, so location of mate partners and oviposition sites should be the main function of olfactory. While most olfactory genes related to recognition of pheromone and host volatiles distribute in insect antennae, therefore, we only compared the expression between female and male antennae of P. interpunctella, to verify the olfactory-related genes.

Conclusion
In this study, we identified a few olfactory gene families in antennal transcriptome of P. interpunctella, including 29 PintOBPs, 15 PintCSPs, three PintSNMPs, 47 PintORs, nine PintGRs and 14 PintIRs. The identification of antennal olfactory-related proteins in P. interpunctella reinforced our knowledge on the molecular and cellular basis of insect chemoreception. More importantly, our data suggested that new methods could be developed to control this pest by interfering their olfactory perception.
Supporting information S1