Insects were taken from a laboratory colony of the Laboratory of Entomology, Wageningen University and Research, Netherlands. Caterpillars were reared on Brussels sprouts plants, Brassica oleracea var. gemmifera cv. Cyrus, Brussels sprouts at 22 ± 2°C, RH 70 ± 10%, light 14 h: dark 10 h until pupation. After eclosion butterflies are transferred to another cage and fed with 5% sugar water at the same environmental condition. Brussels sprouts plants were provided to collect eggs.

Adult male and female butterflies were reared in the same cage for 3 days to ensure mating. Males and females were collected to dissect antennae for transcriptome analysis. Antennae from five male and five female butterflies were pooled for each biological replicate. Three biological replicates were collected for the antennal transcriptomes of P. brassicae and P. rapae.

Antennal samples were ground in Eppendorf tubes in liquid nitrogen and isolated and purified with ZYMO Quick-RNA/Insect Kit (ZYMO, Irvine, CA, USA). After the manual grinding, 450 μL lysis buffer was added. The RNA samples were isolated and purified according to the manual, followed by dissolving the RNA samples in 20 μL RNase-free H₂O through centrifugation. The final RNA samples were stored at −80°C. After that the RNA concentrations were determined using a DeNovix (Wilmington, NC, USA) spectrophotometer and RNA integrity was determined by denaturing gel electrophoresis. Paired-end sequencing (150 bp) was performed by Illumina NovaSeq6000 platform in the Novogene UK sequencing center (London).

Genomes of four different Pieris species: P. brassicae (GCA_942653925, MPI Chemical Ecology and GCA_905147105.1, Wellcome Sanger Institute), P. rapae (GCA_001856805.1) (Shen et al., 2016), P. napi (GCA_905475465.1), P. macdunnoughi (GCA_905332375.1), as well as the genomes of the outgroup butterflies Heliconius melpomene (GCA_000313835.2), and Danaus plexippus (GCA_009731565.1), Manduca sexta (GCF_014839805.1), Spodoptera litura (GCF_002706865.1) and Phoebis sennae (GCA_001586405.1) were collected to evaluate genome identity. The average nucleotide identity (ANI) was measured using fastANI v1.33 (Jain et al., 2018) using default settings. The detailed information of genomes that were used in this study is shown in Supplementary material 1.

The transcriptome was assembled via a genome-guided assembly method. The raw data of different transcriptomes were filtered by Fastp v0.20.0 (Chen et al., 2018) to remove short and low-quality reads and adaptors, the filtered data was subsequently assessed by FastQC v0.11.9 (Andrews, 2010). The filtered data was mapped to the corresponding genomes of P. brassicae (GCA_942653925) and P. rapae (GCA_001856805.1) with Hisat2 v2.1.0 by reporting alignments tailored for transcript assemblers (Kim et al., 2015), followed by transforming and sorting the mapped reads by SAMtools v1.10 (Li et al., 2009), the genome-guided assembly and quantification was achieved by StringTie v2.1.3b using default parameters (Pertea et al., 2015) and non-guided assembly was achieved by Trinity v2.8.6 (Grabherr et al., 2011). Transcripts per million (TPM) was employed to measure the expression levels of chemoreceptor genes in the antennal transcriptomes. In order to compare the expression of chemoreceptors between P. brassicae and P. rapae, we normalized the expression values of all the chemoreceptors to tubulin and rsps20, based on previous studies (Shakeel et al., 2015; Cusumano et al., 2018; Yin et al., 2020; Yang et al., 2021).

The protein sequences of ORs were collected from eight lepidopteran species: Bombyx mori (Tanaka et al., 2009), M. sexta (Koenig et al., 2015), S. littoralis (Walker et al., 2019), S. exigua (Zhang et al., 2018), H. melpomene (Heliconius Genome, 2012), Plutella xylostella (Engsontia et al., 2014), D. plexippus (Zhan et al., 2011), and Helicoverpa armigera (Liu et al., 2014), as queries to do homolog searching against the genomes of P. brassicae and P. rapae. The protein sequences of GRs were collected from nine lepidopteran species, i.e., B. mori (Wanner and Robertson, 2008), M. sexta (Koenig et al., 2015), S. littoralis (Walker et al., 2019), S. exigua (Zhang et al., 2018), H. melpomene (Briscoe et al., 2013), P. xylostella (Engsontia et al., 2014), and D. plexippus (Zhan et al., 2011; Briscoe et al., 2013), H. armigera (Liu et al., 2014) and Cydia pomonella (Walker et al., 2016), as queries for homolog searching. The protein sequences of IRs were collected from a total of 62 lepidopteran species, originating from three different datasets: IR sequences of 29 lepidopteran species were taken from the study by Liu et al. (2018) and IRs of 32 lepidopteran species were collected from the study by Yin et al. (2021), in addition, we also included the sequences of S. littoralis (Walker et al., 2019) as queries to search against the genomes of P. brassicae and P. rapae to find the candidate orthologs. Another butterfly species Ph. sennae (GCA_001586405.1) was introduced as an outgroup of Pieris to analyze the gene duplication cases.

Exonerate v2.2.0 (Slater and Birney, 2005) was applied to search against the genomes of P. brassicae (GCA_942653925) and P. rapae (GCA_001856805.1) with the collected queries to identify chemoreceptors in the two species, using the protein2genome model. Three best matched results were outputted for each query. The search results were sorted by scores. The targets sharing over 100 amino acids in a row and locating closely on the same scaffold were considered as the same receptor gene. The sequences with the highest scores were extracted, and the redundant alignment results and sequence lengths less than 100 amino acids were removed. The output results were collected as queries to further search in both genome-guided and non-guided transcriptome assemblies. The transcriptome searching results were then filtered by removing redundant and short results which are shorter than 100 amino acids. The transcripts were aligned to determine if some of them are from the same genes when they share 50 amino acids in a row. If any, the overlapped sequences will be assembled. The output transcripts and those outputs having no transcriptome targets were collected to search in NCBI database by blastp to determine if the outputs are full-length sequences. If not, the first high-similarity targets will be collected to repeat the previous steps until no longer sequence can be acquired. The identified chemosensory genes were also filtered with pfam domains (7tm_6 for ORs, Lig_chan for IRs, and 7tm_7 for GRs) to ensure the identified sequences are chemoreceptors. Diamond v0.9.25 (Buchfink et al., 2015) was later employed to identify potential chemoreceptors in the antennal transcriptomes by searching against these assemblies with an e-value of 10^–5 using the identified chemoreceptors from the genomes as queries.

The protein sequences of identified chemoreceptors were collected to search against the genomes of S. litura, Papilio glaucus (GCA_000931545.1), Ph. sennae (GCA_001586405.1), Vanessa tameamea (GCF_002938995.1), Delias pasithoe (GCA_010014985.1), and Calepholis nemesis (GCA_002245505.1) in addition to the P. brassicae and P. rapae genomes, to identify the potential gene duplication cases using Exonerate v2.2.0 (Slater and Birney, 2005). The searching was performed under the protein2genome model and only the best matched target for each query was reported.

The protein sequences of ORs, IRs and GRs from the seven lepidopteran species B. mori, M. sexta, S. littoralis, S. exigua, H. melpomene, P. xylostella, and D. plexippus were collected for phylogenetic analysis. The collected chemoreceptor sequences together with those of P. brassicae and P. rapae were aligned by MAFFT v7.475 (Katoh and Standley, 2013). The alignment output was imported into FastTree v2.1.11 (Price et al., 2010) for phylogenetic tree construction under the default settings with the maximum likelihood method. Olfactory receptors were rooted to Orco genes, while IRs were rooted to the iGluR clade and GRs were rooted to the sugar receptor clade. The output of FastTree was visualized and presented by the iTOL online tool v6.5.8 (Letunic and Bork, 2019).

The protein sequences of the identified chemoreceptors from the genomes were collected and imported to The MEME Suite version 5.4.1 (Bailey et al., 2015). The settings of discovery mode, sequence alphabet and site distribution were set to default, the number of motifs that should be found was set to 10. The gene structure information of the chemoreceptors was obtained by searching against the corresponding genomes (GCA_942653925 and GCA_001856805.1) with Exonerate v2.2.0 (Slater and Birney, 2005). The output of motif searching, gene structure as well as phylogenetic tree of chemoreceptors of the two Pieris species were integrated using TBtools v1.09876 (Chen et al., 2020). Transmembrane domains of all Piers chemoreceptors that we identified were predicted by TOPCONS (Tsirigos et al., 2015) and TMHMM 2.0 (Möller et al., 2001).

Average nucleotide identity was measured to evaluate the closeness of the phylogenetic relationship among butterflies and especially for P. brassicae and P. rapae. The P. brassicae genome of (GCA_942653925) has a very high identity with the latest released genome version of P. brassicae (GCA_905147105.1) by the Wellcome Sanger Institute of 99.14%, confirming the validity of both genome assemblies. P. brassicae shares relatively high similarities with the other three species belonging to the genus Pieris: 86.95, 86.84, and 86.01% with P. macdunnoughi, P. napi, and P. rapae, respectively, indicative of a very close relationship among the four Pieris species. The ANI decreases with increasing phylogenetic distance within the Lepidoptera (Table 1). For instance, the P. brassicae genome is less identical with D. plexippus (75.57%), H. melpomene (75.34%), S. litura (75.32%), and M. sexta (75.12%). Ph. sennae also belongs to the Pieridae and has a slightly higher ANI (76.23%) with the four Pieris species than those species belonging to other families.

All identified ORs (see Supplementary materials 2, 3) are named after D. plexippus because the ORs of P. brassicae and P. rapae share a relatively close relationship with D. plexippus according to phylogenetic analysis (Figure 1A). A total of 60 ORs were identified in both P. brassicae and P. rapae (Figure 1A) including the conserved odorant receptor coreceptor (Orco). The clustering of PbraOrco and PrapOrco with the Orco genes of other lepidopteran species was supported by high bootstrap values (Figure 1B). The reported pheromone receptors of moths such as B. mori, M. sexta, S. littoralis, S. exigua, and P. xylostella also cluster (Figure 1C), and show four HmelORs and one DpleOR as potential novel pheromone receptors. No pheromone receptors of P. brassicae and P. rapae were found in the pheromone receptor clade. Some gene duplication cases were found in the phylogenetic tree as well, such as the OR7/8/9 clade which is formed by the OR7, OR8, and OR9 of the two species with high bootstrap values.

Thirty-one PbraIRs and 34 PrapIRs (see Supplementary materials 2, 3) were identified from the P. brassicae and P. rapae genome, respectively. PbraIRs and PrapIRs were named after H. melpomene since HmelIRs are closer to most PbraIRs and PrapIRs in the phylogenetic tree compared to IRs from other species (Figure 2A). Sixteen out of the 31 PbraIRs and 17 out of the 34 PrapIRs were annotated in the antennal IRs clade (Figure 2A, yellow and green in the phylogenetic tree). Moreover, five PbraIRs (PbraIR8a, PbraIR25a, PbraIR25a, PbraIR76b, and PbraIR93a) and six PrapIRs (PrapIR8a, PrapIR25a1, PrapIR25a2, Prap76b1, PrapIR76b2, and PrapIR93a) were identified as ionotropic receptor co-receptors (Ircos) (Figure 2B). PrapIR76b1 and PrapIR76b2 share a very high similarity of more than 90%, their corresponding DNA sequences have a length of 2,072 bp on the same scaffold with the same direction, although there is only one IR76b gene found in P. brassicae genome. PrapIR25a1 and PrapIR25a2 share a lower identity of 76% with a 1,509 bp interval and, PbraIR25a1 and PbraIR25a2 share a 75% similarity with a 10,008 bp interval on the same scaffold with the same direction. Fifteen out of 31 PbraIRs and 17 out of 34 PrapIRs were annotated in divergent IRs clade in blue in the phylogenetic tree (Figure 2). The orthologs of PrapIR7d.4 and PrapIR100e were not found in the P. brassicae genome. Thirty-two out of 34 PrapIRs have been identified previously (Liu et al., 2018), the two newly identified genes are PrapIR4 and PrapIR25a2.

Thirty-nine GRs (see Supplementary materials 2, 3) were identified in both the P. brassicae and P. rapae genomes. All identified GRs were named after H. melpomene because most PbraGRs and PrapGRs clustered with HmelGRs (Figure 3A), just like the IRs presented above. Based on former publications (Wanner and Robertson, 2008; Yang et al., 2020), three PbraGRs and PrapGRs were annotated as CO₂ receptors (Figure 3A, in yellow), three PbraGRs and PrapGRs were annotated as fructose/inositol receptors (Figure 3A, in blue) and nine PbraGRs and PrapGRs were annotated as sugar receptors (Figures 3A, B in pink) and a sinigrin receptor of P. rapae (Yang et al., 2021) was indicated by the blue arrow in Figure 3A.

Based on the chemoreceptor sequences identified from the genomes of the two Pieris species, we compared the full-length gene structures and analyzed the conserved protein motifs and the arrangements of motifs. We found that most chemoreceptors have a comparable gene structure with other orthologs and with the same number of introns but sometimes with different lengths. For example, PrapOR58 and PbraOR60 have much longer introns than their orthologs (Figure 4). Interestingly, we found that most of IRs in the divergent IRs clade annotated in Figure 3 have no intron at all, such as IR7d.3, IR100c, IR100d, IR100e, IR100f, IR87a, IR7d.2, IR7d.2.1, IR7d.2.2, and IR143; three IRs, i.e., IR85a, IR100a, and IR7d.4, have only one very short intron. Notably, the orthologs having no or only one short intron share highly similarly conserved motifs and arrangements. The other IRs in the divergent IRs clade, i.e., IR1.1, IR1.2, IR2, and IR4, have multiple introns similar to those in antennal IRs clade (Figure 5). However, we also found some exceptions in GRs that showed different number of introns, specifically PbraGR1 has 8 introns while PrapGR1 has 9 introns (Figure 6).

Multiple motifs (sequences see Supplementary material 4) are found in the chemoreceptor protein sequences. Motifs 1, 2, 3, 4, and 5 are the most abundant conserved motifs and are also found in most ORs at similar position of the sequences in both species (Figure 4). According to our transmembrane domain prediction results (Supplementary material 5), we found that ORs exhibit six to seven transmembrane domains. Interestingly, the location of Motif 1 is highly consistent with that of the last transmembrane (TM) domain 7, location of Motif 6 overlaps with TM 6, and Motif 5 covers TM 5 in most ORs. The other conserved motifs are located in intracellular or extracellular loops.

In the IRs Motif 8 is present in all sequences, while Motif 2 is found in most IR sequences. In contrast, Motif 6 is only present in sequences of IR25a1/2 of both species (Figure 5), which is located in amino-terminal domain (ATD). According to the transmembrane predictions, IRs have three to four transmembrane domains. Here, locations of identified conserved motifs are highly consistent with known domains in IRs (Benton et al., 2009; Croset et al., 2010). Motifs 7, 9, and 8 are TM 1, TM 2, and TM 3, respectively. Motif 4 overlaps with the pore domain (P). In addition, Motif 2 was found to locate in ligand-binding domain (LBD) S1 and, Motif 1, and Motif 3 were found to locate in LBD S2.

Compared to ORs and IRs which have four to nine conserved motifs, GRs have fewer conserved motifs in their sequences; for instance, no conserved motif has been found in GR6 and over half of the GRs have only 1–2 conserved motifs (Figure 6). The most common motifs 1 and 2 are only found in about half of the GR sequences. Transmembrane domain prediction results showed that most GRs have six to eight transmembrane domains. Locations of Motif 1 and Motif 2 match well with TM 7 (or the last transmembrane domain) of GRs. Motif 3 and Motif 8 covers TM 5 and TM 4, respectively while the two conserved motifs did not strictly match with the transmembrane domains of GRs.

Based on the identification and phylogenetic analysis of chemoreceptors, we found some gene duplication cases. Clade OR18/OR26/OR27, clade OR48/OR49/OR50, and clade OR7/OR8/OR9 are positioned closely on the same scaffold and share high similarities of protein sequences in both Pieris species. Similar cases can also be found in IRs and GRs. IR25a1/IR25a2 and IR100d/IR100f of both species and PrapIR76b1/PrapIR76b2 are placed closely on the same scaffold and show a very high sequence identity. This was also found for GR8.1/GR8.2/GR8.3. Interestingly, while most chemoreceptors mentioned above have the same number of exons, OR49/OR50 have eight exons and are different from OR48 which has ten exons.

In order to clearly elucidate the potential gene duplication and microsyntenic pattern of chemoreceptors in Pieris species, we selected OR7/OR8/OR9 (Figure 7A) and IR25a1/IR25a2 (Figure 7B) as targets to extract their relative positions and gene structures in P. brassicae, P. rapae, P. napi, and P. macdunnoughi. We found that the four selected Pieris species have these target genes, and these paralogs are closely located on the same scaffold in the same order and orthologs share the same number of exons and even share similar exon lengths. In contrast, the ortholog of OR7 is not present in the Pieris outgroup species Ph. sennae and no IR25a duplication was found in any of the selected non-Pieris lepidopteran species (S. litura, Pa. glaucus, V. tameamea, Ph. sennae, D. pasithoe, and C. nemesis).

We compared the chemoreceptor identification results of the de novo assembly and genome-guided assembly and found the results are identical, although de novo assembly did not always reveal the complete sequences. Sixty PbraORs and 60 PrapORs are identified in the adult antennal transcriptome of P. brassicae and P. rapae, respectively, Orco is the highest expressed gene among chemoreceptors in the transcriptomes. However, out of the 60 genes, the expression levels of 25 ORs were found to be significantly different between the transcriptomes of the two Pieris species. Twenty-two of the 25 differently expressed ORs showed higher expression levels in the P. rapae antennal transcriptome, OR27, OR47, and OR54 had higher expression levels in the P. brassicae antennal transcriptome (Figure 8).

Twenty-nine and 31 IRs were identified from the P. brassicae and P. rapae adult antennal transcriptomes, respectively. Ircos PbraIR8a, PbraIR25a1, PbraIR93a and those IRs in the antennal IRs clade PbraIR41a, PbraIR75d, PbraIR75q.1, and PbraIR85a were highly expressed in the antenna of P. brassicae. Out of these, 12 IRs were found differently expressed between the two transcriptomes. Seven (IR1.2, IR7d.3, IR8a, IR25a1, IR25a2, IR75p.2, and IR100c) out of the 12 IRs were found to be higher expressed in the P. rapae antennal transcriptome, such as IR8a, IR25a1, and IR25a2 than P. brassicae and five of the 12 IRs (IR1.1, IR41a, IR64a, IR75q.1, and IR93a) had higher expression levels in the P. brassicae antennal transcriptome. In addition, IR76b1, IR76b2, and IR 100d/e/f were not found in the antennal transcriptome of P. brassicae; some of these genes (IR76b1, IR76b2, and IR100e) were only found in the P. rapae genome. IR76b, IR100a, and IR143 are exclusively found in the P. brassicae transcriptome (Figure 8).

In total, 33 GRs were identified from both the P. brassicae and the P. rapae adult antennal transcriptome, although several genes were only identified in one of the two species. Compared to ORs and IRs, GRs have comparatively lower expression in the adult antennal transcriptomes of both species according to normalized TPM values. The sinigrin receptor GR63 was highly expressed in both antennal transcriptomes of P. brassicae and P. rapae. Ten GRs were found differently expressed between the two transcriptomes. Six (GR3, GR8.1, GR8.2, GR8.3, GR9, and GR16) of these ten GRs showed higher expression in the P. rapae antennal transcriptome, the other four GRs (GR8.4, GR30, GR61, and GR63) showed higher expression in the P. brassicae antennal transcriptome. Notably, four PbraGRs (PbraGR7, PbraGR42, PbraGR44.1, and PbraGR44.2) and four PrapGRs (PrapGR28, PrapGR29, PrapGR68.1, and PrapGR68.2) were not detected in the corresponding antennal transcriptomes. These undetected GRs are displayed in white cells in Figure 8.