Larvae of M. separata were collected in Xinxiang, Henan Province, China. The colony was reared on an artificial diet in the laboratory and maintained under the conditions of 27 ± 1°C, 75 ± 5% relative humidity, and a 14-h/10-h light/dark cycle. Pupae of different sexes were separated in glass Petri dishes before eclosion. Adult moths were provided with sucrose solution 10% (v/v). Brains, antennae, wings, legs and abdomens of unmated moths were collected 2–4 days after eclosion, immediately frozen in liquid nitrogen, and stored at −70°C for RNA extraction.

Total RNA extracted from brains of approximately 600 adult males and females independently were used to construct separately three female and three male cDNA libraries. The libraries were sequenced using the PE100 strategy on the Illumina HiSeq™ 2000 platform (Illumina, San Diego, CA, United States) at Novogene Bioinformatics Technology Co. Ltd. (Beijing, China). Briefly, mRNA was purified from total RNA using magnetic beads with Oligo (dT), and then, it was fragmented into short fragments after adding fragmentation buffer. First-strand cDNA was synthesized using random hexamer primer and M-MuLV reverse transcriptase (RNase H). Subsequently, the second-strand cDNA was synthesized using DNA polymerase I and RNase H. The double-stranded cDNA was purified using the AMPure XP system (Beckman Coulter, Beverly, MA, United States). NEBNext Adaptors having a hairpin loop structure were ligated to prepare for hybridization after the adenylation of the DNA fragments’ 3′ ends. Library fragments were purified using the AMPure XP system for selecting preferentially cDNA fragments of 150–200 bp. Then, the selected fragments were used as templates for PCR amplification. PCR products were also purified using the AMPure XP system, and library quality was assessed on an Agilent 2100 Bioanalyzer and with a Q-PCR system.

A de novo transcriptome was assembled using the paired-reads mode with default parameters using the short-read program Trinity (Grabherr et al., 2011). Trinity outputs were clustered using TGICL (Pertea et al., 2003). The consensus cluster sequences and singletons made up the final unigene dataset. The generation of unigenes was performed using BLASTx and BLASTn programs against the public databases, with an E-value threshold of 10^–5. GO terms were extracted from the best hits obtained from BLASTx against the NR database using the Blast2GO program (Conesa et al., 2005). A GO functional classification of all the unigenes was performed using WEGO software.¹ KOG and KEGG annotations were performed using Blastall software against the KOG² and KEGG³ databases, respectively.

Candidate unigenes encoding putative OBPs and CSPs were selected on the basis of the NR annotation results in the remote server. All the candidate chemosensory genes were further manually checked using the BLASTx program. The open reading frame (ORF) of each candidate unigene was predicted using the ORF finder tool.⁴ The putative signal peptides of OBP and CSP protein sequences were predicted using SignalP 4.1.⁵ In addition, all the candidate genes were compared with previously reported sequences using the BLASTn program (with an E-value threshold of 10^–5) to identify novel OBPs and CSPs genes (Du et al., 2018). These genes were named in accordance with gene naming rules by adding a suffix with a number to indicate the descending order of their coding region lengths.

Multiple alignments of amino acid sequences were performed using the online prediction website MAFFT.⁶ The phylogenetic trees were constructed using the maximum-likelihood method with a bootstrap analysis of 1,000 replicates and the JTT with Freqs. (+F) Substitution Model using MEGA5.2 (Tamura et al., 2011). The phylogenetic trees were visualized using FigTree v1.4.3.⁷ OBPs data sets contained 29 candidate OBPs from M. separata, and 150 from other Lepidopteran moths, including B. mori (Gong et al., 2009), H. armigera (Liu et al., 2012), Helicoverpa assulta (Chang H. et al., 2017), Spodoptera exigua (Liu et al., 2015), Heliothis virescens (Vogel et al., 2010), and Spodoptera litura (Gu et al., 2015). CSPs data sets contained 16 putative CSPs from M. separata, and 72 from other Lepidopteran moths, including B. mori (Gong et al., 2009), H. assulta (Chang H. et al., 2017), H. armigera (Zhang J. et al., 2015), H. virescens (Picimbon et al., 2001), Agrotis ipsion (Gu et al., 2014), S. litura (Zhang Y. N. et al., 2015), and S. exigua (Liu et al., 2015). The amino acid sequences used in the phylogenetic analyses are listed in Supplementary Materials 1, 2.

To confirm the expression profiles of the identified OBPs and CSPs genes, semi-quantitative PCR (RT-PCR) was performed. Total RNA was isolated from brains, antennae, wings, legs, and abdomens of 50–60 adults and extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s instructions. Single-stranded cDNA templates were synthesized from 1 μg of total RNA from various tissue samples using the FastKing gDNA Dispelling RT SuperMix (TianGen, Beijing, China). Specific primers of predicted OBPs and CSPs genes were designed using Premier 5.0 (Supplementary Material 3 and Supplementary Table 1). PCR reactions were carried out using equal amounts of cDNA (200 ng) template. The β-actin (GenBank Acc. GQ856238.1) of M. separata was selected as the reference gene to test the integrity of the cDNA templates and also the expression quantification of the target genes. The PCR was performed in a Mastercycler^® (Eppendorf, Hamburg, Germany) under the following conditions: 94°C for 5 min, 25–33 cycles (depending on the expression level of each gene) of 94°C for 30 s, 56°C for 30 s, and 72°C for 30 s, followed by a final extension at 72°C for 10 min. PCR products were analyzed on 1.0% agarose gels and visualized after staining with ethidium bromide.

The RT-qPCR analysis was conducted using an ABI QuantStudio3 (Applied Biosystems, Carlsbad, CA, United States). The specific RT-qPCR primers were designed using Beacon Designer 8.13 (PREMIER Biosoft International, CA, United States) (Supplementary Material 3 and Supplementary Table 2). Two reference genes, β-actin (GenBank Acc. GQ856238.1) and glyceraldehyde-3-phosphate dehydrogenase (gapdh) (GenBank Acc. HM055756.1) were used to normalize target gene expression. The amplification efficiencies of the target and reference gene primers were evaluated using a four-fold serial dilution of cDNA templates from adult antennae. Reactions for each sample (20 μl) consisted of 10 μl of SuperReal PreMix Plus (TianGen, Beijing, China), 0.5 μl of each primer (10 μM), 0.4 μl of Rox reference dye, 1 μl of sample cDNA, and 7.6 μl of sterilized ultrapure water. Amplification conditions were an initial denaturation at 95°C for 3 min, followed by 40 cycles of 95°C for 15 s, and a single step for annealing and extension was performed at 60°C for 30 s. The PCR products were heated to 95°C for 15 s, cooled to 60°C for 1 min, heated to 95°C for 30 s, and cooled to 60°C for 15 s to determine the dissociation curves. The RT-qPCR reaction of each sample was performed in three technical replicates and three biological replicates. Then, we used the relative quantitation method (2^–ΔΔCT) (Livak and Schmittgen, 2001) to evaluate quantitative variation. Transcript amounts were standardized to 1 using the sample from adult male brain. Data were analyzed using Data Processing System software version 9.5 (Tang and Zhang, 2013). A one-way analysis of variance with Tukey’s multiple comparison test was performed to analyze differences in gene expression levels among multiple samples, and p < 0.05 was considered to be significant.

Six adult brain cDNA libraries, three for females and three for males of M. separata, were constructed and sequenced using the Illumina HiSeq™ 2000 platform. As a result, 61,283,994, 59,876,048, and 68,196,054 raw reads were produced from the three separate female brain samples; and 46,814,292, 60,931,890, and 43,177,818 raw reads were produced from the three separate male brain samples. After trimming the adaptor sequences, contaminating sequences, and low quality sequences, 56,689,466, 54,841,746, and 62,601,386 clean reads of the three separate female brain samples, and 43,297,968, 56,087,138, and 40,004,380 clean reads of the three separate male brain samples, remained for the following assembly (Supplementary Material 3 and Supplementary Table 3). Subsequently, all the clean reads were assembled together and generated 132,516 unigenes with lengths ranging from 201 to 28,894 bp, with a mean length of 579 bp.

Homology searches querying the 132,516 unigenes against other insect species were performed using the BLASTx and BLASTn programs, with the E-value cut-off of 1.0E^–5. In total, 27,594 unigenes (20.82%) had BLASTx hits in the NR database, and 12,499 unigenes (9.43%) had BLASTn hits in the NT database. Among the annotated unigenes, 4445 (3.35%) were annotated in all of the databases [NR, NT, KO (KEGG ontology), SwissProt, protein family (PFAM), GO, and KOG], whereas 35,484 (26.77%) were annotated in at least one database (Supplementary Material 3 and Supplementary Table 4). The analysis showed that most M. separata protein sequences were orthologs of proteins in B. mori (33.5%), Danaus plexippus (15.5%), and Plutella xylostella (14.1%) (Figure 1).

According to the GO category analysis, only 21,188 (15.99%) assembled unigenes corresponded to different functional groups. Because one unigene can align to multiple GO categories, 54,623 (41.22%) unigenes were assigned to biological process, 33,526 (25.30%) to cellular component, and 23,545 (17.77%) to molecular function. In the molecular function category, the terms of binding and catalytic activity were the most represented. In the cellular component terms, cell and cell part were the most abundant. In the biological process category, cellular process, metabolic process, and single-organism process were most abundant (Figure 2A).

In the obtained KOG functional annotation, 12,672 unigenes were categorized into 26 functional groups (Figure 2B). “General function prediction only” was the largest group (2764, 21.81%), followed by the “Signal transduction mechanism (1594, 12.58%) and Posttranslational modification, protein turnover, chaperones” (1353, 10.68%) groups, and “Cell motility” (25, 0.20%) was the smallest group.

In the M. separata brain transcriptomes, 29 OBPs were annotated on the basis of the TBLASTN results. Among them, 16 OBPs contained intact ORFs, with lengths ranging from 133 to 334 amino acids (Table 1). Based on the numbers and locations of the conserved cysteines, the OBPs were classified into three categories, Classic, Pluc-C, and Minus-C OBPs families. Seven full-length OBPs (GOBP2, OBP6, OBP9, OBP11–13, and OBP15) had the typical six conserved cysteines and spacing, forming the Classic OBPs family. Three full-length OBPs (OBP4, -5, and -10) belonged to the Pluc-C OBPs family, having additional two, three, and six cysteines located downstream of the conserved C6. The remaining six full-length OBPs (OBP1–3, -8, -14, and -16) belonged to the Minus-C OBPs family. OBP1 lacked conserved cysteine C1; OBP2 had none of the typical six conserved cysteines; OBP3 only had conserved cysteines C1 and C6; and OBP8, -14, and -16 lacked the conserved cysteines C2 and C5 (Figure 3). Compared with our earlier identified OBPs in M. separata (Du et al., 2018), many of them shared high sequence identity levels, ranging from 80 to 100%. However, OBP2, -3, -17, and -23 shared a no more than 36% sequence identity (Supplementary Material 3 and Supplementary Table 5).

Sixteen transcripts encoding candidate CSPs were identified. Among them, 14 CSPs contained intact ORFs, with lengths ranging from 106 to 290 amino acids (Table 2). These identified full-length CSPs proteins included a signal peptide and four highly conserved cysteine profiles (C1-X6-C2-X18-C3-X2-C4, where X represents any amino acid) (Figure 4). The CSP15 and CSP16 amino acid sequences were incomplete due to the lack of a 3′ or 5′ terminus. Compared with previously identified CSPs (Du et al., 2018), most of them shared high sequence identity levels, ranging from 74 to 100%. However, there were two CSPs, CSP9 and CSP12, that shared less than a 41% sequence identity (Supplementary Material 3 and Supplementary Table 5).

All the putative OBPs clustered with at least one lepidopteran ortholog in the phylogenetic tree (Figure 5). The identified GOBP and PBP protein sequences in brain clustered into the GOBP and PBP clades in the phylogenetic tree, respectively. OBP5 and OBP7 shared a 68.28% sequences identity and clustered together. OBP11 and OBP18 only shared a 29.66% similarity, but they also clustered together. With OBP13, they formed a clade containing SlitOBP25. All the candidate CSPs proteins clustered with at least one lepidopteran ortholog in the phylogenetic tree (Figure 6). CSP3 and CSP15 shared a 67.97% sequences identity, and formed a clade with AipsCSP7. CSP7 and CSP16 shared a 45.60% sequences identity and clustered together.

The RT-PCR expression profiles indicated that the majority of OBPs genes were expressed in the antennae. OBP5, -7, -11, -13, -16, -18, -21, and -24 were detected in the brain. Among them, OBP5, -7, -11, and -13 showed male brain-biased expression, whereas OBP21 showed female brain-biased expression. OBP16, -18, and -24 were expressed in both female and male brains. OBP5, -7, and -11 were also detected in the antennae and abdomens. OBP13 could be detected in the abdomens. OBP16 and OBP24 were detected in the antennae, abdomens, and legs. OBP18 could be detected in all the tissues, and OBP21 was detected in the legs and wings (Figure 7). The RT-PCR expression profiles indicated that most CSPs genes were expressed in all the examined tissues (Figure 7). Three CSP genes, CSP2–4, were expressed in female brain. CSP8 was detected in male brain, and CSP10–12, and CSP15 were detected in both female and male brains (original gel images of RT-PCR can be found in Supplementary Material 4).

To confirm the RT-PCR results, RT-qPCR was performed to measure quantitatively the expression levels of six OBPs (OBP3–5, -9, -10, and -16) and 2 CSPs (CSP3 and CSP4) genes in the various tissues (Figure 8). The RT-qPCR results were mostly consistent with the RT-PCR results. They confirmed that OBP5, OBP16, CSP3, and CSP4 were expressed in brain, and further revealed that OBP10, which was not detected by RT-PCR, was also expressed in brain. The expression levels of OBP5 and OBP10 in brain were not significantly different between the sexes (p > 0.05), whereas that of OBP16 was three times higher in male brain than in female brain (p < 0.05). The expression levels of CSP3 and CSP4 were 40 and 24 times higher in female brain than in the male, respectively (p < 0.05).