Adult C. dieckmanni were collected in a hazelnut orchard in the Sandao Forestry Farm, Mudanjiang City, Heilongjiang Province, China (129°39′02″E, 44°07′04″N). The insects were fed with hydroponic C. heterophylla shoots in a cage (80 cm × 40 cm × 40 cm). Males and females were separated and dissected under a dissecting microscope to collect the antennae, proboscis, head, legs, wings, prothorax, and abdomen. Tissue samples were separately stored in Axygen tubes and frozen at −80°C. A total of 200 pairs of antennae were dissected from female and male adults for transcriptome sequencing, and there were three biological replicates for female and males, including six samples (DM1, DM2, DM3, DF1, DF2, and DF3).

Total RNA was isolated from the antennae of female and male adults of C. dieckmanni, using a Trizol reagent kit (Invitrogen, Carlsbad, CA, United States) based on the manufacturer’s instructions. The integrity of RNA quality was assessed using an Agilent 2100 bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and examined by RNase free agarose gel electrophoresis.

After the total RNA examination was qualified, the mRNA was enriched with magnetic beads with Oligo (dT). We added fragmentation buffer to break the mRNA into short fragments. We then used mRNA as a template for amplification and added random hexamers for reverse transcription to synthesize the first strand cDNA. We added buffer, dNTPs, and reagents, such as DNA polymerase I and RNase H, to synthesize the second-strand cDNA and purified the double-stranded cDNA with AMPure XP beads. The purified double-stranded cDNA was end repaired, poly (A) was added, and the cDNA was connected to Illumina sequencing adapters. We used AMPure XP beads to select the size of the fragments. PCR amplification was performed, and the PCR products were purified with AMPure XP beads to obtain the final cDNA library. After the library was constructed, Qubit 2.0 was used to determine the concentration, and the library was diluted to 1.5 ng/μl. Agilent 2100 was used to detect the insert size of the library. To ensure the quality of the library, the qPCR method was used to accurately quantify the effective concentration of the library. Qualified libraries of female and male adult antennae were sequenced using Illumina HiSeqTM 4000 by Gene Denovo Biotechnology Co. (Guangzhou, China). All the raw sequence data were saved in the NCBI Sequence Read Archive (SRA) database under BioProject accession No. PRJNA769532.

To ensure analysis quality, the original data obtained by sequencing were converted into sequence data by base calling, which contains joints, duplicates, and low sequencing quality. The reads were filtered to avoid affecting sequence assembly and subsequent analysis. We removed reads containing adaptors, reads with a proportion of unknown nucleotides (N) greater than 10%, and reads with low quality (q-value ≤ 20) over 40% bases. After obtaining the clean reads, we calculated the Q20, Q30 and GC content of the clean reads, and used Trinity to assemble the transcriptome from scratch. Trinity is a short reads assembly software (Grabherr et al., 2011).

The results of Trinity assembly were unigenes. To obtain the functional information of the unigenes, the unigene sequences were subjected to basic annotations, including protein function annotations, COG/KOG function annotations, pathway annotations and gene ontology (GO) annotations. For unigene annotation, at the E-value threshold of 1e-5, through BLASTx, all unigenes were compared in the NCBI non-redundant protein database (Nr), manually annotated, and reviewed in the Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, and the Clusters of Orthologous Groups of proteins (COG/KOG) database. According to the best alignment results, unigene sequence directions were determined and protein function annotations were obtained. We used Blast2GO software for GO annotation and unigenes Nr result annotation (Conesa, et al., 2005). We selected the top 20 unigenes with the highest scores and no less than 33 High-scoring Segment Pair (HSPs) hits for Blast2GO analysis. WEGO software was used to classify the unigenes (Ye, et al., 2006).

The predicted conserved domains and open reading frames (ORFs) were manually validated for all candidate chemosensory genes (OBPs, CSPs, ORs, SNMPs, GRs, and IRs). Sequences with errors (insertions/gaps/deletions in homopolymer regions) were deleted. The putative N-terminal signal peptides of ligand-binding proteins (OBPs and CSPs) were searched using the SignalP v4.1 program. The transmembrane domains (TMDs) of ORs, IRs, and GRs were searched using TMHMM Server version 2.0. The translation of nucleotide sequences was performed using the Expasy program.

To validate the annotation of candidate chemosensory genes and the identification of homologous sequences, we conducted a phylogenetic analysis of C. dieckmanni including closely related Coleoptera species and other model insects. Insects selected for analysis were those whose chemosensory genes have been well studied (Supplementary Table S6). A multiple sequence alignment was performed using the MAFFT program, and the results were accessed using the Jalview software for editing and selecting the color scheme. Amino acid sequences of chemosensory genes of selected insects were obtained from the NCBI database, and Maximum likelihood (ML) phylogenetic trees were built based on the Jones-Taylor-Thornton (JTT) model in MEGA-X. Initial trees for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the JTT model, then selecting the topology with the superior log likelihood value. The tree was drawn to scale with branch lengths measured in the number of substitutions per site (Jones et al., 1992; Kumar et al., 2018). Bootstrap support values were based on 100 replicates.

To identify chemosensory DEGs, we first used RPKM to estimate the expression of chemosensory genes in the antennae of male and female adults (Mortazavi et al., 2008). These genes were selected as reference for the expression profile of genes in the antennae of male and female adults. Gene expression levels were determined, according to RPKM values. In addition, DEGs were identified by Benjamini-Hochberg-corrected p values (Guo and Rao, 2008). As a key indicator, the false discovery rate (FDR) was used to screen in multiple hypothesis testing for gene selection. In comparison, genes with multiple changes ≥ 2 and FDR < 0.05 were identified as significant DEGs (Ran et al., 2020).

The 31 genes in female and male transcriptomes were validated by qRT-PCR in different tissues and these genes, compared with all candidate genes, were highly expressed based on RPKM. Antennae (480) were collected from female and male adults. RNA was extracted using Trizol reagent and cDNA was synthesized using a PrimeScriptTM RT Reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan). We designed gene-specific primers using Primer3. The amplification efficiency was calculated for the genes and two candidate internal reference genes (actin and EF1-a). Based on the results, actin and EF1-a were selected as the reference genes for qRT-PCR. PCR was performed using Mx3000P QPCR (Agilent, United States) and ChamQ universal SYBR qPCR Master Mix (TaKaRa). The PCR program was as follows: 95°C for 30 s; 40 cycles of 95°C for 5 s and 60°C for 30 s; and 65–95°C at 0.5°C/5 s. A melting curve was generated. Each reaction had three independent biological replicates, and the relative gene expression level was obtained by the 2^−∆∆Ct method (Guo and Rao, 2008). The relative gene expression in antennae was analyzed using SPSS v19.0 software. An independent sample t-test was performed to compare the gene expression levels between male and female antennae. Prism v6.0 was used to generate bar graphs.

A total of 42.5 Gb raw data were obtained after filtering. The transcriptomes of three female individuals yielded 49,095,002 (97.56%, DF1), 44,906,190 (97.60%, DF2), and 49,073,722 (97.47%, DF3) clean reads. The transcriptomes of three male individuals yielded 55,983,482 (97.52%, DM1), 49,406,426 (97.53%, DM2), and 41,989,468 (97.44%, DM3) clean reads. The Q20 value was >98%, and the Q30 value was >95% (Supplementary Table S1.1). The clean reads of six samples of females and males were assembled to yield a total of 43,060,197 unigenes. There were 44,024 unigenes longer than 200 bp (N50 = 1,861), including 42,971 unigenes from female samples (DF) and 43,190 unigenes from male samples (DM). The maximum length was 27,588 bp, the minimum length was 201 bp, and the average length was 978 bp (Supplementary Figure S1, Supplementary Table S1.2).

The unigene sequences were aligned using BLASTx against the protein databases NR, SwissProt, KEGG, and KOG (e-value < 0.00001). The protein with the highest similarity to a given unigene was selected for protein function annotation, and 21,803 of 44,024 (49.52%) unigenes were successfully annotated. The numbers of reference gene sequences annotated in NR, Swiss-Prot, KEGG, and KOG were 21,733, 12,094, 7,502, and 11,403, respectively, (Supplementary Figure S2, Supplementary Table S1.2). The sequence with the best alignment (lowest e-value) of a unigene, which was selected as homologous in the NR database (if there was a tie, the first sequence was selected). The species of the homologous sequence was determined, and the number of homologous sequences of each species was counted. The sequences with high similarity to those in transcripts were searched in the NR protein database by comparison with BLASTx. The analysis revealed a large number of homologous genes in Coleoptera. The insect species with the highest homology (40.93%) was Dendroctnous ponderosae, followed by Anoplophora glabripennis (14.13%), Tribolium castaneum (5.01%), and Lasius niger (4.68%), which had 8,895, 3,071, 1,088, and 1,018 matched transcripts, respectively (Supplementary Figure S3).

Based on the NR annotation, the Blast2GO (Conesa et al., 2005) software was used to obtain the GO annotations of unigenes. The WEGO software (Ye et al., 2006) was used for GO functional classification of all unigenes to draw a whole picture of the gene functions of this species. These transcripts were assigned to three categories: molecular function (1,808 unigenes), biological process (3,675 unigenes), and cellular component (2,544 unigenes; Supplementary Figure S4). To accurately evaluate the annotation information, the transcripts were aligned to the KOG database, and a total of 11,403 transcripts were assigned to 25 functional categories (Supplementary Figure S5).

Of the chemosensory gene families, a total of 87 expressed chemosensory proteins were found in C. dieckmanni: 23 odorant-binding proteins, 15 chemosensory proteins, 2 sensory neuron membrane proteins, 15 odorant receptors, 23 ionotropic receptors, and nine gustatory receptors (Information for genes in Supplementary Tables S2.1–S2.6). The sequences of chemosensory genes identified in C. dieckmanni transcriptomes are shown in Supplementary Table S5.

The analysis of transcriptome data revealed 15 CSPs, which had four conserved cysteine residues and the featured sequence C1-X6-7-C2-X18-C3-X2-C4 (X is any amino acid) of typical insect CSPs (Figure 3). The genes of 15 CSPs had ORFs ranging from 231 to 963 bp and encoding 76–320 amino acids, of which the genes of seven CSPs (CdieCSP1, CdieCSP2, CdieCSP4, CdieCSP5, CdieCSP6, CdieCSP7, and CdieCSP8) had incomplete ORFs, and eight CSPs (CdieCSP1, CdieCSP2, CdieCSP4, CdieCSP5, CdieCSP8, CdieCSP10, CdieCSP12, and CdieCSP13) lacked signal peptides (Supplementary Table S2.2). A phylogenetic tree was built using the 15 CSPs obtained from the transcriptomes and CSPs from five Curculionoidea species (Dendroctonus adjunctus, D. armandi, Lissorhoptrus oryzophilus, Dendroctonus ponderosae, and Tomicus yunnanensis) (Figure 4).

Genes with complete ORFs encoding CdieSNMP1 (SNMP1) and CdieSNMP2 (SNMP2) were identified based on BLASTx and cluster analysis. The ORF was 1,695 bp for CdieSNMP1 and 1,575 bp for CdieSNMP2, indicating that both were encoded by full-length genes. CdieSNMP1 and CdieSNMP1 had six and three transmembrane regions as well as 74.04 and 56.84% homology, respectively, (Supplementary Table S2.3). Phylogenetic analysis showed that CdieSNMP1 and CdieSNMP2 clustered with the SNMP1 and SNMP2 of other Curculionoidea species, respectively, and were most closely related to D. ponderosae, which confirmed the above results (Figure 5).

A total of 15 ORs were identified in the analysis of transcriptome data, including one odorant receptor coreceptor (Orco). The alignment of ORs against the NR database showed 36.02–94.40% homology and 78–1,449 bp ORFs encoding 25–482 amino acids. Given the sequence characteristics of ORs, 1–9 transmembrane domains were predicted (Supplementary Table S2.4). A phylogenetic tree was established along with the ORs identified in other Curculionoidea species (D. ponderosae, Ips typographus, and Sitophilus oryzae; Figure 6). Based on the amino acid sequence analysis, the ORs from C. dieckmanni was contained in six groups of beetles (Bin et al., 2017; Mitchell et al., 2020). Among them, CdieOR15 clustered with the highly conserved ORco in a branch. CdieOR14 clustered in group 1, CdieOR12 clustered in group 3, CdieOR10 clustered in group 5, three candidate ORs (CdieOR2, CdieOR7, and CdieOR9) clustered in group 7, and four candidate ORs (CdieOR5, CdieOR6, CdieOR8, and CdieOR11) clustered in group 2.

We identified 23 IRs from C. dieckmanni antennae transcriptomes. A phylogenetic tree was built along with the IRs identified in other Curculionoidea species (Cylas formicarius, D. ponderosae, Ips typographus, L. oryzophilus, and S. oryzae; Figure 7). Phylogenetic analysis showed that CdieIR2, CdieIR4, CdieIR15, and CdieIR17 were clustered in IR75; CdieIR10 clustered in IR8a; CdieIR11 clustered in IR25a; CdieIR13 clustered in IR76b; and CdieIR21 clustered in IR41a. These were highly conserved. Related genes were identified based on other Curculionoidea species. The genes encoding nine IRs (CdieIR1, CdieIR8, CdieIR11, CdieIR12, CdieIR13, CdieIR19, CdieIR20, CdieIR22, and CdieIR23) had ORFs similar to or longer than 1,600 bp, and the ORF of CdieIR7 gene was 1,431 bp. The lengths of these nine IRs were comparable to the full lengths of such genes, and the genes of the remaining 14 IRs (CdieIR2, CdieIR3, CdieIR4, CdieIR5, CdieIR6, CdieIR7, CdieIR9, CdieIR10, CdieIR14, CdieIR15, CdieIR16, CdieIR17, CdieIR18, and CdieIR21) had partial ORFs of such genes. The number of transmembrane domains of the 23 IRs ranged from one to six (Supplementary Table S2.5).

From the transcriptome data, we identified 9 GR candidate transcripts with 1–11 transmembrane domains and homology ranging from 29.16 to 61.07%. The genes of 5 GRs (CdieGR1, CdieGR5, CdieGR6, CdieGR8, and CdieGR9) had ORFs longer than 1,000 bp (Supplementary Table S2.6). A phylogenetic tree was constructed to analyze the functional classification of GRs in C. dieckmanni transcriptomes based on the GRs identified in other Curculionoidea species (C. formicarius, D. ponderosae, I. typographus, L. oryzophilus, and S. oryzae; Figure 8). Phylogenetic analysis showed that three GRs (CdieGR4, CdieGR7, and CdieGR8) clustered in the gustatory receptor for sugar receptors, while CdieGR1 clustered in the gustatory receptor for CO₂ receptors.

The RPKM analysis showed that a total of 355 genes were differentially expressed in male and female antennae, among which, 69 were significantly male-biased and 286 were significantly female-biased (Figure 9, Supplementary Tables S3, S4), including only three chemosensory genes. Two candidate odorant-binding proteins (CdieOBP19 and CdieOBP20) were expressed at significantly higher levels in female antennae than in male antennae (FDR < 0.05) (Supplementary Table S4). One candidate odorant receptor protein (CdieOR13) was expressed at significantly higher levels in female antennae than in male antennae (FDR < 0.05) (Supplementary Table S4).

To validate and further understand the functions of the chemosensory genes in adult female and male C. dieckmanni antennae, we studied the expression profiles of these genes using qRT-PCR. The primers are listed in Supplementary Table S7. The qRT-PCR results showed that all chemosensory genes were expressed in both male and female antennae, except CdieOBP19, which was only expressed in female antennae. Three OBPs (CdieOBP9, CdieOBP20, and CdieOBP21) were significantly female-biased, one OBP (CdieOBP8) was significantly male-biased, one CSP (CdieCSP15) was significantly female-biased, two ORs (CdieOR13 and CdieOR15) were significantly female-biased, and no SNMPs, IRs, or GRs showed differential expression between males and females (Figure 10).