Diversification and selection pattern of CYP6B genes in Japanese Papilio butterflies and their association with host plant spectra

Papilio host plant data

Host plant data for nine Papilio species (P. dehaanii, P. maackii, P. xuthus, P. machaon, P. helenus, P. memnon, P. macilentus, P. protenor, and P. polytes) were collected from the public database InsectInDB (http://insect-plant.org). Since there was an inflated number of Citrus cultivars or related species in the lists compared to wild Rutaceae plants, the number of host plant of each Papilio species were biased when the species feed on Citrus plants and did not highlight their diet breadth. Therefore, we performed principal component analysis (PCA) on host plant data from Rutaceae specialist Papilio species to identify potential differences among their host plant spectra. The PCA was performed with prcomp() function in R software (R Core Team, 2019) to acquire PC1 and PC2 scores as indicators of the host plant spectrum of each Papilio species. We excluded P. machaon from this analysis because it feeds exclusively on Apiaceae.

Rutaceae furanocoumarin profiling and comparison with Papilio host plant spectra

To investigate the furanocoumarin profiles of Rutaceae plant species, we collected undamaged leaves from 13 Rutaceae plant species from wild and cultivar plants listed as Papilio hosts in the host plant database (Citrus depresa, C. junos, C. limon, C. trifoliata, C. unshiu, Orixa japonica, Phellodendron amurense, Skimmia japonica, Toddalia asiatica, Zanthoxylum ailanthoides, Z. armatum, Z. piperitum, and Z. schinifolium; http://insect-plant.org). The collected fresh leaves were immediately frozen at −20 °C and freeze-dried for further chemical analysis. We used liquid chromatography–electrospray ionization mass spectrometry to quantify furanocoumarins in the sampled leaves. We ground 20 mg freeze-dried leaves using metal balls in 2-mL tubes and added one mL 80% methanol and 2.5 µM lidocaine as an internal standard. We centrifuged the samples at 9,000 rpm for 3 min and analyzed the collected supernatants. We used a bridged ethyl hybrid (BEH) C18 column (Acquity C18; 1.7 × 2.1 × 100 mm, Waters, Milford MA, USA) for reversed-phase LC and set soluble A as 0.1% HOOH water and soluble B as 0.1% HCOOH acetonitrile. The flow rate was 0.3 mL/min, using a program of 0.5% B (0–1 min), 0.5–80% B (1–12 min), 80–99.5% B (12–15 min), and 99.5–0.5% B (15–20 min). We extracted peaks with m/z ratios of 100–400 and identified furanocoumarins with eight standards (angelicin; CAS 523-50-2, bergapten; CAS 484-20-8, imperatorin; CAS 482-44-0, isobergapten; CAS 482-48-4, isoimperatorin; CAS 482-45-1, isopimpinellin; CAS 482-27-9, psoralen; CAS 66-97-7, and xanthotoxin; CAS 298-81-7). The acquired peaks and retention times of each standard were used to identify and quantify furanocoumarins in the samples. We used the Xcalibur software (Thermo Fisher Scientific, Waltham, MA, USA) to extract peaks and identify furanocoumarins; we removed peaks lower than 0.001. Each plant species was analyzed in triplicate.

We compared the total amount of detected furanocoumarins in each host plant among plant species. Since different furanocoumarins can exert different toxicity against herbivores, we also compared furanocoumarin diversity among plant species using the chemical complexity index (CCI, Becerra, Noge & Venable, 2009); the CCI was calculated based on two Shannon indices of chemical diversity, one based on presence/absence and the other based on relative concentration. We calculated these two Shannon indices based on furanocoumarin profiles of each plant species and CCI was acquired by totaling these two indices (Becerra, Noge & Venable, 2009). For each Papilio species, we calculated the average level of detected furanocoumarin and the CCIs of their host plants. We performed linear regression analyses to compare these values with the PC1/2 scores, which showed the host plant spectrum of each Papilio species. This analysis allowed us to examine the correlation between differences in Papilio host plant spectra and furanocoumarin profiles among host plants.

Larval sampling for transcriptome analyses

We collected eggs or fertilized female butterflies from wild populations of Papilio species in Japan (Table S1). Fertilized females were kept with suitable host plants for egg laying. Hatched neonate larvae were reared with their host plants until the second instar. Five individuals of each Papilio species were dissected and their guts were extracted for transcriptome analysis. We extracted RNA using the RNeasy Mini Kit (Qiagen, Hilden, Germany) and pooled one unit of RNA from each of five individuals for sequencing. Genomic DNA in the samples were digested with the TURBO DNA-free Kit (Thermo Fisher Scientific, Waltham, MA, USA). After quality confirmation using the Agilent 2100 Bioanalyzer, samples were used for library preparation. We performed 100-bp pair-end sequencing using the Illumina HiSeq 4000 system.

De novo transcriptome assembly, CYP6B gene extraction, and phylogenetic analysis

We controlled the quality of the raw reads using the Trimmomatic ver. 0.32 command line tool, with the settings LEADING:10, TRAILING:10, SLIDINGWINDOW:4:20, and MINLEN:40 (Bolger, Lohse & Usadel, 2014). After trimming, the read quality was verified with FastQC software then de novo assembly was performed with the Trinity ver. 2.1.1 software (Grabherr et al., 2011; Haas et al., 2013). We extracted the longest isoforms from the assembly using the Trinity command “get_longest_isoform_seq_per_trinity_gene.pl” then performed tblastn (Altschul et al., 1990; Camacho et al., 2009) on the assembled contigs using CYP6B genes of Papilio species from previous studies as queries (accession nos. AAB06742.1, AAB06743.1, AAK69477.1, AAK69478.1, AAK69497.1, AAK69499.1, AAK69500.1, AAK69503.1, AAK69504.1, AAK69494.1, AAK69495.1, AAK69496.1, AAK69498.1, AAK69501.1, AAK69505.1, and AAB06741.1) with the e-values set at 0.0001. After extracting hits from the tblastn search, sequences shorter than 300 bp were eliminated. We translated the extracted contigs to amino acid sequences and aligned them with reference sequences using the mafft tool (Katoh & Standley, 2013), with the –auto option. Lepidopteran CYP sequences found in the P450 database (http://drnelson.uthsc.edu/cytochromeP450.html) were used as additional reference sequences in this analysis. We reconstructed a maximum likelihood (ML) phylogeny using the IQ-TREE software with 1000 bootstrap replicates (Nguyen et al., 2015). Based on the resulting phylogeny, we excluded genes that were not included in the CYP6B clade and assigned sequences with >  55% amino acid identity to the reference as CYP6B genes. To observe the patterns of CYP6B gene duplication and loss among tested Papilio species, we conducted gene-tree species-tree reconciliation analyses with NOTUNG software (Stolzer et al., 2012). Lower supported nodes (<80% in bootstrap values) were rearranged and the gene tree were reconciled along with the species tree generated from the transcriptome data (see Species phylogeny reconstruction section below).

To see the relationships between host plant spectrum and CYP6B duplication patterns among Papilio species, we also analyzed correlations between the observed CYP6B numbers and host plant spectrum (PC1/2 scores and average host plant CCI or furanocoumarin amount) of each Papilio species by linear regression.

As an ad hoc analysis, we also analyzed available genome data from P. machaon (NCBI: GCA_001298355) and P. momnon (NCBI: GCA_003118335) to confirm the absence of a CYP6B gene with a signature of positive selection (see Results). We performed the analyses described above using the assembled genome and extracted CYP6B genes. Unaligned intron regions of the extracted genes were trimmed out in the phylogenetic analysis. We performed the same phylogenetic analysis again, including CYP6Bs extracted from the genomes.

Species phylogeny reconstruction

BUSCO single-copy gene sets were extracted and used to reconstruct Papilio species phylogeny from the transcriptome data by running BUSCO (Seppey, Manni & Zdobnov, 2019) on each transcriptome assembly. When duplicate hits occurred, we extracted the longer contig as the representative. We extracted BUSCO genes that were found in all species, aligned the sequences using the mafft tool, and concatenated the alignment for phylogenetic analyses. The concatenated alignment was de-gapped with using the TrimAl software (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) and IQ-TREE was used to reconstruct ML phylogenies of the nine Papilio species from the concatenated sequences, with substitution model selection and 1000 bootstrapping iterations. RNA sequencing data from four pierid butterflies (P. brassicae, P. canidia, P. melete, and P. napi) were used as outgroups (EBI Accession numbers: ERX2829492–ERX2829499, ERX3552761).

Detection of positive selection

We performed branch-site model tests to examine patterns of positive selection among CYP6B gene family members in Papilio species using the codeml program implemented in the PAML software package (Yang, 2007). We excluded outgroups and reconstructed a ML tree (nucleotide) for CYP6B using IQ-TREE. All the major internal branches having >70% node support according to the bootstrap values were tested. We used model 2, with NSsites = 2, and ran an alternative model, which allows varied the non-synonymous substitution (dN) to synonymous substitution (dS) ratios (dN/dS) across sites and lineages, as well as a null model with a fixed dN/dS ratio (fixed_omega = 1). We compared the results of the two models using the likelihood ratio test (LRT) with a chi-square distribution to detect significant differences between the alternative and null models. P values were adjusted with false discovery rates. We performed subsequent Bayes empirical Bayes analysis (0.95 cut-off) to identify sites potentially under positive selection once the alternative model was selected by the likelihood ratio test. The gene tree topology with lower node supports might affect the results, therefore, we also performed the same branch-site model tests using a gene tree rearranged and reconciled by NOTUNG along with the species tree. The same branches were selected for the test if the branch was still conserved after NOTUNG gene tree rearrangement and reconciliation processes to confirm the results of branch-site model tests based on the ML tree.

Host plant spectrum associated with host plant furanocoumarin profiles

We performed principal component analysis to compare the host plant spectra of eight Rutaceae-feeding Papilio species (Fig. 1). Although the host plant spectra overlapped among Papilio species, P. memnon and P. polytes tended to use more Citrus species, whereas P. dehaanii, P. protenor, and P. macilentus were more reliant on Orixa japonica, Skimmia spp., or Zanthoxylum spp. We used these PC1/2 scores as indicators of the host plant spectra of these Papilio species in later analyses.

The chemical analysis showed that Orixa japonica and Skimmia japonica had higher levels of furanocoumarins than those of Citrus species (Fig. 2, Table S2). Bergapten was dominant in Orixa japonica, and isoimperatorin was dominant in Skimmia japonica. These findings were consistent with those of previous studies (Atkinson, Boyd & Grundon, 1974; Zobel & Brown, 1990). We detected lower levels of furanocoumarins in the leaves of Citrus spp. than in those of other species. Although higher furanocoumarin levels were observed in Citrus peels or juice in previous studies (Barreca et al., 2011; Dugrand et al., 2013), lower levels in Citrus leaves were also reported in a previous study (Durand-Hulak et al., 2015). Based on this analysis, we calculated CCIs to assess the furanocoumarin profile complexity of each Rutaceae species. Orixa japonica and Zanthoxylum ailanthoides had higher CCI values (Table S2).

We compared Papilio host plant spectra with furanocoumarin profiles of the host plants. The average level of detected furanocoumarin was higher in host plants of P. macilentus and P. maackii and lower in those of P. memnon and P. polytes (Table S3). There were no apparent differences in the furanocoumarin CCI among the host plants of each Papilio species (ANOVA; P = 0.969) (Table S3). We found a significant negative correlation between the average level of detected furanocoumarin in host plants and PC1, which showed the host plant spectrum of each Papilio species (Fig. 3A). However, we did not find a significant relationship between the average CCI of the host plants of each Papilio species and the PC axes (Fig. 3B). Thus, the level of furanocoumarins but not the complexity differed among host plants of each Papilio species.

No specific diversification patterns of CYP6B correlating to patterns of host plant spectrum or species phylogeny

The transcriptome assembly statistics are shown in Table S4. We extracted CYP6B-related genes from the de novo transcriptome assembly of each Papilio species and performed phylogenetic analysis. The extracted genes included 33 CYP6B genes from nine Papilio species. Figure 4A shows an ML phylogeny of CYP6B with reference sequences. We also reconstructed a tree of these Papilio species using 858 extracted BUSCO genes (Fig. 4B, Fig. S1B) that were found across all transcriptome data and aligned for phylogenetic analysis. The length of this alignment was 567,026 bp after gap elimination.

Generally, we did not detect species- or lineage-specific diversification patterns of CYP6B (Fig. 4 and 4B). The gene tree–species tree reconciliation analysis showed that the diversification of CYP6B genes among Papilio species occurred in the most ancestral Papilio species of those tested species, followed by subsequent gene loss in each lineage (Fig. 4B). We detected no apparent diversification pattern associated with the host plant spectra (Fig. 4A). P. machaon, which exclusively uses Apiaceae, not Rutaceae, did not exhibit a divergent CYP6B clade, which was also observed in the Rutaceae-feeding species. We also compared CYP6B gene number of each species and their host plant spectrum and this did not show any significant correlations (Fig. S2).

Signature of positive selection on CYP6B branches among Papilio species

The branch-site model tests with both CYP6B ML tree and reconciled CYP6B tree showed that branch1, which was closer to furanocoumarin-inducible CYP6B references (P. canadensis CYP6B18 and P. glaucus CYP6B4 in Fig. 4A), displayed a signature of positive selection (Fig. 4A, Fig. S1A, Table 1, Table S5). Branch 1 involved CYP6Bs from most of the Papilio species analyzed, except for P. machaon and P. memnon. None of the species had duplicated CYP6Bs associated with this branch in our transcriptome data. At ML tree-based analysis, we found evidence of positive selection on 10 sites at this branch and 8 of them were also found in the same analysis with reconciled CYP6B tree (Table 1, Fig. S1). These positively selected sites included substrate recognition sites 3 (SRS3) and SRS6, which may have a strong impact on the substrate specificity of CYP6Bs (Li, Schuler & Berenbaum, 2003) (Table S5). As an ad hoc analysis, we also searched for these particular CYP6B genes within the genomes of P. machaon and P. memnon and did not find these genes within this particular branch (Fig. S4). At ML tree-based analyses, we also found a signature of positive selection on branch 6 (Fig. 4A, Table 1), which included CYP6Bs from P. xuthus, P. helenus, and P. memnon. One site was also located at SRS2 (Fig. S4) and found to be under positive selection at this branch (Table 1). However, the significance at this branch disappeared in the analysis with reconciled tree (Fig. S3).