A total of three species used in this study, N. flavidorsalis, C. castaneipars, and S. postornata, were collected from the coniferous and broad-leaved mixed forest by light trapping in Mt Tiantong, Zhejiang Province, China (29°49′N, 121°48′E) (Supplementary Figure 1). A number of five male and five female moths were segregated into different cages (40 cm × 40 cm × 60 cm) containing 10% (w/v) sugar solution until the use for microscopy. Antennae from other males and females were excised and stored in RNAlater (Ambion, Austin, TX, United States). Then, all samples were taken back indoors and held at –80°C until RNA isolation. The specimens were identified by Prof. Chun-sheng Wu (Institute of Zoology, Chinese Academy of Sciences, China), and DNA barcoding was also used in species identification. Their voucher specimens are preserved at the Entomological Museum of Capital Normal University, China, and the number was Lep/Lim/191027 to Lep_Lim_191119.

The adults of the slug moth were first anesthetized by freezing at −20°C for 5 min. Antennae were then dissected and immersed in a freshly prepared fixative solution containing 2% paraformaldehyde and 2.5% glutaraldehyde mixed with a phosphate-buffered solution (PBS) (pH 7.4) for 24 h at 4°C. The antennae were kept at 10% KOH at 80°C for 30 min to remove the scales. Subsequently, the antennae were dehydrated using a graded ethanol series (30, 50, 70, 80, 90, 95, and 100%) followed by the critical-point drying (Leica EM CPD 030, Tokyo, Japan). The dried specimens were carefully glued onto SEM stubs and sputter-coated with gold (Eiko IB-5 Ion Coater, Tokyo, Japan; 45 s, 20 mA). The preparations were examined using a Hitachi SU-8010 cold field emission scanning electron microscope (Japan) under 3 kV voltage.

Given that different terminological terms based on the morphological characters have been used on the antennal sensilla of Lepidoptera, this study primarily followed the study of Schneider (1964) and Jefferson et al. (1970). The antennal sensilla of N. flavidorsalis, C. castaneipars, and S. postornata were identified, counted, and measured according to the previous measurement method (Onagbola and Fadamiro, 2008; Ivanov and Melnitsky, 2016). Means were based on the measurements (μm) from at least 20 photomicrographs of individual sensilla of the identical type. The types and abundance of sensilla were compared between species and genders and analyzed with the two-way ANOVA with R. The level of significance in all tests was set at 0.05.

Narosoideus flavidorsalis and S. postornata were selected for antennal transcriptome analysis. Total RNA was extracted from twelve adults’ antennae from each species using TRIzol reagent (Invitrogen, Carlsbad, CA, United States) and the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. RNA quantity and purity were checked using the NanoDrop 8000 (Thermo, Waltham, MA, United States). A total amount of 3 μg RNA per sample was used for cDNA library construction. All samples had RIN values above 8. The RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). The cDNA library construction and subsequent Illumina sequencing of samples were performed at Novogene Corporation (Beijing, China). The cDNA libraries were generated using the TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, United States). Random hexamer primers were used to synthesize the first-strand cDNA, then synthesizing the second-strand cDNA using buffer, dNTPs, RNase H, and DNA polymerase I at 16°C for 1 h. The remaining overhangs were converted into blunt ends via exonuclease or polymerase activities and removed the enzymes. After end repair, A-tailing, and the ligation of adaptors, the products were amplified by PCR and quantified precisely using the Qubit DNA Br Assay Kit (Q10211; Invitrogen, Carlsbad, CA, United States). The library fragments were purified using the MinElute Gel Extraction Kit (Qiagen, Hilden, Germany). The resulting cDNA library preparations were then sequenced on an Illumina HiSeq-2000 platform. The quality of the sequences generated from the PE 200 bp, and all mate-pair libraries were assessed using FastQC (Brown et al., 2017).

To ensure the accuracy of sequence assembly, the clean data were obtained from raw reads through the following steps: filtered out the reads with adapters, deleted the reads with uncertain bases more than 10%, then removed low-quality (the bases with sequencing error rates greater than 1% are more than 50% in the read) and adaptor sequences by Trimmomatic¹. At the same time, the Q20, Q30, GC-content, and other related information of the clean data were calculated. All downstream analyses were based on the clean data with high-quality reads. Transcriptome assembly was carried out with program trinity, in which all the parameters were set to their defaults (Grabherr et al., 2011). The transcripts longer than 200 bp were performed by NCBI BLASTx searches to databases (Nr, Pfam, Swissprot, KOG, and KEGG); the given transcripts were functionally annotated as the retrieving proteins or nucleic acid with the highest sequence similarity. The blast results were then imported into the Blast2GO pipeline² for GO annotation (Conesa et al., 2005; Götz et al., 2008).

All candidate chemosensory genes (OBPs, CSPs, SNMPs, ORs, GRs, and IRs) in N. flavidorsalis and S. postornata and their open reading frames (ORFs) were manually verified by BLASTx searches against custom-made databases and the non-redundant nucleotide collection at NCBI. The ORFs were identified, and the annotation was confirmed by the additional BLAST searches³. Transmembrane domains of ORs, IRs, and GRs were predicted using the default parameter of TMHMM2.0 and TMPred, and the N-terminal signal peptide of the candidates, OBPs and CSPs, was predicted by SignalP4.0 (Quevillon et al., 2005).

For further verification of the candidate chemosensory genes and identification of orthologs, phylogenetic analyses were conducted among two slug moths and other related Lepidoptera species, such as Bombyx mori, Helicoverpa armigera, and Holcocerus hippophaecolus. The available amino acid sequences of chemosensory genes identified in these species were downloaded manually. Maximum likelihood trees were reconstructed using the predicted OR, IR, OBP, CSP, and SNMP protein sequences and orthologs in other species of Lepidoptera and model insects (Drosophila melanogaster) to analyze the characteristics of olfactory genes in two species of family Limacodidae.

Amino acid sequences were aligned with MAFFT (version 7.308) (Katoh and Standley, 2013), and the corresponding maximum-likelihood trees were constructed in IQ-TREE (version 2.1.7) using best-fitting substitution model (GTR + I + G) (Trifinopoulos et al., 2016). The tree structure and node support reliability were evaluated by the bootstrap analysis with 1,000 replicates. The phylogenetic trees were colored and arranged in FigTree (version 1.4.2)⁴.

Ecological niche modeling (ENM), a widely accepted method to visualize the distribution patterns (Peterson and Soberón, 2012), was adopted to compare the area of suitable optimal habitats of N. flavidorsalis and S. postornata. Meanwhile, the habitat suitability of the host plants of these two species was also predicted. Pyrus sorotina was selected as the dominant host plant of N. flavidorsalis. A number of three plants, Juglans regia, Populus simonii, and Camellia sinensis, were chosen as representative host plants of euryphagous S. postornata. The occurrence records for niche modeling analysis were collected from the online dataset Global Biodiversity Information Facility (GBIF) and published documents. The present bioclimatic variables were downloaded from the WordClim website⁵. The maximum entropy (MaxEnt) approach was employed to predict the species distribution using the ENMTools version 1.0.4 package in the R platform with MaxEnt.jar version 3.4.4 (Warren and Seifert, 2011). A combination of six MaxEnt feature classes [linear (L), linear quadratic (LQ), hinge (H), linear quadratic hinge (LQH), linear quadratic hinge product (LQHP), and linear quadratic hinge product threshold (LQHPT)] was tested to optimize model parameters for the calibration. We trained that the MaxEnt models were trained with the sub-sample method, in which 1/3 of the presence points were set aside for model testing evaluation. The model was trained by replicating each model ten times, and the final output used was the average of all runs. We selected the model with a combination of feature classes and the regularization multiplier that had the lowest akaike information criterion (AICc) value, which could describe the model fit and complexity (Shin et al., 2021).

The antennae structure of the three species is broadly in line with other lepidopteran insects, which consisted of a scape, a pedicel, and a flagellum. Based on the antennal data from male samples, there was no significant difference in the pedicel and total antennae length among the three species, but only in the side-branches size (Table 1). There was no significant difference in antennae length among the three female species, similar to the male data. Still, a significant difference was detected in the length and width of flagellomeres among the three species (Table 2). The antennae of all three species bear a pair of sexual dimorphic antennae, the female ones were filiform, and the male ones were bipectinate (Supplementary Figure 2). The number of flagellomeres of these three species was different, which was consistent with the difference in antennae length. Meanwhile, the number of flagellomeres of the males was slightly more than that of the females (Supplementary Table 1).

A total of nine types of sensilla were identified on the antennae of three slug moths, in which five types each having two types of sensilla: sensilla chaetica (SCh I, SCh II), Böhm’s bristles (BB I, BB II), sensilla styloconica (SSt I, SSt II), sensilla coeloconica (SCo I, SCo II), and sensilla furcatea (SFu I, SFu II). In these types of sensilla, SB, SSq, and SCo I were detected only in N. flavidorsalis, while SCh II and SFu II were detected only in S. postornata. In general, N. flavidorsalis has more types of sensilla than the other two species (Table 3).

Sensilla Trichoderma (STr) was the most widespread sensilla in the three species, densely distributed on the ventral side of the flagellum spindle and lateral branches of both male and female antennae. The surface facial structure of STr among the three species was different: N. flavidorsalis was irregular, and the others were arranged in an oblique line (Supplementary Figure 3). Sensilla basiconica (SB) is located vertically in a shallow pit with a wrinkled proximal base (Supplementary Figure 4). SCh mainly found on the end of the antennal axis and lateral branches. Two subtypes of SCh (I and II) were identified based on the pattern of these grooves: SCh I were longitudinally striated from base to the end, while in SCh II, the grooves of the stripe gradually changed from longitudinal line to the imbricated texture (Figure 1). BB sensilla, which were found in clusters at the bases of the scape and pedicel, were quite long with a smooth surface and sharp end. The difference between the BB I and BB II was whether the end of the bristles was bifurcated or not (Figure 1). Sensilla squamiformia (SSq) was similar to the antennae scale structure, with no porous structure and conspicuous longitudinal stripes (Supplementary Figure 5). SSt with a terminal conical protrusion (without a pore) were on the surface of a decorative pattern of a cylindrical protrusion from the antennal surface. The two subtypes were found based on their structure: SSt I has a long columnar and a small cone, and SSt II has a short columnar with a longer conical body (Supplementary Figure 6). SCo were irregularly scattered on the surface of the flagellum, which consists of multiple longitudinal grooves and no stoma distribution. It can be divided into two subtypes according to whether the base is in the circular cavity of the epidermal bulge: SCo I was covered by the socket at the bottom, and SCo II was wholly exposed at the base (Figure 2). The sensilla uniporous peg (SUp) base was stuck into the convex epidermal with a smooth surface; the top was blunt and perforated (Figure 1). The SFu gradually narrowed from the base to the end and bifurcated at the end. There were two subtypes in S. postornata. The SFu II generally had longer poles and shorter bifurcation than SFu I (Supplementary Figure 7).

There were significant differences in several sensilla types between female S. postornata and female N. flavidorsalis, which includes sensilla trichoidea, sensilla coeloclnica (SCo), SSt, BB. The differences between species C. castaneipars and the other two species were not noticeable. There was no significant difference of SCh I among species, but there was a considerable difference in sensilla length among species and sex; the sensilla from females is significantly longer than males (Table 4). The morphological differences of sensilla among species were mainly reflected in females. There were significant differences in the length of sensilla BB I between S. postornata and the other two species, which was significantly longer than others in females (Table 4). There were significant differences in the SstI between the female S. postornata and the other female samples, such as the length, basal width, and the size of conical extremist of sensilla (Table 4).

The transcriptomic sequence data were generated using the antenna cDNA library and illumina HiSeq2500 sequencing platform. A total of 76,763,082 and 101,749,890 raw reads were obtained from N. flavidorsalis and S. postornata, respectively. After removing adaptor sequences, low-quality sequences, and N-containing sequences, approximately 73.3 million and 99.0 million clean reads were retained. The assembly of all clean reads together led to the generation of 148,264 and 198,318 unigenes (Table 5). The clean reads for N. flavidorsalis and S. postornata. have been deposited in the NCBI SRA database under the accession numbers SRR15330236 and SRR15330240.

We used the unigenes assembled in the transcriptome as the queries in BLASTx searches of the GO database, and 11,543 unigenes of N. flavidorsalis and 15,902 unigenes of S. postornata were annotated. All of the unigenes were divided into three categories: molecular function, cellular component, or biological process according to the biological processes and functional annotations. GO annotation indicated that the antennal transcriptomes of these two species were highly similar concerning GO terms (Supplementary Figure 8). In the biological process terms, single-organism, cellular, and metabolic occupied the majority of both differentially expressed unigenes. Cell and cell parts were the most abundant for all unigenes in the cellular component terms. In the molecular function category, catalytic activity, binding, and transporter activity had huge preponderance; however, unigenes in “signal transducer activity” and “chemoattractant activity” were also present (Supplementary Figure 8).

There were more putative OBPs identified in S. postornata than in N. flavidorsalis; however, the numbers of putative CSPs and SNMPs have no significant difference between these two species. A number of fourteen and eighteen OBPs, including common odor-binding protein (GOBPs) and sex pheromone-binding protein (PBPs), were identified in the antennae transcriptome of N. flavidorsalis and S. postornata, separately. The sequence identities of the OBPs with other lepidopteran insects ranged from 40 to 92% in the NCBI database (Supplementary Table 2). A phylogenetic tree of the OBPs was constructed based on the orthologs from Drosophila melanogaster (Dmel) and the three lepidopteran species, which are Bombyx mori (Bmor), Helicoverpa armigera (Harm), and Lobesia botrana (Lbot) (Figure 3A). The phylogenetic analysis demonstrated that the lepidopteran PBP and GOBP sequences were highly conserved and clustered into three lineage-specific clades according to their different functions. Meanwhile, other OBPs showed an extremely divergent trend. There were two GOBPs (Spos_OBP15 and Spos_OBP16) and one PBP (Spos_OBP18) identified in S. postornata, while only one GOBP (Nfla_OBP10) was detected in N. flavidorsalis. In each Limacodidae species, 14 unigenes were identified as CSP genes (Supplementary Table 3). All CSPs shared high sequence identities to known lepidopteran CSPs with an average of 67% identity. Of the 14 unigenes in N. flavidorsalis, 12 contained full-length ORFs encoding 107–152 amino acid residues. Unlike N. flavidorsalis, 9 out of 14 unigenes containing full-length ORFS were detected in S. postornata. A phylogenetic tree of the CSP was constructed based on the orthologs from D. melanogaster (Dmel), Helicoverpa armigera (Harm), B. mori (Bmor), and Eogystia hippophaecolus (Ehip). The amino acid identities between the orthologous CSPs in the two moths were relatively high, and most of the CSPs appeared in pairs on the dendrogram (Figure 3B). Only four SNMPs (3 with complete ORFs) were identified in N. flavidorsalis. In addition, five SNMPs were detected in S. postornata, but only one with complete ORF (Supplementary Table 4). As expected, SNMPs were clustered into two branches with other SNMP1 and SNMP2 orthologs from other lepidopteron (Supplementary Figure 9).

Unlike the results of OBPs, more ORs and IRs were detected in N. flavidorsalis than in S. postornata. We identified transcripts encoding 76 putative ORs in N. flavidorsalis, among which 33 likely represented full-length genes. In S. postornata, we identified 61 candidate OR genes comprising 32 full-length genes (Supplementary Table 5). Sex pheromone receptors (PRs) and co-receptor (ORco) were marked with orange and blue lines in the phylogenetic tree (Figure 2), in which both contain members from Lepidoptera and Drosophila (Dmel). This result shows that 15 putative PR genes of S. postornata were clustered into the PR subfamily and one into the ORco subfamily. Meanwhile, 27 putative PR genes and 1 putative ORco genes of N. flavidorsalis were assigned to PR and ORco subfamilies, respectively. The phylogenetic analysis showed a separate branch of Lepidoptera in OR, and the genes of PR and ORco subfamilies were relatively conservative (Figure 2A).

Ionotropic gene family prediction results showed that there were 27 members in the IR family of N. flavidorsalis, among which 17 genes were full length (the length of ORF ranged from 137 to 1,057 aa), while there were 22 members of S. postornata, in which there were 10 with full length (the size of ORF ranged from 107 to 1,059 aa) (Supplementary Table 6). In the phylogenetic tree reconstructed with IR genes, the 16 IR gene member clusters were marked with different colors. The putative IRs in N. flavidorsalis were clustered into 11 known branches, while the putative IRs in S. postornata were classified into 13 branches. Meanwhile, other IR members from these two species were clustered together, which indicated that the IR family was highly conservative (Figure 2B).

Family branches were classified and labeled according to the functions identified by the previous studies (Jiang et al., 2015; Hu et al., 2016; Xu et al., 2016). For example, bMOR_GR7-9, BMOR_66-67 were sugar receptors, Harm_GR9 and 13 were GR_43a receptor (fructose receptor), Dmel_28a and Dmel_32 were bitter receptors, Dmel_64a and Dmel_64F were carbohydrate receptors, and Harm_1 ∼ 3 were CO2 receptors. Spos_GR1, Spos_GR6, Spos_GR8, Spos_GR9, Nfla_GR2, Nfla_GR7, and Nfla_GR10 were clustered to CO2 receptor branches, while Spos_GR4 and Spos_GR5 were clustered to GR43a (Figure 2C and Supplementary Table 7). These results showed that the evolution of different taste receptors was relatively conservative, while neither bitter receptor nor sugar receptor was found in the two species.

The results of ENM showed areas suitable and unsuitable for the monophagous N. flavidorsalis and polyphagous S. postornata and the area of low to high habitat suitability for the corresponding host plants (Figure 4). The MaxEnt models using several feature classes (L, LQ, H, LQH, LQHP, and LQHPT) determine how predictor variables are transformed for each species or species combination. Based on the value of AUC_DIFF and delta AICc, the LQHPT was chosen as the best model. Generally, the suitable distribution area of the S. postornata was much larger than that of the N. flavidorsalis, which were remarkably consistent with the known distribution of these two species. Meanwhile, the suitable habitat distribution of N. flavidorsalis was uniform with the area of high suitability of its host plants, which indicated that the sympatric distribution of herbivorous insects and their host plants was obvious.