In 2015, we collected about 200 D. punctatus pupae in Quanzhou, Guilin City, Guangxi province, China, and reared them in our laboratory at 26 ± 2°C, 50 ± 10% relative humidity, and a 16 h light: 8 h dark photoperiod. Male and female insects representing four different physiological conditions were prepared for transcriptome sequencing as follows: newly emerged (within 5 h after emergence, unmated; eclosion), calling females and corresponding males, mating status (copulating), and after mating status (post-coitum). Male and female insects were kept in two different insect rearing cages placed in close proximity to each other and separated by only two layers of screen cloth, so that males could sense the female sex pheromones. Antennae from 15 female and male D. punctatus specimens at each stage were cut off and immediately frozen in liquid nitrogen. Insect antennae from each group were divided into three equal parts, as three biological replicates. Thus, in total, we constructed 24 libraries for RNA-seq (four conditions for male and female insects respectively, and three replications for each status).

As previously described (Zhang et al., 2014a, 2017), total RNA samples were extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and treated with RNase-free DNase I (TaKaRa, Dalian, Liaoning, China). Subsequently, RNA purity, concentration, and integrity were checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA), a Qubit® RNA Assay Kit and a Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and the RNA Nano 6000 Assay Kit on the Bioanalyzer 2100 system (Agilent Technologies, CA, USA), respectively.

Duplex-specific-nuclease normalized cDNA was synthetized using 3 μg total RNA samples (Zhulidov et al., 2004; Bogdanova et al., 2008). The RIN values of all samples were > 8. We prepared sequencing libraries using an Illumina TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA, USA), and added four index codes to identify sequences from each sample. To preferentially select cDNA fragments of 200 bp length, we purified the libraries using the AMPure XP system (Beckman Coulter, Beverly, MA, USA). PCR (10 cycles) was performed to enrich for the two-end adaptor ligated DNA fragments. Finally, the products were purified using an AMPure XP system and quantified on an Agilent Bioanalyzer 2100.

Index-coded samples were clustered using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) on a cBot Cluster Generation System, then sequencing performed on an Illumina Hiseq 2500 platform, according to the manufacturer's instructions.

Raw sequencing data were filtered to remove reads containing adapter sequence, reads with > 10% N (uncertain bases), and sequences with error rates > 1% for more than 50%, using self-written Perl scripts, to obtain clean data. Then we calculated the Q20, Q30, GC-content, and sequence duplication level of the clean data, which were subsequently used for downstream analyses. Clean data sequences were compared with the NT database to determine whether they were polluted. Trinity (vesion:trinityrnaseq_r20131110) was used to perform transcriptome assembly (Grabherr et al., 2011). TGICL software was used to reduce redundancy (Pertea et al., 2003). The raw data from our experiments have been deposited in the NCBI SRA database under the accession number SRP102206 (Bioproject accession number PRJNA374901). We assessed the transcriptome assembly using benchmarking universal single-copy orthologs (BUSCO) based on the percentage of sequences aligned with highly conserved protein sequences, (Simão et al., 2015).

First, transcript sequences were searched using BLAST against the NR, SWISSPROT, KEGG, and KOG databases, with a cut-off value of 1e-5, and the highest sequence similarity targets selected for functional annotation of the transcripts. Next, Blast2GO was used to perform GO annotation of the transcripts (Conesa et al., 2005; Götz et al., 2008). Finally, the molecular function, biological process, and cellular component of the genes were assigned (Ashburner et al., 2000; Krieger et al., 2004).

Based on our previous research (Zhang et al., 2017), we further identified the chemosensory genes in D. punctatus. Previously identified chemosensory genes were confirmed in our new transcriptome database using tBLASTx searches, and the complete sequences of some previously identified partial genes obtained. We further identified some new chemosensory genes by contig tBLASTx searches. The open reading frames (ORFs) of possible genes were verified by additional BLAST searches (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Newly identified olfactory genes were submitted to GenBank; the updated accession numbers are listed in Table S1. Maximum likelihood (ML) and neighbor-joining (NJ) phylogenetic trees of chemosensory genes were constructed using MEGA5 with 1,000 bootstrap replications (Tamura et al., 2011). MEGA's model test was used to select the best model for ML tree construction. Dendrograms were colored in Adobe illustrator (Adobe Systems). Motif analysis of the predicted intact ORFs of chemosensory genes was performed using the MEME online server (Version 4.12.0.) (Bailey et al., 2009) http://meme-suite.org/tools/meme. For OBP and CSP, the motif discovery parameters were: minimum width = 6, maximum = 10, maximum motifs to find = 8; for ORs they were: minimum width = 15, maximum = 50, maximum motifs to find = 8.

To measure the gene expression levels in transcriptomes, we used the FPKM (fragments per kilobase of exon per million fragments mapped) criteria (Trapnell et al., 2010). Three biological replicates were sequenced for each D. punctatus status, and the mean FPKM value and standard error obtained from the three replicates. Differentially expressed genes (DEGs) between different mating status insects were calculated using DESeq (http://bioconductor.org/packages/release/bioc/html/DESeq.html) (Anders and Huber, 2010), based on the reads of each unigene. Unigene expression levels and DEGs were normalized following the compatible-hits-norm model (Bullard et al., 2010). DEGs were screened to identify those generating q-values ≤ 0.05 using the false discovery rate (FDR) method (Noble, 2009).

GO Enrichment analysis of DEGs was carried out using GOstat (Beißbarth and Speed, 2004), with p-values approximated using Chi-square tests, with all annotated genes used as the background. Short Time-series Expression Miner (STEM, vision 1.3.11) was used to analyze the expression trends of some DEGs (Ernst and Bar-Joseph, 2006). FPKM values were log2 transformed and imported into the software.

qPCR was carried out to validate the RNA-Seq data, similar to our previous report (Zhang et al., 2014a,b, 2017). qPCR primers (Table S2) were designed based on cDNA sequences. RT-PCR was performed to test whether qPCR primers could amplify the correct products. Beta-actin was used as the housekeeping gene. T-easy clones containing the tested genes were constructed as reference genes to construct qPCR standard curves. Amplification efficiencies of all primers tested were 90–100%. Real-time PCR was carried out in a Roche LightCycler 480 (Stratagene, La Jolla, CA, USA). The PCR cycles were as follows: 2 min at 95°C, 40 cycles of 20 s at 95°C, 20 s at 58°C, and 20 s at 72°C; finally, melting curve analysis (58 to 95°C) was performed to evaluate the specificity of the PCR products. Ct values were calculated using the Roche qPCR software (version 1.5.1) with the second derivative method. Three independent biological reactions were completed for each insect status, along with three technical replicates for each reaction. Gene expression levels tested by qPCR in female and male D. punctatus antennae (relative to that of the actin gene) were compared with the transcriptome expression data (FPKM), as illustrated in Figure S1.

Transcriptomic sequence data were generated from antenna cDNA libraries from D. punctatus adults at different mating stages using Illumina HiSeq™2500 technology. We acquired 204.30 Gbp of clean sequence data in 1,634,361,960 clean reads. After assembly, 110,760 unigenes were obtained, with an N50 of 2,380 bp. Approximately 80% of unigenes were >500 bp, with a maximum length of 54,680 bp (Figure S2). We evaluated the completeness and accuracy of transcriptome assembly using BUSCOs, and the results demonstrated 98.1% complete BUSCOs (C: 98.1 [S: 83.8%, D: 14.3%], F: 1.0%, M: 0.9%).

BLAST analysis indicated that D. punctatus transcriptome unigenes were most similar to amino acid sequences from three other Lepidoptera species: B. mori (16,769 hits with E-values <1e-5), Danaus plexippus (7826 hits with E-values <1e-5), and P. xylostella (7034 hits with E-values <1e-5) (Figure 1). These three species accounted for ~75% of hits.

Function distribution, determined by Gene Ontology (GO) analysis, indicated that D. punctatus genes were primarily enriched for binding or catalytic activity (Figure S3A), and followed by transporter and structural molecule activity. KOG classification indicated that genes involved in signal transduction occupied an important position (Figure S3B).

GO enrichment analysis of DEGs revealed that waves of gene expression changes occurred in the antennae of D. punctatus during the mating process (Figure S4). Overall, olfactory detection genes, particularly olfactory receptors and odorant binding related genes, exhibited dramatic differences between sexes and also during the mating process of male and female insects. To further elucidate the characteristics of gene expression in insects of different mating status, we analyzed the expression trend of DEGs with short time-series expression miner (STEM). Thirteen significant profiles were obtained, four of which were related to chemosensory genes (Figure 2A). GO enrichment of the profiles indicated that profiles 41 and 26 included the most chemosensory associated genes (Figure 2B). Of the four chemosensory related profiles, profile 7 continually declined and included genes expressed more highly in female antennae; profile 41 continually rose, and included genes expressed at higher levels in male antennae; profiles 26 and 27 fluctuated in the eight D. punctatus groups, and included genes with expression levels that oscillated during the mating process of this insect.

GO enrichment and STEM analysis indicated that chemosensory genes may be very important in the mating process. Thus, detailed analyses were performed to determine the characteristics of the olfactory-related gene families identified from the transcriptomes of D. punctatus in different mating states.

In our previous work, we identified a considerable number of D. punctatus chemosensory genes (Zhang et al., 2017). Here, after further effort, the complete sequences of many of the partial gene sequences identified previously were acquired, including six OBPs (NCBI accession numbers KY225481–KY225486), one CSP (KY225487), 23 ORs (KY225488–KY225510), one GR (KY225519), and eight IRs (KY225529–KY225536). Some new genes were also identified, including eight ORs (KY225511–KY225518), nine GRs (KY225520–KY225528), and five IRs (KY225537–KY225541).

The correlation between the expression levels of chemosensory genes and mating status was examined in detail, and the expression levels of OBPs (Figure S5), CSPs (Figure S6), ORs (Figure S7), GRs (Figure S8), and IRs (Figure S9) determined. Further analysis indicated that the chemosensory genes exhibited different expression levels in insects in different mating states, with three different patterns identified (Figure 3). First (Type I), some genes were more strongly expressed in male than female antennae. In general, these genes were upregulated during calling or mating, and downregulated after mating. There were nine OBPs (Figure 3, Figure S5A); ten CSPs (Figure S6A); and nine ORs (Figure S7A) in this category. Second (Type II) were genes expressed at higher levels in female than male antennae. These genes were generally upregulated during calling or mating, and downregulated after mating (Figure 3). This category consisted of four OBPs (indicated as red in Figure S5B), two CSPs (Figure S6B), and six ORs (Figure S7B). Third (Type III) were genes expressed at higher levels in female than male antennae, and continually upregulated over time (Figure S3). This third category contained three OBPs (green in Figure S5B) and nine ORs (green in Figure S7B). To further confirm the expression level of olfactory genes, the ORs (Figure S1A) and OBPs (Figure S1B) that were belong to type I were selected for qPCR verification of the transcriptome expression data, and the results indicated that the transcriptome expression data were credible (Figure S1).

Next, phylogenetic analyses of chemosensory genes identified from D. punctatus were performed. A phylogenetic tree of the identified OBPs was constructed (Figure 4A). Unsurprisingly, PBP1, PBP2, and two GOBPs were grouped together, and three of these four genes were expressed at higher levels in male than female antennae. Interestingly, other OBPs more strongly expressed in male antennae were all grouped with an OBP that was preferentially expressed in female antennae. Phylogenetic analysis indicated that CSPs expressed at higher levels in male antennae (indicated by solid circles) were dispersed into two subclasses in the tree, while CSPs expressed more strongly in female antennae (filled triangles) were grouped separately (Figure 4B). Interestingly, the male biased ORs were almost all clustered in a single branch (Figure 4C).

To further analyze the characteristics of D. punctatus ORs, we performed phylogenetic analysis including ORs from two sister species of D. punctatus, Dendrolimus houi, and Dendrolimus kikuchii (Zhang et al., 2014a), and four other Lepidopteran species, including B. mori, D. plexippus, M. sexta (Grosse-Wilde et al., 2011), and Cydia pomonella (Bengtsson et al., 2012). The results permitted several observations (Figure 5): first, the co-receptor Orco was identified in D. punctatus and was conserved among these moths; second, the ORs from D. punctatus generally formed small subgroups together with those of D. houi and D. kikuchii, and sometimes with B. mori and M. sexta; third, the sex pheromone receptors from B. mori, M. sexta, D. plexippus, and C. pomonella formed clade in the tree (labeled “sex pheromone receptors” in Figure 5); however, none of the ORs from the three Dendrolimus species were clustered in this group; finally, a group of ORs from D. punctatus, D. houi, and D. kikuchii (labeled “Dendrolimus Specific Odorant Receptors” in Figure 5) formed a subgroup that included no receptors from the other moths, this is unusual in the Lepidoptera, and the specific functions of these ORs require further investigation.

To further understand the sequence characteristics of the chemosensory genes in D. punctatus, we performed motif-pattern analysis of OBPs, CSPs, and ORs using the MEME server. OBP motif analysis revealed eight groups (Figure 6). Motif 1 was contained in all OBPs. All GOBPs and PBPs shared four motifs 1–4, while the eighth motif was exclusive to PBPs, and the seventh motif was only found in GOBPs. OBPs in the third and fourth groups that contain the fifth motif were all minus-C OBPs, indicating that Motif 5 is a characteristic of minus-C OBPs. The eighth group included two plus-C OBPs.

Motif analysis of CSP sequences indicated that they were relatively conserved, with the majority containing the same motif pattern, with some exceptions (Figure 7). Interestingly, two CSPs (CSP2 and CSP15) that contained different motif patterns with respect to the others were those expressed at higher levels in female antennae.

OR motif analysis indicated that the majority of sequences could be separated into nine groups (Figure 8), depending on their motif patterns. We designated the first five groups as class 1, and the sixth to ninth groups as class 2, as the conserved motifs of class 1 were concentrated at the 5′ end of the genes, while those of class 2 were at the 3′ end. Comparative analyses indicated that male biased ORs (Figure 3, Type I) all belong to class 1, except for OR51, while female bias ORs (Figure 3, Types II and III) all belong to class 2, other than OR62. Overall, the motif patterns and expression biases of these genes indicated their functional differentiation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.