The insects H. parallela were obtained from the Institute of Plant Protection, Chinese Academy of Agricultural Sciences, Beijing, China. They were reared in plastic containers (100×50cm) with damp soil (20% moisture) at 25°C, 80% RH, and on a 12-hL:12-hD photoperiod and supplied with fresh elm, Ulmus parvifolia Jacquin, leaves until RNA extraction (Yi et al., 2018). Various tissues, including antennae, heads, thoraxes, abdomens, legs, and wings, were dissected separately and immediately thrown in liquid nitrogen. These tissue samples were frozen at −80°C until use. Total RNAs of all tissues were extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer’s protocol. The RNA integrity and purity were examined by 1.2% agarose electrophoresis and with a NanoDrop^™ spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States).

Fifty male or female antennae were used for the RNA extraction and transcriptome analysis. The poly (A) mRNA was separated from 20μg of total RNA and purified with Oligo d(T) magnetic beads, and fragmented into short fragments by fragmentation buffer. The first-strand of cDNA was synthesized using a random hexamer primer with these mRNA fragments as templates. Next, the second-strand of cDNA was synthesized by adding DNA polymerase I, dNTPs, and RNase H and purified with a QIAquick PCR purification kit (Qiagen, Hilden, Germany), resolved with EB (ethidium bromide) buffer for end reparation and single nucleotide A (adenine) addition to the 3' end of the cDNA. Next, the short fragments were connected with sequencing adapters. After that, fragment sizes were assessed by agarose gel electrophoresis, and the appropriate fragments were subjected to PCR amplification and sequencing (Illumina HiSeq^™ 2000, San Diego, CA, United States) by the Beijing Genomics Institute (BGI) sequencing company (Shenzhen, China). After sequencing, image deconvolution and quality value calculations were performed using the Illumina GA pipeline 1.3 (Huang et al., 2012). Then, the low quality reads (with >50% of nucleotides for which the Phred Quality Score Q was less than or equal to 5) were filtered out to generate clean reads. Transcriptome de novo assembly was carried out with short reads assembling program – Trinity (version 20130225;¹ Grabherr et al., 2011). Firstly, clean reads with a certain length of overlap were combined to form longer contiguous sequences (contigs). Then the clean reads were mapped back to contigs; with paired-end approaches it was able to detect contigs from the same transcript as well as the distances between these contigs. Next, Trinity connected the contigs, and sequences that could not be extended were obtained. These result sequences of Trinity were defined as ‘unigenes’ by the BGI Company (Shenzhen, China; Xiao et al., 2015; Xie et al., 2015).

The unigenes were firstly aligned to protein databases using BLASTx, including those from the database of Nr, COG,² Swiss-Prot,³ and the KEGG⁴ with a cut-off E-value of <10⁻⁵, and they were also aligned to the Nt database using BLASTn with a cut-off E-value of <10⁻⁵. With Nr annotation, Blast2GO (v2.5.0) was used to obtain gene ontology (GO) annotation of the unigenes. After getting GO annotation for every unigene, Web Gene Ontology Annotation Plot software was used to complete GO functional classification for all unigenes. With the KEGG database, the metabolic pathway annotation for unigenes was also performed. The above detailed steps are similar to a previous study (Huang et al., 2012).

Candidate CXE genes were chosen from the transcriptome data. Further, all candidate CXEs were manually checked by the BLASTx program at the National Center for Biotechnology Information (NCBI). The open reading frames (ORFs) of candidate CXE genes were predicted using the ORFfinder program.⁵ Then, the theoretical pI and Mw of deduced CXE proteins were calculated using the ExPASy tool with average resolution.⁶ Next, the putative N-terminal signal peptides of deduced CXE proteins were predicted using the SignalP 5.0 server, and the organism group was set to Eukarya.⁷ Multiple sequence alignment of identified CXE sequences was generated using the online Clustal Omega program.⁸

The amino acid sequences for constructing phylogenetic trees were obtained from Coleoptera, Lepidoptera, Hymenoptera, and Diptera species, including Tribolium castaneum (Tcas), Leptinotarsa decemlineata (Ldec), Bombyx mori (Bmor), A. polyphemus (Apol), S. littoralis (Slit), S. exigua (Sexi), Spodoptera litura (Slitu), Drosophila melanogaster (Dmel), Anopheles gambiae (Agam), and Apis mellifera (Amel). Their accession numbers from Genbank are listed in Supplementary Table 2. The amino acid sequences were aligned using ClustalX 2.0 software.⁹ The unrooted neighbor-joining (NJ) trees of candidate CXEs with full-length ORFs were constructed using the MEGA 7.0 software with the p-distance model.¹⁰ Gaps/missing data were treated as partial deletion with a site coverage cut-off=95%. Node support was assessed using a bootstrapping procedure based on 1,000 replicates. Some CXEs functionally characterized as ODEs or restrictively expressed in olfactory sensilla were marked with black dots in the phylogenetic trees from species A. polyphemus (Apol-PDE and Apol-ODE; Ishida and Leal, 2002, 2005), D. melanogaster (DmelEst6, DmelJHE and DmelJHEdup; Chertemps et al., 2012, Steiner et al., 2017, Hopkins et al., 2019), S. exigua (SexiCXE4, 10, 14; He et al., 2014a,b, 2015), S. littoralis (SlitCXE7, SlitCXE10 and Slit-EST; Durand et al., 2010a, 2011, Merlin et al., 2007), Sesamia nonagrioides (Snon-EST; Merlin et al., 2007), and P. japonica (Pjap-PDE; Ishida and Leal, 2008). The generated phylogenetic trees were colored and arranged using Figtree 1.42 software¹¹ (Yi et al., 2018). CXEs were mainly divided into secreted enzyme, intracellular enzyme, and neuro-signaling enzyme clades according to the classification system of CXEs described previously (Durand et al., 2010b).

Fifty male or female antennae and ten male or female heads, thoraxes, abdomens, legs, and wings were used as a biological sample for the RNA isolation and qPCR analysis. RT-qPCR was performed on a StepOne Plus Real-time PCR System (Applied Biosystems, Foster City, CA, United States) using SYBR Premix ExTaq II (Tli RNaseH Plus; Takara, Dalian, China). Candidate GAPDH, actin, and 18s rRNA were selected to evaluate the suitability as internal reference genes, and GAPDH had the most stable expression. So GAPDH was chosen as the reference gene for qPCR analysis. Primers of 20 CXEs and the reference gene (GAPDH) were designed using the Primer Premier 5 software, and are listed in Supplementary Table 3. The primer efficiencies of each gene were calculated and the primers with efficiency values ranging from 0.95 to 1.05 were selected for further experiments. RT-qPCR was performed under the following conditions: 30s at 95°C, 40cycles of 5s at 95°C, 10s at 55°C, and 34s at 72°C. This program was followed by a melting temperature analysis: 95°C for 15s, 60°C for 1min, increasing 0.3°C per min, and 95°C for 15s. Each reaction was run in triplicate from three biological replicates. The expression levels of these genes were calculated using the 2^-ΔΔCt method.

Significant differences were analyzed by t-test (between sexes) and one-way ANOVA (between different tissues), followed by Tukey’s HSD multiple comparisons test using SPSS software version 13.0. Expression profiles were created using the software Prism 6.0 (GraphPad Software, CA, United States; Yi et al., 2019). The significant differences between sexes were marked with asterisks (^*, p<0.05; ^**, p<0.01, NS, no differences). The expression levels among different tissues of each gender followed by the different lowercase or uppercase letters were significantly different.

Using the Illumina HiSeq2000 sequencing platform, 34,706 unigenes in total were obtained in the antennal transcriptomes and analyzed by searching against Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide database), Swiss-Prot, KEGG, COG (Cluster of Orthologous Groups of proteins), and GO databases. As to the results, significant matches were found against 18,312 (52.76%) unigenes in the Nr database, 8,524 (24.56%) unigenes in the Nt database, 14,071 (40.54%) unigenes in the Swiss-Prot database, 6,879 (19.82%) unigenes in the KEGG database, 6,508 (18.75%) unigenes in the COG database, and 8,919 (25.70%) unigenes in the GO database. A total of 19,025 (54.82%) unigenes were successfully annotated in at least one database, while 15,681 (45.18%) unigenes had no matching sequences in any of these databases (Table 1). The sequencing data were available at the NCBI Sequence Read Archive database with accession number SRP233063.

In order to explore the role of CXEs in the process of odorant signal inactivation, it is necessary to identify candidate CXEs related to odorant degradation by blasting against the Nr database. A total of 20 candidate CXEs were identified with complete ORFs. The length of deduced HparCXEs was between 535 and 854 amino acid residues (Table 2). The predicted theoretical isoelectric point (pI) of HparCXEs was from 4.29 to 9.09 with the predicted molecular weight (Mw) of 60.6kDa to 98.4kDa. Furthermore, all HparCXEs except for HparCXE1, 11, and 16 were predicted to have the putative N-terminal signal peptide (SP), indicating that they could be secretory proteins. Multiple sequence alignments of candidate HparCXEs showed that they had the conserved pentapeptide motif of insect CXEs, Gly-X-Ser-X-Gly (“X” represents any residue) of typical CXE proteins except for HparCXE5, 9, and 20 (Figure 1). In addition, conserved oxyanion hole-forming residues (Gly, Gly, and Ala) were also found, which were thought to stabilize the transition states of the hydrolysis reaction (Durand et al., 2010b). The oxyanion hole-forming residues of four HparCXEs (HparCXE2, 3, 7, and 15) and nine HparCXEs (HparCXE1, 4, 5, 8, 11, 14, 16, 17, and 18) consisted of “G-G-A” and “G-G-G,” respectively (Figure 1). Additionally, the HparCXEs possessed the three conserved catalytic residues: serine (S), glutamate (E), and histidine (H), except HparCXE9 and 20 (lacking S), HparCXE7 (lacking E) and 19 (lacking E and H).

In order to identify the phylogenetic relationships between CXEs of H. parallela and those of other insect species, the NJ tree of H. parallela CXEs was constructed together with CXEs from several Coleoptera, Lepidoptera, and Diptera species (Figure 2). The phylogenetic analysis grouped all 20 HparCXEs into three families: the secreted enzyme family (clade II, III, V-VII: 8 HparCXEs), the intracellular enzyme family (clade I, IV: 10 HparCXEs), and the catalytically inactive, neuro-signaling family (VIII-XI: 2 HparCXEs). HparCXE1, 7, 8, 11, and 12 were clustered into the α-esterase clade in the intracellular enzyme family, well known for their involvement in detoxification of insecticide/xenobiotics and digestion of food esters (Durand et al., 2010b). This clade included some functionally characterized ODEs, such as SlitCXE10 from S. littoralis (Durand et al., 2010a) and SexiCXE10 from S. exigua (He et al., 2015). HparCXE13 shared a close relationship with juvenile hormone esterases (JHEs) of other species in the secreted enzyme family clade, with more than 90% bootstrapping support with JHEs of other species. In this clade, DmelJHE and DmelJHEdup have been functionally characterized and have high degrading activities against ester plant volatiles (Steiner et al., 2017; Hopkins et al., 2019). HparCXE6, 9, and 10 were clustered into the cuticular/antennal esterase clade in the secreted enzyme family. HparCXE2, 4, 15, and 17 were clustered into the β and pheromone esterase clade in the secreted enzyme family, together with three functionally characterized ODEs, Apol-PDE (Ishida and Leal, 2005), Pjap-PDE (Ishida and Leal, 2008), and DmelEst6 (Chertemps et al., 2012). HparCXE3, 5, 14, 16, and 18 belonged to coleopteran xenobiotic metabolizing enzymes in the intracellular enzyme family. Only HparCXE19 and HparCXE20 were clustered into the neuro-signaling enzyme family (Figure 2).

Tissue expression profiles of all 20 CXEs were determined by RT-qPCR analysis (Figure 3). Five putative HparCXEs (HparCXE6, 10, 11, 13, and 16) were significantly expressed higher in the male antennae than in the female antennae, whereas only one HparCXE (HparCXE18) showed significantly higher expression in the female antennae than in the male antennae. HparCXE4, 11, 16, 17, 18, 19, and 20 were expressed higher in both male and female antennae than in other tissues, whereas other HparCXEs were widely distributed in heads, thoraxes, abdomens, legs, and wings. HparCXE5, 6, 14, and 15 were expressed at a higher level in legs than in other tissues in both sexes. HparCXE1 and HparCXE12 had higher expression levels in both female and male wings and legs than in other tissues. HparCXE9 was significantly expressed in both male and female heads.