The colony of L. erysimi utilized in this study was established in 2020, based on the field populations from Xinfeng in the Jiangxi province of China (Kuang et al., 2023). The population was continuously maintained on Chinese cabbage Shanghaiqing (Brassica rapa var. chinensis) without exposure to any insecticides under controlled conditions (27°C–28°C, 60%–65% RH, 14:10 L:D photoperiod).

Samples of L. erysimi were collected at five different stages, which include 1st instar nymph, 2nd instar nymph, 3rd instar nymph, 4th instar nymph, 1-day apterous and alate adults. The apterous aphids anesthetized on ice were dissected under a stereomicroscope to collect various tissues of L. erysimi, such as antenna, head, leg, and cuticle. Total RNA extraction was performed according to the manufacturer’s protocol of Trizol reagent (Sigma, St. Louis, MO, United States). The purity and concentration of RNAs were determined using a NanoDrop One^C spectrophotometer from Thermo Fisher Scientific (Waltham, MA, United States). The first-strand cDNA synthesis was synthesized with 500 ng of purified RNA, utilizing the highly effective EasyScript^® One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TransGen Biotech, Beijing, China). The resulting cDNA was properly stored at −20°C until ready to be used.

To isolate the antennal-biased genes in L. erysimi, the high-throughput transcriptome data sets were retrieved from our former study (GenBank accession number PRJNA947784), which included conducting Illumina sequencing on the antennae and bodies (excluding antennae) of adult apterous aphids, transcriptome de novo assembly, as well as functional annotation of the unigenes (Kuang et al., 2023). The differential gene expression analysis was carried out using the antennal and body transcriptome data as described by Yu et al. (2023). Briefly, the transcript abundances were determined by RSEM (version 1.2.12); DESeq2 (version 1.4.5) was utilized to identify differentially expressed genes (DEGs) between samples, and a gene was considered differentially expressed if the corrected p-value was ≤ 0.05.

The antennae-biased ODE genes with FPKM ≥10 were selected from the gene repertories obtained through comparative transcriptome analysis. Candidate ODEs were confirmed using the BLASTX algorithm, and their open reading frames (ORFs) were predicted using the ORF finder tool (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). The amino acid components, theoretical isoelectric points (pIs), and molecular weights (MWs) of ODE genes were calculated using ExPASy (http://web.expasy.org/protparam/). The deduced protein sequences were submitted to the SignalP 5.0 server (https://services.healthtech.dtu.dk/services/SignalP-5.0/) for the prediction of the signal peptide sequences and their corresponding cleavage sites. Conserved domains were predicted with CDD-BLAST (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and InterProScan (https://www.ebi.ac.uk/interpro/) servers.

The dataset submitted for phylogenetic analysis comprised the candidate CXE, CYP, GST, and UGT genes of L. erysimi in this study, their homologs from Acyrthosiphon pisum (Ramsey et al., 2010), Aphis craccivora (Yang et al., 2021), Aphis gossypii (Pan et al., 2018), Myzus persicae (Pan et al., 2019), Nilaparvata lugens (Vontas et al., 2002; Zhou et al., 2013), Diaphorina citri (Yu and Killiny, 2018; Tian et al., 2019; Wu et al., 2020; Kuang et al., 2022), Anopheles gambiae (Ding et al., 2003), and Bombyx mori (Yu et al., 2008); in addition to some well-identified ODE genes, including the studied antennal CXEs (Ishida and Leal, 2005; 2008; Durand et al., 2010; Durand et al., 2011; He et al., 2015; Wei et al., 2021), CYPs (Keeling et al., 2013; Chiu et al., 2019a), GSTs (Rogers et al., 1999; Li et al., 2018; Liu et al., 2021; Xia et al., 2022), and UGTs (Wang et al., 1999; Bozzolan et al., 2014). The protein sequences were first aligned using CLUSTAL_X version 1.83. The joint unrooted phylogenetic tree was constructed with MEGA11 using the neighbor-joining method (Tamura et al., 2021). Branch support was evaluated through the bootstrap method which consisted of 1,000 replicates. The phylogenetic tree was visualized using iTOL web tool (https://itol.embl.de/).

The quantitative real-time PCR (qRT-PCR) reactions were conducted in a 20 μL volume that comprised of 4 μL of diluted cDNA, 0.4 μM of each primer, and 10 μL of PerfectStart^® Green qPCR SuperMix (TransGen Biotech, Beijing, China). Reactions were performed on a Roche LightCycler 96^® system (Roche Diagnostics, Mannheim, Germany) with the following thermal program: initial denaturation for 10 min at 95°C, followed by a 40-cycle two-step amplification profile of 95°C for 5 s and 60°C for 30 s. Two reference genes, actin (GenBank accession number OQ626608) and GAPDH (GenBank accession number OQ626609), were employed to standardize the quantity of cDNA added to the PCR reactions. The relative expression of each ODE gene was analyzed using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). Reactions were performed in triplicate, and the gene-specific primers are listed in Table 1.

The statistical differences of ODE gene expression levels among different developmental stages and tissues were analyzed using analysis of variance (ANOVA), followed by Tukey multiple comparison test. The statistical analysis was conducted using GRAPHPAD PRISM software (version 6.0; GraphPad Software Inc., La Jolla, CA, United States) and a p-value of ≤ 0.05 was set as the threshold for statistical significance.

Comparative analyses of the antennal and body transcriptomes in this study provide useful information to identify the antennae-abundant and/or antennae-biased genes. A total of 8,932 differentially expressed genes (DEGs) with a Q value ≤ 0.05 were identified, and 4,797 DEGs were significantly upregulated in the antennae (Figure 1). Among the antennae-abundant DEGs, one CXE, seven CYPs, two GSTs, and five UGTs were identified by blasting against the Nr database. The candidate ODEs were designated according to gene names of the top blast hits in NCBI. Detailed information on these ODE enzymes is shown in Table 2.

One putative CXE gene, LeCXE6, showed a high level of expression in the antennae, with FPKM values over 10 times higher than in the rest of the body. LeCXE6 encodes a 564 amino-acid protein, with a signal sequence cleavage site predicted between Gly-18 and Phe-19 at the N-terminus. The predicted protein has a theoretical molecular mass of 62.73 kDa and an isoelectric point of 5.57, as determined using ProtParam tool in Expasy server (Table 2). Conserved domain and sequence alignment analysis revealed that LeCXE6 contained the typical motif of carboxylesterase family, including a conserved pentapeptide (Gly-X-Ser-X-Gly) and an oxyanion hole (Gly-Gly-Ala; Supplementary Figure S1). Phylogenetic analysis showed that LeCXE6 fell into the beta-esterase clade and was closely clustered with the well-studied odorant-degrading enzymes, PjapPDE and ApolPDE (Figure 2).

Seven DEGs abundant in L. erysimi antennae were identified to be CYPs by blasting against the Nr database. All candidate LeCYPs were found to contain full-length ORFs without any predicted signal peptide sequences. Their antennal RPKM values ranged from 29.83 to 1,220.17, representing more than a ten-fold increase in comparison to the body group. Among them, LeCYP6a13 was the most abundant CYP gene in antennae with a value exceeding 1,200, followed by LeCYP18a1 and LeCYP6a2 (Table 2). Phylogenetic analysis showed that the selected CYPs were well categorized into four subclasses, namely, CYP2, CYP3, CYP4, and mitochondrial CYP. The CYP3 class comprised of five antennal LeCYPs, including LeCYP6a13, LeCYP6a2, LeCYP6k1, LeCYP6a14.1, and LeCYP6a14.2. LeCYP18a1 was classified as a member of the CYP2 clan, while LeCYP4c1 was categorized into the CYP4 clade (Figure 3).

Two DEGs abundant in antennae were identified to be GSTs. Both LeGST1 and LeGST transcripts had full-length ORFs, encoding proteins that are 157 and 198 amino acids, respectively. It is noteworthy that LeGST1 showed an expression pattern that was particularly abundant in antennae, with an FPKM value of 93.81 that exceeded 18-fold higher than in the body (Table 2). Conserved domain analysis showed LeGST1 had a MAPEG (membrane-associated proteins in eicosanoid and glutathione metabolism) domain and was similar to the microsomal GST1, while LeGST had one GSH binding site (G-site) in the N-terminus and one hydrophobic substrate binding pocket (H-site) in the C-terminal region (Supplementary Figure S2). The phylogenetic analysis revealed that eight subclasses, namely, Microsomal-, Delta-, Epsilon-, Omega-, Sigma-, Theta-, Zeta-, and the unclassified-GST, were well clustered in their respective phylogenetic branches. LeGST1 was classified under the Microsomal-GST subclass, and LeGST was classified as a member of Sigma-GST (Figure 4).

A total of five antennal LeUGT genes were identified in the utilized transcript set. All LeUGT transcripts had full-length ORFs, encoding proteins ranging from 513 to 542 amino acids. Signal peptides were predicted in all candidate LeUGTs, with the exception of LeUGT2B13. FPKM analysis showed that LeUGT2B7 was the most antennae-abundant UGT with a value of 404.61, which was >20-fold higher than in the body (Table 1). Multiple alignments revealed the UGT motif signature sequence, (FVA)-(LIVMF)-(TS)-(HQ)-(SGAC)-G-X (2)-(STG)-X (2)-(DE)-X (6)-P-(LIVMFA)-(LIVMFA)-X (2)-P-(LMVFIQ)-X (2)-(DE)-Q, was situated at the C-terminus of LeUGTs (Supplementary Figure S3). Phylogenetic analysis showed that the candidate LeUGTs were grouped into three distinct subclades, with each subclade including several homologs from other aphid species. Specifically, LeUGT2B7, LeUGT2B13, and LeUGT2C1.2 were categorized into the UGT344 clade; LeUGT2C1.1 was clustered in the UGT343 subclade, and LeUGT1-7 was found to be a member of the UGT351 subclade (Figure 5).

The developmental and tissue expression profiles of LeCXE, LeCYP, LeGST, and LeUGT genes were analyzed using qRT-PCR. Developmental expression data showed that the candidate ODE genes were consistently detected throughout the various developmental stages of L. erysimi, spanning from the first instar nymph to the adult stage. Notably, LeCYP4c1, LeCYP6a2, LeCYP6a13, LeCYP6a14.2, LeCYP18a1, LeUGT2B7, and LeUGT2B13 exhibited significantly higher expression levels in alate aphids compared to apterous and nymph aphids (Figure 6). Tissue expression analysis revealed that LeCYP6a14.1 and LeGST were highly expressed in both antenna and gut tissue, while the remaining 13 ODE genes displayed antennae-enriched expression profiles. In particular, the antennal expression levels of LeCYP6a13, LeCYP6k1, LeCYP6a14.2, LeGST1, LeUGT2B13, and LeUGT2C1.2 were >10-fold higher than in other tissues; LeCXE6, LeCYP4c1, LeCYP6a2, LeCYP18a1, LeUGT2B7, and LeUGT2C1.1 exhibited more than four times higher expression in antennae compared to non-olfactory tissues such as the head, leg, gut, and cuticle (Figure 7).