Identification of Chemosensory Genes Based on the Transcriptomic Analysis of Six Different Chemosensory Organs in Spodoptera exigua

Insects have a complex chemosensory system that accurately perceives external chemicals and plays a pivotal role in many insect life activities. Thus, the study of the chemosensory mechanism has become an important research topic in entomology. Spodoptera exigua Hübner (Lepidoptera: Noctuidae) is a major agricultural polyphagous pest that causes significant agricultural economic losses worldwide. However, except for a few genes that have been discovered, its olfactory and gustatory mechanisms remain uncertain. In the present study, we acquired 144,479 unigenes of S. exigua by assembling 65.81 giga base reads from 6 chemosensory organs (female and male antennae, female and male proboscises, and female and male labial palps), and identified many differentially expressed genes in the gustatory and olfactory organs. Analysis of the transcriptome data obtained 159 putative chemosensory genes, including 24 odorant binding proteins (OBPs; 3 were new), 19 chemosensory proteins (4 were new), 64 odorant receptors (57 were new), 22 ionotropic receptors (16 were new), and 30 new gustatory receptors. Phylogenetic analyses of all genes and SexiGRs expression patterns using quantitative real-time polymerase chain reactions were investigated. Our results found that several of these genes had differential expression features in the olfactory organs compared to the gustatory organs that might play crucial roles in the chemosensory system of S. exigua, and could be utilized as targets for future functional studies to assist in the interpretation of the molecular mechanism of the system. They could also be used for developing novel behavioral disturbance agents to control the population of the moths in the future.

Insects have a complex chemosensory system that accurately perceives external chemicals and plays a pivotal role in many insect life activities. Thus, the study of the chemosensory mechanism has become an important research topic in entomology. Spodoptera exigua Hübner (Lepidoptera: Noctuidae) is a major agricultural polyphagous pest that causes significant agricultural economic losses worldwide. However, except for a few genes that have been discovered, its olfactory and gustatory mechanisms remain uncertain. In the present study, we acquired 144,479 unigenes of S. exigua by assembling 65.81 giga base reads from 6 chemosensory organs (female and male antennae, female and male proboscises, and female and male labial palps), and identified many differentially expressed genes in the gustatory and olfactory organs. Analysis of the transcriptome data obtained 159 putative chemosensory genes, including 24 odorant binding proteins (OBPs; 3 were new), 19 chemosensory proteins (4 were new), 64 odorant receptors (57 were new), 22 ionotropic receptors (16 were new), and 30 new gustatory receptors. Phylogenetic analyses of all genes and SexiGRs expression patterns using quantitative real-time polymerase chain reactions were investigated. Our results found that several of these genes had differential expression features in the olfactory organs compared to the gustatory organs that might play crucial roles in the chemosensory system of S. exigua, and could be utilized as targets for future functional studies to assist in the interpretation of the molecular mechanism of the system. They could also be used for developing novel behavioral disturbance agents to control the population of the moths in the future.
The acquisition, bioinformatics analysis, and expression pattern of putative chemosensory genes are the crucial steps to explore the exact roles of several key genes in the insect chemosensory process. The development of modern molecular biology techniques and experimental equipment, such as high-throughput sequencing, has created more efficient, inexpensive, and higher accuracy technologies than what has been traditionally utilized (McKenna et al., 1994;Picimbon and Gadenne, 2002;Xiu et al., 2008;Liu et al., 2012). These have been successfully applied in the identification of insect chemosensory genes, including many moth species, such as Spodoptera littoralis (Legeai et al., 2011), Sesamia inferens , Helicoverpa armigera (Liu et al., 2014b), Plutella xylostella , and Ectropis grisescens (Li et al., 2017).
The beet armyworm, Spodoptera exigua Hübner (Lepidoptera: Noctuidae), is a major agricultural polyphagous pest that causes significant economic losses to many crops worldwide (Xiu and Dong, 2007;Acín et al., 2010;Lai and Su, 2011). To date, only partial chemosensory genes of S. exigua have been identified, including several OBPs (Xiu and Dong, 2007;Zhu et al., 2013;Liu et al., 2015b), CSPs (Liu et al., 2015b) and a few chemosensory receptor genes (Liu et al., 2013(Liu et al., , 2014a(Liu et al., , 2015b. This is much lower than other moth species from which chemosensory genes have been obtained from transcriptomic data of chemosensory organs. These limited gene resources impede our interpretation of the chemosensory molecular mechanism of S. exigua. To obtain greater olfactory and gustatory gene resources, we utilized the six major olfactory and gustatory organs (female antennae: FA, male antennae: MA, female proboscises: FPr, male proboscises: MPr, female labial palps: FLP, and male labial palps: MLP) of S. exigua adults in the present study. We first built a genetic database of genes that were expressed in the six chemosensory organs of S. exigua using an Illumina HiSeq TM 4000 sequencing platform and completely identified 159 genes (110 genes were newly obtained) as being potentially involved in the chemosensory system. To postulate the functions of these identified genes, we performed phylogenetic analyses of all genes and investigated SexiGRs expression patterns using quantitative real-time polymerase chain reaction (qPCR). Our results showed that several of the genes had differential expression in olfactory organs compared to gustatory organs that might play different and crucial roles in the chemosensory system of S. exigua, and could be utilized as targets for future functional studies (using the heterologous expression system of Xenopus oocytes or Escherichia coli in vitro and with genetic modification by the CRISPR/Cas9 editing system in vivo) to assist in the interpretation of the molecular mechanism of the system.

Insects Rearing and Tissue Collection
S. exigua larvae were purchased from Keyun Biology Company in Henan province, China. As we previous studies (Zhang et al., 2017a), we used same rearing conditions and methods to rear the insect. For transcriptome sequencing, 200 female antennae (FA), 200 male antennae (MA), 300 female proboscises (FPr), 300 male proboscises (MPr), 300 female labial palps (FLP), 300 male labial palps (MLP), 30 female abdomen (FAb), and 30 male abdomen (MAb) were collected from 3-day-old unmated adults. For the tissue distribution analysis, 100 FA, 100 MA,200 FLP,200 MLP,200 FP,and 200 MP for each replicate experiment were collected under the same conditions. All these organs were immediately frozen in liquid nitrogen and stored at −80 • C until use.
cDNA Library Preparation, Clustering, and Sequencing Sample total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). cDNA library preparation and Illumina sequencing were carried out by Novogene Bioinformatics Technology Co., Ltd. (Beijing, China). The 1.5 µg total RNA per sample was used as input material for the RNA sample preparations, and sequencing libraries were generated using NEBNext R Ultra TM RNA Library Prep Kit for Illumina R (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-) (NEB, USA). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I (NEB, USA) and RNase H (NEB, USA). Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3 ′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 150∼200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 µL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq TM 4000 platform and paired-end reads were generated.

Differential Expression Analysis
Firstly, the read counts were adjusted by edgeR 3.0.8 program package through one scaling normalized factor for each sequenced library. Then, the differential expression analysis of two samples was performed using the DEGseq 1.12.0 R package (Wang et al., 2010). P-value was adjusted using q-value (Storey, 2003). q < 0.005 & |log2(foldchange)|>1 was set as the threshold for significantly differential expression.

RNA Isolation and cDNA Synthesis
Total RNA was extracted using the MiniBEST Universal RNA Extraction Kit (TaKaRa, Dalian, China), following the manufacturer's instructions, in which we used DNase I to digest sample DNase to avoid genomic DNA contamination. The RNA quality was assessed spectrophotometrically (Biofuture MD2000D, UK). Single-stranded cDNA templates were synthesized from 1 µg total RNA obtained from various tissue samples using the PrimeScript TM RT Master Mix (TaKaRa, Dalian, China) according to the manufacturers' instructions.

Quantitative Real-Time PCR (qPCR) Analysis
According to the minimum information for publication of qPCR experiments (Bustin et al., 2009) and our previous studies (Zhang et al., 2017a), we performed the qPCR assay of tissue distribution of SexiGRs in ABI 7300 (Applied Biosystems, Foster City, CA, USA) by using 2×SYBR Green PCR Master Mix (YIFEIXUE BIO TECH, Nanjing, China) as the manufacturer's instructions. Briefly, the reaction programs were 10 min at 95 • C, 40 cycles of 95 • C for 15 s and 60 • C for 1 min. The qPCR primers (Table S2) were designed using Beacon Designer 7.9 (PREMIER Biosoft International, CA, USA). Then, the relative expression levels of SexiGRs mRNA were calculated based on the Ct-values of target gene and two reference genes SexiGAPDH (glyceraldehyde-3-phosphate dehydrogenase) and SexiEF (elongation factor-1 alpha) by using the Q-Gene method in Microsoft Excel-based software of Visual Basic (Muller et al., 2002;Simon, 2003), the qPCR data are listed in Table S3. To ensure the reliability of the results, we carried out three biological replications for each sample and three technical replications for each biological replication.

Statistical Analysis
Data (mean ± SE) from various samples were subjected to oneway nested analysis of variance (ANOVA), followed by the least significant difference test (LSD) for comparison of means using SPSS Statistics 22.0 (SPSS Inc., Chicago, IL, USA).

Overview of Transcriptomes From the Six Organs
We used next-generation sequencing to sequence the six cDNA libraries constructed from the chemosensory organs (FA, MA, FPr, MPr, FLP, and MLP) of S. exigua adults based on the Illumina HiSeq TM 4000 platform and acquired 65.81 (from 10.60 to 11.90) giga base reads. After clustering and redundancy filtering, we finally obtained 144,479 unigenes and 266,645 transcripts with FIGURE 1 | Distribution of unigene size in the S. exigua transcriptome assembly.
FIGURE 2 | Percentage of homologous hits of the S. exigua unigenes to other insect species. The S. exigua transcripts were searched by Blastx against the Nr protein database with a cutoff E-value 10 −5 . Species that have more than 1% matching hits to the S. exigua transcripts are shown.
Frontiers in Physiology | www.frontiersin.org  a N50 length of 2,177 base pair (bp) and 1,552 bp, respectively ( Table 1). Statistics showed that 59.22% of the 144,479 unigenes were greater than 500 bp in length (Figure 1). The number of reads, unigenes, and transcripts were higher than most other insects based on transcriptome studies. In total, 60,373 unigenes were matched to entries in the National Center for Biotechnology Information (NCBI) non-redundant (NR) protein database (http://www.ncbi. nlm.nih.gov/protein) by a BLASTX homology search with a cut-off e-value of 10 −5 . The highest match percentage (37.40%) was identified with sequences of Bombyx mori followed by sequences of Danaus plexippus (15.60%), P. xylostella (13.20%), Homo sapiens (4.30%), and H. armigera (1.40%; Figure 2).    Based on methodology described in our previous studies Li et al., 2015), we applied Blast2GO to classify the functional groups of all unigenes. The results showed that only 29.29% (42,331) of the 144,479 unigenes could be annotated based on the sequence homology, with this proportion similar to that found in other insects (Gu et al., 2013;Zhang et al., 2013;He et al., 2017). One possible reason for this might be that a great amount of S. exigua unigenes belong to non-coding or homologous genes without a gene ontology (GO) term. In addition, the GO annotation of S. exigua unigenes displayed similar classification to the unigenes of chemosensory organs from other moth species (Grosse-Wilde et al., 2011;Zhang et al., 2013;Cao et al., 2014;Xia et al., 2015). For example, unigenes of S. exigua during biological processes were predicted to be mostly enriched in three sub-categories: cellular, metabolic, and singleorganism processes. There was also expected to be similarity in the cellular components (e.g., cell, cell part, and organelle) and molecular function categories (binding, catalytic, and transporter activity; Figure 3), indicating that some unigenes in these sub-categories might play important roles in the chemosensory behavior of moths.

Differentially Expressed Genes (DEGs)
To investigate the DEGs among different organs, we compared each organ pair-wise within each sex against all other organs (Figure 4). Gene expression dynamics can be reflected by upor down-regulation among the six different organs by pairwise comparisons. The results showed that there were a number of DEGs between different organs and different sexes, and the number of DEGs was highest in FPr vs. FLP (6,029 genes in total: 4,050 up-regulated genes and 1,979 downregulated genes), followed by MA vs. FPr (5,127 genes in total: 1,928 up-regulated genes and 3,199 down-regulated genes), and MPr vs. FLP (4,033 genes in total: 2,513 up-regulated genes and 1,520 down-regulated genes). This indicates that these DEGs, especially in the gustatory vs. olfactory organs, provide substantial genetic sources that are important for studying the differential mechanism of gustatory vs. olfactory  Table S1. This tree was constructed using PhyML based on alignment results of ClustalX. organs in S. exigua. Additionally, they provide some important target genes to analyse the functions of expressed sex-specific genes to reveal sex differences in chemosensory mechanisms in the future.

Identification of Putative Chemosensory Genes
Based on sequence similarity analyses and characteristics of insect chemosensory genes from previous studies (Xu et al., 2009;Croset et al., 2010;Zhou, 2010;Zhang et al., 2011;Ray et al., 2014), such as the conserved C-pattern of OBPs and CSPs, and the conserved transmembrane structure and motifs of chemosensory receptors (ORs, IRs, and GRs), we totally identified 159 putative genes from the transcriptomic data of S. exigua chemosensory organs that belonged to five insect chemosensory gene families. These included 24 OBPs, 19 CSPs, 64 ORs, 22 IRs, and 30 GRs (Tables 2, 3). The number of putative chemosensory genes of S. exigua identified in the present study was higher than that in other moth species where the same family genes had been identified by analysis of the transcriptome of specific organs. This included H. armigera (143 genes: 34 OBPs, 18 CSPs, 60 ORs, 21 IRs, and 10 GRs) (Liu et al., 2014b), H. assulta (147 genes: 29 OBPs, 17 CSPs, 64 ORs, 19 IRs, and 18 GRs) (Xu et al., 2015), and P. xyllostella (116 genes: 24 OBPs, 15 CSPs, 54 ORs, 16 IRs, and 7 GRs) . We found that the amount of transcriptomic data of these three different moth species was less than that of S. exigua in the present study, which suggests that the large amount of transcriptomic data could help us obtain more insect chemosensory genes.

OBPs
We obtained a complete set of 24 different unigenes encoding putative OBPs in S. exigua (Table 2), of which 3 were newly identified. Sequence analysis revealed that 23 sequences were predicted to have full-length open reading frames (ORFs) and encoded 118-239 amino acids, but only 3 of the 23 SexiOBPs did not have signal peptide sequences ( Table 2). The phylogenetic analysis showed that all 24 SexiOBPs were clustered in an OBP tree with Manduca sexta, B. mori, and Athetis lepigone  Table S1. This tree was constructed using PhyML based on alignment results of ClustalX. (Figure 5), including 5 SexiOBPs (SexiPBP1-3, SexiGOBP1-2) clustered into the PBP/GOBP subfamily. The results suggest that these SexiOBPs belonged to the insect OBP family and should have the corresponding functions of the insect OBP (Poivet et al., 2012;Jeong et al., 2013;Pelosi et al., 2014;Liu et al., 2015a). The two new SexiOBPs (SexiOBP-N1 and SexiOBP-N3) encoded protein with high identities (97 and 99%) to OBPs in Spodoptera litura, respectively, indicating that SexiOBP-N1 and SexiOBP-N3 might have conserved functions in the two closely related species, such as recognizing the same host plant volatiles (Li et al., 2013;Gu et al., 2015). Therefore, they can be considered as target genes to simultaneously prevent and control these two pests (S. exigua and S. litura) in the future.

CSPs
Nineteen putative genes encoding CSPs were acquired in S. exigua based on the analysis results from the transcriptomes of the six chemosensory organs, of which four were newly attained ( Table 2). Among the 19 SexiCSPs, 18 had full length ORFs with 4 conserved cysteines in the corresponding position and a predicted signal peptide at the N-terminus. The constructed insect CSP tree using protein sequences from S. exigua, M. sexta, B. mori, and A. lepigone (Figure 6) indicated that all 19 SexiCSPs were distributed along various branches and each clustered with at least 1 other moth ortholog. Thus, we inferred that these SexiCSPs should have a similar chemosensory function in insects, especially moths (Lartigue et al., 2002;Campanacci et al., 2003;Zhang et al., 2014). Similar to SexiOBPs, we also found three of the four new SexiCSPs (SexiCSP-N2, SexiCSP-N3, and SexiCSP-N4) encoded proteins with high identities (99 and 100%) to CSPs in S. litura. This showed that they were very similar, maybe even the same CSPs, and might play the same role as OBPs in the two moths. In future studies, we intend to use the combination of in vitro (Jin et al., 2014;Zhang et al., 2014) and in vivo (Zhu et al., 2016;Dong et al., 2017;Ye et al., 2017) methods to explore the exact function of these conserved OBPs and CSPs in the two closely related species. In addition, we plan to study the exact functions of all the unknown functional OBPs and CSPs of S. exigua, which will help us define the odorant binding spectrum of each gene. This will provide potential behavioral disturbance agents to control the moths by using reverse chemical ecology methods .  Table S1. This tree was constructed using PhyML based on alignment results of ClustalX.

ORs
Sixty-four different unigenes encoding putative ORs were identified by analyzing the transcriptome data of S. exigua, of which 57 were newly obtained ( Table 3). A total of 28 out of 64 SexiORs contained full-length ORFs that encoded 351 to 473 amino acids with various transmembrane domains (TMD). The phylogenetic analysis showed that all 64 SexiORs were clustered in an OR tree with B. mori, D. plexippus, and H. armigera, with each clustering having at least one other moth ortholog (Figure 7). In accordance with previous studies (Liu et al., 2013), we also identified a chaperone and higher conserved insect OR-SexiOrco (Krieger et al., 2005;Nakagawa et al., 2005;Xu and Leal, 2013;Missbach et al., 2014) and four pheromone receptors (SexiOR6, 11, 13, and 16) (Table 3, Figure 7), which suggests that our sequencing and analysis methods were reliable. The results of the phylogenetic and sequence homology analyses showed that we were able to obtain the fifth PR gene of S. exigua, SexiOR59. Liu's research (Liu et al., 2013) found that only two PRs (SexiOR13 and SexiOR16) showed higher electrophysiological responses to the three sex pheromone components (Z9, E12-14:OAc, Z9-14:OAc, and Z9-14:OH) of S. exigua; however, no PRs displayed specific or higher response to the fourth pheromone component Z9, E12-14:OH. Therefore, further studies are required to determine whether SexiOR59 can respond highly or not to Z9, E12-14:OH or other pheromone components. Additionally, other researchers have found that several non-PR ORs could respond to host plant volatiles, such as SlitOR12 of S. litura (Zhang et al., 2015c), EpstOR1, and three from Epiphyas postvittana (Jordan et al., 2009). Therefore, some ORs of the 58 non-PR ORs in S. exigua might play a similar role in the chemosensation of the volatiles in host plants.
To better infer the potential functions of these SexiGRs, we applied the qPCR method to investigate the expression profiles of all SexiGRs in six chemosensory organs (FA, MA, FPr, MPr, FLP, and MLP) and two non-chemosensory organs (Female abdomen, FAb and Male abdomen, MAb) (Figure 10).  Table S1. This tree was constructed using PhyML based on alignment results of ClustalX. The results showed that the organ with the highest SexiGRs expression was FPr (28 genes), followed by MPr (27 genes), FLP (25 genes), and MLP (22 genes), indicating that SexiGRs mainly exist within the gustatory organs, not the olfactory or non-chemosensory organs. This explains why the numbers of GR based on the antennae or non-gustatory organs transcriptome of other insects (Liu et al., 2014b;Xu et al., 2015) are lower than the SexiGRs in the present study. Additionally, we found 4, 16, 11, and 1 SexiGR genes that were highly expressed in the antennae, proboscises, labial palps, and abdomen of S. exigua, respectively, and some genes also showed differences in sex expression, which suggests that SexiGRs not only plays a pivotal role in gustatory processes (Jiang et al., 2015;Poudel et al., 2015), but might also be involved in olfactory (Agnihotri et al., 2016;Poudel et al., 2017) and other physiology processes (Xu et al., 2012;Ni et al., 2013). These results indicate that the proboscises and labial palps play more important roles in the taste perception process of than the olfactory organs do, which provides an important reference for future study of the taste perception mechanism in S. exigua as well as in other moths.
In conclusion, 159 genes encoding putative chemosensory genes were obtained by analyzing six chemosensory organs of S. exigua. Our approach proved to be highly effective for the identification of chemosensory genes in S. exigua, for which genomic data are currently unavailable. As the first step toward understanding gene functions, we conducted a comprehensive phylogenetic analysis of these genes and investigated all SexiGRs expression patterns, most of which were highly expressed in gustatory organs. The present study greatly improves the gene inventory for S. exigua and provides a foundation for future functional analyses of these crucial genes.

AUTHOR CONTRIBUTIONS
Y-NZ conceived and designed the experimental plan. Y-NZ, J-LQ, M-YL, and X-XX performed the experiment. Y-NZ, X-YZ, TX, and LS processed and analyzed the experiment data. J-WX and C-XL provided important suggestions to help modify the manuscript. Y-NZ wrote the manuscript.