The L. decemlineata adults were collected from potato fields in Xinjiang Province, China. Male and female adults were separated, not considering the ages or mating status. The antennae were pulled off with tweezers grasped at the very root of the antennae. The separated antennae were stored in RNAlater (Ambion, Austin, TX, USA) and taken to the Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing. After removing the residual RNAlater, the stored antennae were crushed with a vitreous homogenizer. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. The RNA was dissolved in RNase-free water and the integrity and quantity of RNA was determined by gel electrophoresis and Nanodrop ND-2000 spectrophotometer (NanoDrop products, Wilmington, DE, USA). Residual gDNA in total RNA was removed by DNase I (Promega, Madison, WI, USA) before cDNA library construction.

Five micrograms of total RNA extracted from approximately 200 antennae of adult male or female adults were sent to Beijing Genome Institute (Shenzhen, China) for construction of cDNA library and sequencing. Briefly, mRNA was isolated and fragmented into 200–700 nt pieces. Random hexamers were used for first-strand cDNA synthesis. Then the second-strand cDNA was synthesized using RNase H and DNA polymerase I. The resulting double-stranded cDNAs were treated with T4 DNA Polymerase and T4 Polynucleotide Kinase for end-repairing and dA-tailing. After that, they were ligated to sequencing adaptors with barcode using T4 DNA ligase. Finally, fragments with around 200 bp length were collected by 2% agarose gel electrophoresis and purified with QiaQuick GelPurify Kit (Qiagen, Hilden, Germany), and used as templates for PCR amplification. The libraries were pair-end sequenced using PE90 strategy on Illumina HiSeq™ 2000 (Illumina, San Diego, CA, USA) at the Beijing Genome Institute. The male and female libraries were sequenced in one lane then raw-reads were sorted out by barcodes.

Raw reads from each library were filtered to remove low quality reads and the sequence reads containing adapters and poly-A/T tails. The resulting clean reads were assembled to produce unigenes with the short reads assembling program-Trinity using the default parameters (Grabherr et al., 2011). Then the unigenes from the two samples were pooled together and clustered by TGI Clustering Tool (TGICL) (Pertea et al., 2003). The consensus cluster sequences and singletons make up the unigenes dataset.

The annotation of unigenes was performed by NCBI blastx against a pooled database of non-redundant (nr) and SwissProt protein sequences with an e-value cut-off of 1e-5 (Altschul et al., 1997). The blast results were then imported into the Blast2GO for GO Annotation (Conesa et al., 2005). Protein coding region prediction was performed by OrfPredictor (Min et al., 2005) according to the blast result. The signal peptide of the protein sequences were predicted using SignalP 4.0 (Petersen et al., 2011). The transmembrane-domains of annotated genes were predicted using TMHMM Server Version2.0 (http://www.cbs.dtu.dk/services/TMHMM) (Krogh et al., 2001).

The phylogenetic reconstruction implemented for the analysis of OR, IR, OBP, and CSP was performed based on the amino sequences of the candidate olfaction genes and the collected data sets. The OR data set contained OR sequences identified in Coleopteran species (239 from T. castaneum (Richards et al., 2008; Kim et al., 2010), 49 from D. ponderosae (Andersson et al., 2013), 42 from I. typographus (Andersson et al., 2013), and 57 from M. caryae (Mitchell et al., 2012). The IR data set contained 15, 7, and 66 IR sequences from D. ponderosae (Andersson et al., 2013), I. typographus (Andersson et al., 2013) and D.melanogaster (Croset et al., 2010), respectively. The OBP data set contained 46 sequences from T. castaneum (Richards et al., 2008; Kim et al., 2010), 31 sequences from D. ponderosae (Andersson et al., 2013), and 15 sequences from I. typographus (Andersson et al., 2013). The CSP data set contained the 40 sequences from T. castaneum (Richards et al., 2008; Kim et al., 2010), 11 sequences from D. ponderosae (Andersson et al., 2013), and 5 sequences from I. typographus (Andersson et al., 2013). The protein name and accession number of the genes used for phylogenetic tree building are listed in Supplementary Material S1. Amino acid sequences were aligned using MAFFT (https://www.ebi.ac.uk/Tools/msa/mafft/) (Katoh and Toh, 2008). Unrooted trees were constructed by the maximum-likelihood method in FastTree 2.1 software using the default parameters. To estimate reliability of each split in the tree, the local support values were computed based on the Shimodaira-Hasegawa (SH) test (Price et al., 2010) Dendrograms were created and colored in FigTree software (http://tree.bio.ed.ac.uk/software/figtree/).

To illustrate and compare the expression of candidate ORs in male and female antennae, semi-quantitative reverse transcription PCR (RT-PCR) was performed using cDNAs prepared from male and female antennae. L. decemlineata tissue samples were collected for three biological replicates. In each replicate, about two micrograms total RNA were extracted from approximately 100 antennae of male or female adults as mentioned above. Prior to cDNA synthesis, RNA were treated with DNase I to remove trace amounts of genomic DNA. The cDNA was synthesized by First Strand cDNA Synthesis Kit (Fermentas, Vilnius, Lithuania) and was used as a template in PCR reactions with gene-specific primers. Ribosomal protein L31 (LdecRL31) and ribosomal protein S3 (LdecRPS3) were used as controls. Primers were designed using the Primer Premier 5 software (PREMIER Biosoft International). The primer sequences are available in Supplementary Material S2. Taq MasterMix (CWBIO, Beijing, China) was used for PCR reactions under general 3-step amplification of 94°C for 30 s, 53°C for 30 s, 72°C for 30 s. The PCR cycle-numbers were adjusted respectively for each gene. For OR, cycle-numbers ranged from 38 to 40. For high-express-level control genes LdecRL31 and LdecRPS3, cycle-numbers were reduced to 28. PCR products were run on a 2% agarose gel and verified by DNA sequencing. In the negative control, the cDNA template was replaced by water.

Using an Illumina HiSeq 2000 90PE RNA-Seq strategy, a total of 56.75 million and 59.19 million raw-reads were obtained respectively from the libraries of male and female antenna. After removing low quality and adaptor reads, 51.43 million and 52.36 million clean-reads were generated. The total bases of sequence data were approximately 4.63 and 4.71 gigabases from male and female samples, respectively. The clean reads of the L. decemlineata antennal transcriptome were deposited in the NCBI SRA database, under the accession number of SRX974484 (male) and SRX974488 (female). The clean-reads were assembled into 47,808 and 50,605 unigenes separately for male and female. All unigenes were merged and clustered into the final 45,179 unigenes consisting of 12,483 distinct clusters and 32, 696 distinct singletons. The transcript dataset was 32.46 megabases in size and with a mean length of 718 nt and N50 of 1, 116 nt. 10,756 unigenes were larger than 1000 nt in length, which comprised 23.81% of all unigenes (Table 1).

The functional annotations of the unigenes were performed mainly based on the blastx results against the nr database. Through annotation by blastx, 24,880 (55.1%) unigenes matched to known proteins. Among the 24,880 annotated unigenes, 14,463 (58.1%) showed strong homology (e-value smaller than 1e-45), whereas 5897 (23.7%) showed poor matches with e-value between 1e-15 and 1e-5 (Figure 1A). The similarity comparison showed 12,089 (48.6%) unigenes have more than 60% similarity with known proteins (Figure 1B). Blast analysis showed that 70.4% of the annotated unigenes matched with T. castaneum, followed by D. ponderosae (3.8%), Acyrthosiphon pisum (2.6%) and Camponotus floridanus (1.4%) (Figure 1C).

Gene ontology (GO) annotation of the unigene set was obtained using Blast2GO pipeline according to the blastx search against nr. From the 45,179 final unigenes set, a total of 11,704 unigenes were assigned various GO terms. In the molecular function category, the genes expressed in the antennae were mostly enriched to binding activity (e.g., nucleotide, ion, and odorant binding) and catalytic activity (e.g., hydrolase and oxidoreductase). In the biological process terms, cellular, and metabolic processes were the most represented. In the cellular component terms, cell, cell part, and organelle were the most abundant (Figure 2).

Within the L. decemlineata antennal transcriptome, 26 different sequences encoding odorant binding proteins were identified. Sequence analysis identified all but four transcripts (LdecOBP3, LdecOBP4, LdecOBP16, and LdecOBP26) with a full length ORF. The signal peptide, which is a typical structure of OBPs was not found in only two LdecOBPs (LdecOBP3 and LdecOBP26), due to incomplete N-termini. The length of all full-length LdecOBPs ranged from 122 to 255 amino acids. Compared to the ORs, insect OBPs are more highly conserved. The similarity between the LdecOBPs and known OBP of other insects was relatively low. Only seven predicted OBPs (LdecOBP4, LdecOBP9, LdecOBP10, LdecOBP12, LdecOBP13, LdecOBP20, and LdecOBP26) have more than 50% similarity with OBPs from T. castaneum or Batocera horsfieldi (Table 2). In our phylogenetic analysis of the OBPs in different beetles, LdecOBPs are spread across various branches (Figure 3) where they generally formed small subgroups together with OBPs from other three beetles. These splits were strongly supported by high local support values. A species specific branch consisting of 5 OBPs from L. decemlineata (LdecOBP3, LdecOBP6, LdecOBP7, LdecOBP8, and LdecOBP11) that is divergent from OBPs of other insects has been identified; these specific LdecOBPs might have some key species specific function.

The information, including unigene reference, length, and best blastx hit of all the 26 LdecOBPs are listed in Table 2. The sequences of all 26 LdecOBPs are listed in Supplementary Material S3.

Bioinformatic analysis led to the identification 15 different sequences encoding candidate CSPs in the L. decemlineata antennal transcriptome. Sequence analysis identified ten unigenes with a full length ORF with a predicted signal peptide sequence (Table 3).

Compared to OBPs, the conservation of CSPs of different Coleopteran was relatively high. Two thirds (10) of the LdecCSPs had more than 50% similarities with other CSPs (Table 3). The phylogenetic analyses also indicated conservation of Coleopteran CSPs (Figure 4). Most candidate LdecCSPs clustered with orthologs of T. castaneum, D. ponderosae and I. typographus into a separate clade. Only 2 LdecCSPs (LdecCSP2 and LdecCSP3) formed one small subgroup. Only one sequence-LdecCSP15 had low local support value unable to clearly demonstrate their phylogenetic positions.

The information, including unigene reference, length, and best blastx hit of all the LdecCSPs are listed in Table 3. The sequences of all 15 LdecCSPs are listed in Supplementary Material S3.

We found three SNMPs (LdecSNMP1-3) in our transcriptome. Two of them were predicted to have full-length ORF. Both LdecSNMP1 and LdecSNMP2 had more than 50% (51 and 59%) identity with SNMP of D. ponderosae and T. castaneum. LdecSNMP3 had only 40% similarity with SNMP2 of T. castaneum (Table 4).

The information, including unigene reference, length, and best blastx hit of all the three SNMPs are listed in Table 4. The sequences of all three SNMPs were listed in Supplementary Material S3.

The unigenes related to candidate OR were identified by keyword search of the blastx annotation. We identified 37 distinct unigenes that were putative OR genes. Of these, a full-length LdecOrco gene coding 479 amino acids was easily identified because it had intact open reading frames and seven transmembrane domains, which are characteristic of typical insect ORs. The 36 predicted incomplete ORs were of short length and only three of them contained a deduced protein longer than 300 amino acids. The deduced protein length of 24 ORs were even shorter than 200 amino acids.

The blastx results showed that the identities of these predicted ORs with known insect ORs is quite low. Only six predicted ORs (LdecOrco, LdecOR6, LdecOR8, LdecOR10, LdecOR33, and LdecOR36) have greater than 50% identity with ORs from T. castaneum. Even the LdecOrco had only 86% identity with the Orco from T. castaneum. Phylogenetic analysis was performed with ORs from T. castaneum, D. ponderosae, I. typographus and M. caryae. The results once again suggest high divergence of the OR genes (Figure 5). The branch of Orco was easily detected as it has a high degree of identity. All of the other LdecORs were distributed in different branches of the phylogenetic tree. A species-specific branch was identified consisting of four ORs from L. decemlineata (LdecOR17, LdecOR22, LdecOR25, and LdecOR31) that was clearly divergent from other ORs. Four LdecORs (LdecOR16, LdecOR18, LdecOR23, and LdecOR30) showed close relation to OR167 from T. castaneum, and these five ORs formed a distinct subgroup. Most of the splits in the tree were supported by high local support values and only a few splits were not reliable.

Information, including unigene reference, length, and best blastx hit of all 37 OR are listed in Table 5. The sequences are listed in Supplementary Material S3.

The putative IR genes in the L. decemlineata antennal transcriptome were represented according to their similarity to known insect IRs. Bioinformatic analysis led to the identification of ten candidate IRs, all ten sequences are marked as incomplete due to lacking a complete 5′ or 3′ terminus. The insect IRs contained three transmembrane domains (Benton et al., 2009). TMHMM2.0 predicted nine candidate IRs with different numbers of transmembrane domains (Table 6). One candidate IR was deemed to be an IR8a homolog due to its high identity (59%) to DponIR8a. A candidate IR25a homolog was also easily identified. The subgroup of IR75q2 is likely to extend to L. decemlineata, as four transcripts had high identity to IR75q2 homologs from C. pomonella, S. littoralis, and Aedes aegypti. Two IR76b homologs (LdecIR76b.1 and LdecIR76b.2) were also detected. The remaining two LdecIRs have similarity with IR87a and IR93a from D. melanogaster, respectively. In the phylogenetic tree of IRs, all L. decemlineata IR candidates clustered with their ionotropic receptor orthologs into separate sub-clades (Figure 6). Because of the relative high conservation of IRs, all the splits of LdecIRs were strongly supported by high local support values. The information, including unigene reference, length, and best blastx hit of all the ten IRs are listed in Table 6. The sequences of all 20 IRs were listed in Supplementary Material S3.

The expression patterns of the candidate 37 ORs in male and female antennae were analyzed by RT-PCR. Results for all of these genes are listed in Figure 7. The RT-PCR results showed all of the 37 LdecORs expressed in the antennae, but the expression level was quite low. For the control genes LdecRL31 and LdecRPS3, the 28 cycle of amplification was sufficient for detection. Conversely, for all the candidate LdecORs (including LdecOrco), the bands were difficult to detect unless the cycle-numbers increased to 38. One candidate OR- LdecOR6 was detected to expressed only in male antennae. Except LdecOR6, the expressions of all the other candidate ORs were detected in both male and female antennae. The expression of LdecOR5, LdecOR12, LdecOR26, and LdecOR32 was clearly higher in male compared to female, and LdecOR3 and LdecOR29 expressed higher in female.

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.