The Birth-and-Death Evolution of Cytochrome P450 Genes in Bees.

The birth-and-death model of multigene family evolution describes how gene families evolve and diversify through duplication and deletion. The cytochrome P450s are one of the most diverse and well-studied multigene families, involved in both physiological and xenobiotic functions. Extensive studies of insect P450 genes have demonstrated their role in insecticide resistance. Bees are thought to experience toxin exposure through their diet of nectar and pollen, as well as the resin-collecting behavior exhibited by some species. Here, we describe the repertoire of P450 genes in the orchid bee Euglossa dilemma. Male orchid bees form perfume bouquets used in courtship displays by collecting volatile compounds, resulting in exposure to compounds known to be toxic. In addition, we conducted phylogenetic and selection analyses across ten bee species encompassing three bee families. We find that social behavior and resin collection are not correlated with the repertoire of P450 present in a bee species. However, our analyses revealed that P450 clades can be classified as stable and unstable, and that genes involved in xenobiotic metabolism are more likely to belong to unstable clades. Furthermore, we find that unstable clades are under more dynamic evolutionary pressures and exhibit signals of adaptive evolution. This work highlights the complexity of multigene family evolution, revealing that multiple factors contribute to the diversification, stability, and dynamics of this gene family. Furthermore, we provide a resource for future detailed studies investigating the function of different P450s in economically important bee species.

Significance

Gene family evolution is characterized by deletion and duplication, but it is unclear if this is driven by adaptive or neutral evolutionary processes. To investigate the dynamics of gene family evolution we analyze the P450s of ten bee species, a well-studied family involved in both physiological and detoxification roles. We find that genes involved in detoxification are more likely to belong to unstable clades, exhibit dynamic evolutionary pressures, and show signals of adaptive evolution. We do not find evidence for a relationship between resin collection or sociality and the P450 repertoire of a bee species. Our findings reveal that multiple factors contribute to the dynamic evolution of P450s and provide a resource for investigating P450s in economically important bee species.

Introduction

Multigene families arise from the duplication of a common ancestral gene, resulting in groups of genes that share similar sequences, and often functions (Nei and Rooney 2005). To understand how these families evolve, early studies focused on concerted evolution in rRNA, whereby members of a gene family evolve together due to repeated unequal crossover events (Brown et al. 1972). More recently, an alternative model, the birth-and-death model, has been used to explain gene family evolution. In this model, genes evolve independently and expansion and contraction occurs through gene duplication, formation of pseudogenes, and gene deletion (Nei and Rooney 2005). The relative roles of neutral versus adaptive evolutionary forces within the framework of the birth-and-death model has proven more complex (Eirín-López et al. 2012).

Cytochrome P450s are one of the best-studied and most diverse multigene families (Nelson et al. 2013). They are enzymes that use molecular oxygen to change the structure of their substrates, a reaction that has an important role in both physiological endogenous processes, and also ecological and xenobiotic processes (Nelson et al. 2013). P450s are highly diverse and groups have been named based on sequence similarity, with four prominent groups described in insects (Nelson 1998; Tijet et al. 2001; Dermauw et al. 2020). Of these, the CYP3 group is the largest and most dynamic, exhibiting many lineage-specific duplications in insects, often linked to insecticide resistance and xenobiotic metabolism (Feyereisen 2006).

There appears to be a link between the function of P450s and their evolutionary dynamics. In vertebrates, genes involved in xenobiotic detoxification are more likely to be evolutionarily unstable, exhibiting duplications and deletions between species (Thomas 2007). In contrast, those involved in viability are more likely to be stable with a one-to-one orthology between species (Thomas 2007). This pattern has also been found in Drosophila where P450s associated with a role in development are duplicated less than those implicated in detoxification or have unknown functions (Drosophila 12 Genomes Consortium 2007).

Although the link between evolutionary instability in P450s and xenobiotic function has been shown, this does not prove whether selection is acting on these P450s. The expansion of P450 subfamilies involved in xenobiotic functions is often assumed to be due to environmental adaptation driven by natural selection. However, even if an initial duplication is selected for, further duplications that lead to subfamily expansions could simply be due to the self-sustaining process whereby a duplicated gene has twice the likelihood of duplicating again (Feyereisen 2011). In contrast to this neutral viewpoint, at a molecular level, correlation between amino acid replacement and the number of duplications found for a P450 lineage in the Drosophila phylogeny is inconsistent with stochastic models (Good et al. 2014). Furthermore, CYP expansions have been linked to environmental factors, such as specialized diets in Lepidoptera (Calla et al. 2017), implicating adaptive evolutionary forces in P450 expansions. However, the observed pattern of many groups with few genes and few groups with many genes (power-law distribution) does not require an adaptive explanation (Dermauw et al. 2020). Birth-and-death models of gene family evolution are sufficient to explain the pattern, not requiring any further explanation based on the ecology or life-history of the species (Sezutsu et al. 2013).

An interesting group for the study of P450s and detoxification is bees. While providing important ecosystem services as pollinators bees are exposed to a wide range of toxins, both natural and synthetic (Johnson 2015). In particular, bees differ from other insect pollinators in that they are specialized on consuming pollen and nectar during all life-stages, a diet with high levels of potentially toxic flavonoids, and in fact bees have been called “flavonoid specialists” due to their consistent exposure to these compounds (Johnson et al. 2018). Exposure to flavonoids is thought to be higher in perennial eusocial bees due to the concentration of flavonoids when nectar is converted into honey and pollen into beebread for storage (Johnson et al. 2018). Furthermore, resin collection, not only associated with the evolution of sociality but also found in some solitary species, also increases flavonoid exposure (Bankova et al. 1983).

Understanding the molecular basis of toxin sensitivity is important to protect bees, especially from the negative effects of pesticides (Berenbaum and Johnson 2015). The sequencing of the honeybee genome (Apis mellifera) revealed a lower diversity of detoxification genes than expected when compared with other insect species (Claudianos et al. 2006; The Honeybee Genome Sequencing Consortium 2006). Although the beetle Tribolium castaneum and the mosquito Anopheles gambiae have over 100 P450s, the honey bee A.mellifera has <50 (Ranson et al. 2002; Claudianos et al. 2006; Oakeshott et al. 2010; Zhu et al. 2013). This deficit of P450s is not only exclusive to A.mellifera, but is also the case in other bee species, including the solitary leafcutter bees Megachile rotundata and Osmia bicornis (Beadle et al. 2019; Hayward et al. 2019). Furthermore, this reduction is not equal across all P450s. Mitochondrial and CYP2 P450s have mostly been maintained, with the genes present from these groups linked to physiological roles in other insects, such as hormone biosynthesis (Claudianos et al. 2006). In contrast, the CYP4s are greatly reduced, especially those members thought to play a role in xenobiotic detoxification. The CYP3s are also greatly reduced in diversity, but widespread duplication in the CYP6AS family within the CYP3s means that these make up a large proportion of total bee P450S (Claudianos et al. 2006; Berenbaum and Johnson 2015). This may explain why, in spite of these reductions, honey bees are not more sensitive to insecticides that other insects (Hardstone and Scott 2010).

One particular group of bees, the orchid bees, exhibits unique natural history and adaptation that presents additional challenges for detoxification. Male bees collect chemical compounds from different sources, such as orchid flowers and fungi, which they store in specialized leg pouches to concoct a “perfume” bouquet (Dressler 1982). These bouquets are then released during courtship displays and are thought to play a key role in mate choice (Eltz et al. 2005; Pokorny et al. 2017). The chemical bouquets differ between species, but also between individuals of the same species and may provide information to the female regarding mate quality (Eltz et al. 1999; Zimmermann et al. 2009; Weber et al. 2016). Male orchid bees collect a wide range of chemical compounds, many of which are expected to be toxic to bees. Collection of these compounds, therefore, could act as an indicator of male quality, through their ability to handle and detoxify these compounds (Eltz et al. 1999; Arriaga-Osnaya et al. 2017). We may expect adaptation at the molecular level, perhaps in the form of an increased P450 repertoire, or increased expression levels, to be able to detoxify such a wide range of compounds.

Here, we study the molecular evolution of cytochrome P450s in bees and evaluate the hypothesis that orchid bees have an expanded P450 repertoire. We annotate P450s in the genome of Euglossa dilemma and combine these newly annotated P450s with those previously identified in nine other bee species: A.mellifera, Bombus terrestris, Dufourea novaeangliae, Eufriesea mexicana, Habropoda laboriosa, Lasioglossum albipes, M.rotundata, Melipona quadrifasciata, and O.bicornis. These species are from three different families, and include Ef.mexicana, another orchid bee. To determine the evolutionary history of these P450s we carry out phylogenetic analyses and classify clades as stable or unstable based on their history across these species. We measure evolutionary change and search for signals of adaptive evolution in these clades to identify the selection pressures acting on P450s in these bee species. We also investigate the patterns of correlation between the level of sociality of a bee species, whether it collects resin, and its P450 inventory. We aim to shed light both on the evolution of the P450 family and orchid bee biology in relation to detoxification.

Results

Phylogenetic Analysis of Bee P450s

To compile a data set of bee P450s we firstly combined previously annotated P450s for nine species: A.mellifera, B.terrestris, D.novaeangliae, Ef.mexicana, H.laboriosa, L.albipes, M.rotundata, M.quadrifasciata, and O.bicornis (Kapheim et al. 2015; Johnson et al. 2018; Beadle et al. 2019; Hayward et al. 2019) (fig. 1). To this set of P450s, we then added 41 P450s which we identified and annotated in Eg. dilemma with complete protein domains. We also identified four P450s with incomplete protein domains that may represent pseudogenes or may be the result of poor assembly, sequencing errors, or sequencing gaps. We do not include these for the following analyses (compiled in the file incomplete_genes.fa).

Fig. 1.

Phylogeny of cytochrome P450s in bees. The phylogeny was constructed in IQ-TREE (model GTR+F+R4) using five nuclear genes at 10,508 nucleotide sites. Bootstrap values (n = 1,000) are illustrated. For each species its resin-collecting behavior, level of sociality, and number of P450s in each family is shown. Melipona quadrifasciata photo courtesy of Lucas Rubio, O. bicornis photo courtesy of Amelia Bassiti.

To determine orthology of P450s we constructed a phylogeny combining the Eg. dilemma P450s with those previously identified in nine other bee species. The P450s group into the four expected clusters in insects, CYP2, CYP3, CYP4, and mitochondrial (fig. 2; supplementary fig. 1, Supplementary Material online). As expected for bees, the CYP4s are reduced across all species, and the CYP6AS family is expanded. In some clades, the phylogeny of P450s does not reflect the species’ phylogeny. This is likely due to errors in phylogenetic reconstruction in these particular P450 clades rather than complicated patterns of duplicated and deletion between species. No members of the CYP6AS clade containing CYP6AS10 from B. terrestris and A. mellifera are found in either Eg. dilemma or Ef. mexicana, suggesting this is a loss in orchid bees. Overall, Eg. dilemma has a comparable number and distribution of P450s as other previously studied bee species (fig. 1; supplementary fig. 1, Supplementary Material online).

Fig. 2.

Phylogeny of cytochrome P450s in bees. The phylogeny was constructed in IQ-TREE (model JTT+F+R10) using 481 amino acid sequences across ten bee species belonging to three different families. Clades are collapsed to show 40 clades. Bootstrap values (n = 1,000) are illustrated.

Is the Number of P450s in the Genome Associated with Species’ Ecology?

We found that total number of P450s in the genome of each species ranged from 40 in H. laboriosa to 59 in D. novaeangliae, with most of this variation driven by the number of CYP3s found in each species (fig. 1). We investigated two aspects of the ecology of the bee species included in our study: resin collection and level of sociality. We tested for relationships between P450 number and ecology while accounting for phylogeny using a phylogenetic independent contrast (PIC). We did not find a relationship between sociality and the number of P450s (supplementary fig. 2A, Supplementary Material online), the number of CYP3s (supplementary fig. 2B, Supplementary Material online), or the number of CYP6AS genes (supplementary fig. 2C, Supplementary Material online). Furthermore, we did not find a difference between bee species that exhibit substantial resin collection and those which do not in the number of P450s (supplementary fig. 3A, Supplementary Material online), the number of CYP3s (supplementary fig. 3B, Supplementary Material online), or the number of CYP6AS genes (supplementary fig. 3C, Supplementary Material online).

Evolutionary History of P450s in Bees

To investigate the dynamics of gene family evolution and phylogenetic instability we used two approaches: The reconciliation of the gene tree with the species tree to detect duplication and loss events; and the use of birth–death models fitted to a species tree (Yohe et al. 2019).

We used MiPhy to implement the first approach of species tree and gene tree reconciliation to identify clades and estimate phylogenetic instability scores for each clade. Without the “merge singletons” option, MiPhy divided our data into 42 clades including two clades which only contained a sequence from H. laboriosa. For the final analysis, we used the option “merge singletons” to divide our data into 40 clades ranging in their instability score from −0.1 to 25.05 (supplementary fig. 4, Supplementary Material online). The instability scores of the clades fell into two broad groups. Clades with low scores (−0.1 to 3.61) were classified as stable, and those with higher scores (9.12–25.05) were classified as unstable. In total, 27.5% (11/40) clades were classified as stable. The distribution of stable and unstable clades was not equal amongst P450 families, with 52% of CYP3 clades found to be unstable, in comparison with only 0% of CYP2, CYP4 and mitochondrial P450 clades (fig. 2). CYP3s have a higher instability score than the other CYP groups (supplementary fig. 5, Supplementary Material online).

To apply the second approach to our data set we use computational analysis of gene family evolution (CAFE), a program which uses birth–death to model gene gain and loss across a species tree, using the clades previously identified by the MiPhy analysis (Mendes et al. 2020). Clades which are fast evolving are identified by comparing models in which all clades evolve at the same rate with models in which different clades vary in their evolutionary rate. Using this approach, we identified one CYP6AS clade, CYP6AS1, as an outlier when compared with all other families, with respect to its evolutionary rate of family size change (P = 0.029). This was the same clade which had the highest instability in the MiPhy analysis (supplementary fig. 4, Supplementary Material online). As expected, this family also had the strongest assignment to the gamma category of the highest median lambda (rate of evolutionary change), (posterior probability = 1.0). Furthermore, although not detected as outliers relative to the entire data set, seven additional clades were also inferred to be rapidly evolving in size, as indicated by their significant placement in the same fast evolving gamma category (supplementary table 1, Supplementary Material online; posterior probability > 0.95). All of these were classed as unstable by the MiPhy analysis. Although all of the clades identified by CAFE were identified as unstable in the MiPhy analysis, this was not true in reverse, probably as CAFE uses gene counts and does not consider the gene tree, making it unable to distinguish between independent gene duplication events and inheritance of paralogs (Curran et al. 2018).

The exact pattern of duplications and deletions is unclear, but certainly some have happened on a more recent timescale. For example, CYP9DN1 is not found in A. mellifera, B. terrestris, D. novaeangliae, or Ef. mexicana, suggesting this deletion has occurred at least twice as D. novaeangliae is not closely related to the other species. Of the 11 unstable clades, three belong to this category of requiring at least two deletion events to explain the phylogenetic pattern. Other clades could be explained by one deeper internal deletion event, but more species would be needed to determine this. For example, the clade containing A. mellifera CYP6AS10 is missing in both Ef. mexicana and Eg. dilemma, which could be explained by a deletion in all orchid bees, or two separate deletion events, with further sampling necessary to distinguish the two. Duplication events are also a mix of lineage-specific and potentially deeper internal duplications. The genes CYP6AS131–135 and CYP6AS91–95 are expanded specifically in O.bicornis and D.novaeangliae, respectively, suggesting more recent duplication events. In contrast, CYP336 shows a pattern consistent with duplications deeper in the lineage leading to O.bicornis and M.rotundata (supplementary fig. 1, Supplementary Material online).

Do Unstable Clades of P450s Exhibit Increased Evolutionary Change?

To further investigate the evolutionary dynamics of these clades we compared various evolutionary measures. Branch length can be used as a measure of evolutionary change in a phylogenetic tree. We found a correlation between instability and branch length, when gene number is corrected for, suggesting higher rates of evolution than for genes found in stable clades (fig. 3). Unstable clades also have a higher cumulative patristic distance, which includes the length of internal branches also (supplementary fig. 6, Supplementary Material online). One clade, CYP9FX, did not follow this trend, with a low instability score but high CBL and cumulative patristic distance. This is the smallest clade in our analyses with only three genes from D. novaeangliae and L. albipes.

Fig. 3.

Relationship between CBL and the instability of a clade. When CBL is normalized for gene number in a clade, there is a significant correlation between clade CBL and clade instability (Spearman’s correlation, ρ = 0.401, P = 0.0103). CBL = cumulative branch length.

Are There Signals of Adaptive Evolution in Bee P450s?

Branch-Specific Models

Although branch length can give an indication of the amount of divergence between two P450s, it does not provide information on the predominant type of selection that shaped their evolution. To explore both the amount and the type of evolutionary change happening in stable and unstable clades, we compared the rate of nonsynonymous with synonymous mutations (dN/dS) in different lineages. The simplest analysis carried out was a “one-ratio” model for each clade. In these models, each branch in a clade is assumed to have the same dN/dS value. All values obtained were below one when all members of a clade are forced to have the same dN/dS value, suggesting that in general purifying selection is acting. To test if unstable and stable clades are under different selection pressures, we compared their “one-ratio” dN/dS values. We found that all clades are under purifying selection with dN/dS values less than one. However, we detected a correlation between clade instability and dN/dS ratio, with more unstable clades having a higher dN/dS ratio (fig. 4). Furthermore, we found a correlation between the cumulative branch length (CBL) of a clade, the amount of evolutionary change, and the dN/dS value of the clade (fig. 4).

Fig. 4.

Relationship between the dN/dS of a clade, clade instability, and normalized CBL. (A) Correlation between the instability of a clade and its dN/dS ratio in “one-ratio” models (Spearman’s correlation, ρ = 0.567, P = 0.00014). (B) Correlation between Normalized CBL (average branch length per gene in a clade) and dN/dS of the clade in a “one-ratio” model (Spearman’s correlation, ρ = 0.754, P = 1.889e−08). CBL = cumulative branch length.

The “one-ratio” model is restrictive, only allowing a single dN/dS value per clade. To investigate which clades are under more dynamic selection pressures than the simple “one-ratio” model, we compared these models with models under which each lineage can vary in its dN/dS ratio, also called “free-ratio” models (supplementary table 2, Supplementary Material online). As for the “one-ratio” models above, these “free-ratio” models were carried out by modeling each clade individually. The evolutionary history of unstable clades is more dynamic than stable clades. The more complex model has a better fit for more of the unstable clades (100%, 11/11) than the stable clades (69%, 20/29) (Fisher’s exact test, P = 0.0433).

As well as testing selection pressures at a clade level we tested each branch in each clade individually for evidence of selection. We compared “one-ratio” models for each clade with models in which only one branch is allowed to vary in its dN/dS ratio (“two-ratio” models). In these models, we can ask whether the branch of interest has a dN/dS value significantly different from the “background” dN/dS of the clade. Again, we found the same pattern in the resulting dN/dS values. Of all the branches with a dN/dS ratio significantly different from the background of the clade, dN/dS was higher for branches from unstable than stable clades (supplementary fig. 7, Supplementary Material online). Furthermore, all branches with a significant dN/dS >1.1 (excluding branches where dS < 0.01 and dS > 2), indicating positive selection, were from the CYP3 clade. Of these, 70% (7/10) were from unstable CYP3 clades, and 70% (7/10) were from the CYP6AS family (table 1). The three genes which were not found in unstable clades were all in the CYP6AS7 clade. Full list of tests can be found in supplementary data 2, Supplementary Material online.

Table 1

Branches with dN/dS >1 and Significantly Different from Background dN/dS of the Clade in the “Two-Ratio” Tests (Adjusted P Values <0.05)

Name	dN/dS	CYP Clade	CYP Family	Genus
GB40288	2.316	CYP3	CYP6AS1	Apis mellifera
LALB19295RA_B	1.822	CYP3	CYP6AS200	Lasioglossum albipes
MQUA24239RA	1.597	CYP3	CYP336A42	Melipona quadrifasciata
Internal branch	1.597	CYP3	CYP336
MH500637	1.322	CYP3	CYP6AS127	Osmia bicornis
NW_015373885.1C	1.268	CYP3	CYP9FW6	Dufourea novaeangliae
Internal branch	1.243	CYP3	CYP6AS7
MH500655	1.224	CYP3	CYP6AS131	Osmia bicornis
MH500653	1.206	CYP3	CYP6AS132	Osmia bicornis
XM_017902249.1	1.137	CYP3	CYP6AS87	Eufriesea mexicana

Note.—Branches with dS < 0.01 and dS > 2 have been removed.

Site-Specific Models

Branch-specific models assume a consistent dN/dS ratio across all sites in a given gene. We also carried out site-specific tests which allow dN/dS to vary among sites. We found evidence of positively selected sites in 17.5% (7/40) of all clades in the M7 and M8 model, including clades from CYP4 and CYP3 (supplementary table 3 and data 3, Supplementary Material online). For four of these clades, CYP6AS1, CYP6AS7, CYP6AS8, and CYP336, positively selected branches were previously identified in the branch-specific test comparing “one-ratio” and “two-ratio” models. CYP6AS1 and CYP336 are the two clades with highest instability as identified by MiPhy and CYP6AS1 was also identified as the sole outlier family in terms of evolutionary rate by CAFE. We identified positively selected sites at 90% posterior probability for three of these clades: CYP336, CYP6AS7, and CYP6AS8. We found no evidence of positively selected sites in the M1a–M2a comparison.

To determine where on the protein these sites were located, we took advantage of the AlphaFold2 algorithm (Jumper et al. 2021; Mirdita et al. 2021) to create protein structural models for a representative sequence in each clade. In general, the residues are found on surface regions of the proteins. The resides identified in CYP6AS7 and CYP6AS8 are found in the same surface region of the protein (supplementary fig. 8, Supplementary Material online). The residue in CYP6AS7 is found between the F′ and F helices, and the residue for CYP6AS8 is found between the G and G′ helices. This region is not close to the active site and may be involved in protein–protein interactions (Sevrioukova et al. 1999). The residue for CYP336 is found in another region of the structure (supplementary fig. 8, Supplementary Material online). The residue in CYP336 is located on a nonhelical, nonbeta surface area of the protein that is not part of the active site of the protein.

Discussion

Although the role of duplication and deletion in the evolution of multigene families is well-established, controversy remains surrounding whether this evolution is primarily adaptive or neutral. In this study, we demonstrate that genes known to be involved in xenobiotic metabolism (CYP3 family) are more likely to be unstable in the form of more gene duplication and deletion. Furthermore, we find that not only are these genes more likely to be unstable but also under more dynamic evolutionary pressures, and exhibit signals of adaptive evolution. This suggests that both gene duplication and positive selection driving sequence divergence contribute to the diversification of P450s. We do not find evidence for a correlation between the P450 repertoire of a bee species and its biology. Furthermore, our hypothesis that orchid bees would have an expansion of P450s due to their perfume collection behavior was also not supported.

The patterns we identify here are very similar to those previously identified in Drosophila. P450s with developmental functions are more likely to be found in stable clades, and also more likely to be under purifying selection (Drosophila 12 Genomes Consortium 2007; Good et al. 2014). These patterns, however, are not restricted to P450s, and have also been identified in the fatty acyl-coA synthase multigene family in Drosophila, which are involved in both essential physiological and chemical communication functions (Finet et al. 2019). While no single evolutionary model will be able to explain multigene family evolution, we can understand some common underlying principles using the framework of the birth-and-death model (Eirín-López et al. 2012).

Another common pattern found across multigene families is lineage-specific expansions, which are often associated with the evolution of novelty (Lespinet et al. 2002). In the CYP3 clade of P450s in bees, we find extensive evidence of duplications, particularly in the CYP6AS subfamily. It has been suggested that this expansion is related to specialization on a plant-based diet and the ability to metabolize flavonoids, in comparison to solitary carnivorous ancestors such as solitary wasps (Johnson et al. 2018). As previously described by Johnson et al. (2018), we find that not only are CYP6AS subfamilies expanded by duplication, but also show evidence of positive selection, with the majority of positively selected genes detected belonging to this family. CYP6AS3 and CYP6AS4 from A. mellifera were found to be under positive selection by Johnson et al. Here we find CYP6AS1 to be under positive selection in the branch-specific models and evidence for positively selected sites in the clade containing CYP6AS1 and CYP6AS3 in the site-specific models. CYP6AS1, CYP6AS3, CYP6AS4, and CYP6AS10, have all been linked to quercetin metabolism, a plant flavonoid found in honey (Mao et al. 2009). Other CYP6AS genes, as well as CYP9R1, CYP9P1, CYP9S1, and CYP9Q family genes are upregulated in response to quercetin treatment in A. mellifera (Mao et al. 2017). Furthermore, CYP9Q4 in B. terrestris and CYP9Q3 in A. mellifera metabolize the neonicotinoid thiacloprid (Manjon et al. 2018). All of these genes known to be involved in xenobiotic metabolism were identified in our study as members of unstable clades, again linking instability and gene duplication with xenobiotic function.

In contrast to these genes with known xenobiotic function, we find that genes predicted to have physiological roles are mostly found in stable clades under purifying selection. The mitochondrial clade genes CYP302A1, −314A1, AND −315A1, which are orthologs of the Drosophila melanogaster Halloween genes, are thought to be involved in ecdysteroid synthesis and are all found in stable clades in our analysis (Gilbert 2004; Rewitz et al. 2007; Feyereisen 2011). From the CYP2 clade, CYP18A1, thought to be involved in ecdysteroid inactivation, CYP15A1, a predicted juvenile hormone epoxidase, and CYP306A1, a CYP2 member orthologous to a D. melanogaster ecdysteroid 25-hydroxylase, are also found in stable clades (Helvig et al. 2004; Niwa et al. 2004; Claudianos et al. 2006). However, CYP6AS7 and CYP6BD1, despite being stable clades, showed evidence of positive selection in site-specific models. Furthermore, the dichotomy of physiological and xenobiotic genes is not so simple, as both can be found in phylogenetic proximity (Dermauw et al. 2020) whereby members of the same subfamily may have divergent roles. For example, while the CYP6AS subfamily is known to play a role in xenobiotic detoxification, CYP6AS8 and CYP6AS11 of A. mellifera are expressed in mandibular glands, potentially involved in fatty acid signal synthesis (Wu et al. 2017; Dermauw et al. 2020). This highlights the complexity of multigene family evolution which does not always follow generalized predictions.

If P450 function is linked to their evolutionary dynamics, we might expect groups or species with different biology to differ in their P450s (Rane et al. 2019). For example, lineage-specific expansions of insect metallothioneins in Lepidoptera, important for heavy metal detoxification, is thought to be linked to the wide range of host species/tissues used throughout their lifetime (Luo et al. 2020). In this study, we hypothesized that Eg. dilemma would exhibit an expanded repertoire of P450s when compared with other bees, due to its perfume-collecting behavior increasing exposure to different chemical compounds. In contrast, we found that Eg. dilemma has a comparable set of P450s similar to other bees, with no significant expansions detected. This does not rule out an increased capacity for xenobiotic metabolism, which would need to be tested experimentally. Specialization could also occur at the sequence level rather than the number of genes; however, we did not find evidence for positive selection in Eg. dilemma P450s. One alternative is that Eg. dilemma exhibits increased gene diversity at another step in the detoxification pathway. In honeybees, not only are P450s reduced in comparison to other insects, but also glutathione-S-transferases (GSTs) and carboxyl/cholinesterases (CCEs) (Claudianos et al. 2006; Berenbaum and Johnson 2015). Both GSTs and CCEs are known to be involved in xenobiotic metabolism in insects, thus meriting further investigation in Eg. dilemma (Li et al. 2007). Another possibility is that orchid bees do not demonstrate higher abilities to metabolize xenobiotic compounds, but instead reduce the amount of compound that can enter the body, for example, by decreasing cuticle penetration, as reported in A.gambiae (Balabanidou et al. 2016).

We also might expect a reduction in P450 repertoire in species which have a potentially lower xenobiotic exposure due to a specialist lifestyle (Rane et al. 2019). Drosophila sechellia, a specialist island species only found within a narrow ecological niche exhibits an increased loss of P450 genes compared with other Drosophila species (Good et al. 2014). In our data set, the two most specialized species are D.novaeangliae, a specialist, with only one known pollen source (Eickwort et al. 1986), and H. laboriosa, with only a few known pollen sources (Pascarella 2007). Interestingly, D. novaeangliae has the largest P450 repertoire of all species in our analysis and H. laboriosa has the fewest P450s. We therefore do not find a strong link between specialization in diet and P450 repertoire in our bee data set. More generally, we did not find a relationship between P450s and the biology of different bee species, neither at the total number level nor subfamily level. This is in contrast to the previously identified trend that degree of sociality correlates with CYP6AS subfamily size, and that resin-collecting bees have more CYP6AS members (Johnson et al. 2018). This is likely because of the inclusion of O. bicornis in this study, a solitary bee which does not collect resin and yet has a large P450 repertoire, and the inclusion of phylogeny in these analyses. Future studies with more species, in particular species such as specialists, will allow us to untangle the relationship between species’ biology and P450s.

The P450 family in bees can provide us with insights into the evolution of multigene families. In addition to this evolutionary perspective, understanding detoxification in bees is crucial to native bee conservation and the management of agricultural pollination in the face of widespread pesticide use. Bees provide an important ecosystem service through pollination, both of wild plants and economically important crop species (Klein et al. 2007; Vanbergen and Initiative 2013). Although most insecticide assessments are carried out in A. mellifera, they do not always reflect the sensitivity of native bees, with variation between species (Arena and Sgolastra 2014; Franklin and Raine 2019). Toxicity assays require large numbers of bees, difficult to achieve with many native bee species. Understanding the genetic mechanisms underlying insecticide sensitivity may allow us to make predictions based on P450 repertoires and sequences in different species (López-Osorio and Wurm 2020). A recent example of this comes from M.rotundata, a species which exhibits increased sensitivity to neonicotinoids, and also lacks the P450 enzymes known to be involved in neonicotinoid metabolism in other species (Hayward et al. 2019). A combination of comparative genomic studies, alongside toxicity assays and molecular and functional studies of metabolism and insecticide resistance will allow us to extend our understanding to other bee species.

Materials and Methods

Gene Family Annotation

First, we gathered previously annotated P450s for nine species: Apis mellifera, B.terrestris, D.novaeangliae, Ef. mexicana, H.laboriosa, L.albipes, M.rotundata, M.quadrifasciata, and O.bicornis (Kapheim et al. 2015; Johnson et al. 2018; Beadle et al. 2019; Hayward et al. 2019). We combined all P450s described in previous publications and aligned them using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013). We then used the bio3d package in R to compare the sequences, removing duplicates to curate a final data set (Grant et al. 2021), We searched either the National Center for Biotechnology Information (NCBI) or the Hymenoptera Genome Database to provide accession numbers for the included P450s. (Elsik et al. 2018; NCBI Resource Coordinators 2018). For genes not currently included in an annotation but previously manually annotated, we included the scaffold accession number where the sequence can be found (supplementary data 1, Supplementary Material online). We checked each sequence for a complete protein domain using the NCBI conserved domain search (Marchler-Bauer et al. 2015). In some cases, the sequence from the database contained an incomplete domain which had been manually edited in the data set provided by Johnson et al. (2018), and in these cases we used the manually curated sequence and noted this (supplementary data 1, Supplementary Material online).

We then annotated the P450 gene family in the Eg. dilemma genome (Brand et al. 2017) using all previously annotated A. mellifera protein sequences as a reference (Elsik et al. 2018). Firstly, we identified scaffolds in the Eg. dilemma genome containing P450s using TBLASTN (Altschul et al. 1990), and then annotated proteins on these scaffolds using exonerate (Slater and Birney 2005). Gene models were extracted, translated to amino acid sequences, and checked by comparison to A. mellifera proteins. Gene models were manually curated using IGV, checking for splice sites and start and stop codons (Robinson et al. 2011). To check for P450s missed by the initial search, the predicted Eg. dilemma P450s, as well as previously described Ef. mexicana P450s, were used as a query for the initial step of TBLASTN followed by exonerate predictions (Altschul et al. 1990; Slater and Birney 2005). As a further check we used already published RNA-seq data of ovaries and brains of Eg. dilemma (Bioproject PRJNA523381) to improve the annotations (Saleh and Ramírez 2019). We trimmed the reads using TrimGalore! (Martin 2011). We then mapped the genes (two-pass alignment) to the Eg. dilemma genome (Brand et al. 2017) with the manually edited annotation (Edil_v1.0_revised.gff) using STAR (Dobin et al. 2013). We used StringTie (Pertea et al. 2016) to assemble and merge the transcripts to form an annotation which we compared with the edited annotation in IGV (Robinson et al. 2011). To check that the predicted P450s contain complete functional protein domains we used the NCBI conserved domain search (Marchler-Bauer et al. 2015). One gene, Edil_14289, is found at the end of a scaffold in a poorly sequenced area. The exon from the unsequenced area is found at the end of another scaffold and we combined these to form a complete gene in our analyses. All P450s identified were included in a revised version of the current Eg. dilemma annotation (Edil_v1.0_revised.gff) and included in supplementary data 1, Supplementary Material online.

Phylogenetic Analyses

All 481 bee P450 sequences were combined with CYP51 from Mus musculus as an outgroup. The sequences were aligned using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013). Phylogenetic trees were constructed using IQ-TREE with the ModelFinder function to determine the best-fit model (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018). The Newick phylogeny was plotted using MEGA X, and the packages ape, evobiR, and geiger in R version 3.6.2 (Pennell et al. 2014; Blackmon and Adams 2015; Kumar et al. 2018; Paradis and Schliep 2019; R Core Team 2019; Stecher et al. 2020).

To test for correlation between P450 repertoire and bee biology we carried out phylogenetic comparative analyses. The data set used to produce the species-level phylogeny was composed of five nuclear genes: wingless, arginine kinase, opsin (long-wavelength), Nak, Ef1α (F2 copy) (Danforth et al. 2011). We identified copies of the five genes in Apis mellifera from NCBI (NCBI Resource Coordinators 2018). We then used TBLASTN to search the nucleotide data sets on NCBI for the sequences in B.terrestris, D.novaeangliae, Ef. mexicana, H.laboriosa, M.rotundata, and O.bicornis (Altschul et al. 1990). For Eg.dilemma, L.albipes, and M.quadrifasciata, we used TBLASTN and exonerate to find and annotate the proteins in the genome (Altschul et al. 1990; Slater and Birney 2005; Kocher et al. 2013; Kapheim et al. 2015; Brand et al. 2017). The sequences were aligned using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013). Alignments were trimmed using trimAl (-gt 0.6, sites only included where sequence from 6/10 bee species present) (Capella-Gutierrez et al. 2009), and trimmed alignments were then concatenated. We constructed phylogenetic trees constructed using IQ-TREE with the ModelFinder function to determine the best-fit model (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018), and constrained the relationship of corbiculate bees according to a previous phylogeny (Romiguier et al. 2016). We rooted the tree using the midpoint.root function in the phytools package in R (Revell 2012).

Bees were determined as resin collecting or nonresin collecting following Johnson et al. (2018). We designated level of sociality as described previously (1 = ancestrally solitary, 2 = facultative basic eusociality, 3 = obligate basic sociality, and 4 = obligate complex eusociality) (Kapheim et al. 2015). We tested for correlation between number of P450s, CYP3s, or CYP6AS genes and the level of sociality by using Spearman’s rank correlation on PICs. We calculated PICs using the function pic in the ape package (Paradis and Schliep 2019). We tested for a relationship between resin collection and P450 repertoire, using phylogenetic ANOVA, with the function aov.phylo from the package geiger (Pennell et al. 2014).

Clade Stability and Branch Length Analyses

We used MiPhy to identify clades and assess their stability (Curran et al. 2018). MiPhy requires a rooted tree as an input, however, our root did not belong to any of the species present in the species tree. Therefore, we removed the Mus musculus CYP15 from the data set, aligned the sequences using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013), and constructed a phylogenetic tree using IQ-TREE with the ModelFinder function to determine the best-fit model (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018). We rooted the tree using the midpoint.root function in the phytools package in R (Revell 2012). We used default parameter settings along with the option “merge singletons” to ensure that all clades had more than one gene. Clades were named using nomenclature from A. mellifera. For the CYP6AS subfamily which contained multiple clades, we named the clades using the lowest numbered group member from A. mellifera.

To compare instability between the four CYP groups, we used a Welch one-way ANOVA test, followed by a post hoc Games–Howell test using the package rstatix (Kassambara 2021). These tests do not assume homogeneity of variance.

We then applied CAFE to our data set to analyze the evolutionary rates of different clades. This program uses a birth–death process to model gene gain and loss across a species tree to identify fast-evolving clades. The input files required are a species tree and gene counts for each species in each clade, in this case the clades identified by MiPhy. We used the python script provided in the CAFE tutorial with the program r8s to convert the species tree constructed as described above into an ultrametric tree, as required for CAFE input (Sanderson 2003). We calibrated the tree using a previously published time-calibrated phylogeny (Cardinal et al. 2018). We found that three gamma rate categories (k = 3) were the best fit for our data and showed convergence between runs.

For each clade we determined the CBL by adding the terminal branches for each gene within a clade, obtained from the IQ-TREE output. As each clade has a different number of genes, we normalized the CBL by dividing the CBL by the number of genes in the clade, following (Finet et al. 2019). We also calculated the cumulative patristic distance per clade by summing all branch lengths, including internal branches within the clade. Again, we normalized the cumulative patristic distance by dividing by the number of genes in the clade. We compared both CBL and cumulative patristic distance between stable and unstable clades using a t-test. We tested for correlation between both CBL and cumulative patristic distance and the instability of a clade using Spearman’s rank correlation.

Selection Analyses

We split the data set into the previously identified clades and carried out analysis individually on each clade. To construct a codon alignment for selection analyses we first aligned amino acid sequences using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013). We then used PAL2NAL to align the corresponding nucleotide sequences to the amino acid alignment (Suyama et al. 2006). Phylogenetic trees were constructed using IQ-TREE (Nguyen et al. 2015; Kalyaanamoorthy et al. 2017; Hoang et al. 2018).

Gene Conversion

As gene conversion can lead to false positive results when testing for positive selection with phylogenetic analysis by maximum likelihood (PAML), we first tested our data set for evidence of gene conversion (Casola and Hahn 2009). We split the data set by species and aligned the sequences from each species separately using MAFFT v7.453 (-maxiterate 1000, using L-INS-I algorithm) (Katoh et al. 2002, 2005; Katoh and Standley 2013). We tested for gene conversion using GENECONV in RDP5 (Padidam et al. 1999; Martin et al. 2020). We tested the sequences as linear sequences and found no evidence for gene conversion in our data set.

To compare selection pressure on clades which have undergone expansions with those which have not we used codon substitution models implemented in PAML (Yang 2007). Firstly, we carried out “one-ratio” models. For these, we assumed that dN/dS (the ω ratio of nonsynonymous to synonymous substitutions) has one value across the whole clade. We compared the dN/dS value between stable and unstable clades using a t-test, and tested for a correlation between the dN/dS value of a clade and its normalized CBL using Spearman’s rank correlation.

We then carried out “free-ratio” models, in which we allowed ω to vary between lineages. We compared the two models (“one-ratio” and “free-ratio”) for each clade using a likelihood ratio test (LRT). To correct for multiple-testing, we used the p.adjust function in R, with false detection rate (fdr) correction, which controls for the proportion of false positives. To test if the unstable clades were more likely to be explained by a “free-ratio” model then stable clades we used Fisher’s exact test.

Furthermore, to compare the selection pressures on individual genes or groups of genes within the clades, we again compared two models. The first assumed a single ω for the clade (“one-ratio”, as above), the second assumes two ω ratios, one for the lineage of interest, or group of interest, and the second for the rest of the tree. Again, we compared the two models using a LRT, correcting for multiple-testing. To automate the process of testing each branch we used ETE3 framework to implement codeml (–leaves—internals –codeml_param CodonFreq, 3 Nssites, 0 fix_omega, 0 omega, 1 fix_kappa, 0 kappa, 0 cleandata, 1 fix_blength, 1) (Huerta-Cepas et al. 2016).

To test for positively selected sites, we carried out two model comparisons: M1a versus M2a, and M7 versus M8. Model M1a is a nearly neutral model with two classes of sites, one evolving neutrally and the other under purifying selection. M2a is the same as M1a but with an additional site class for positive selection. Model M7 allows for a beta distribution of dN/dS across sites, while M8 has a beta distribution plus an additional site class with dN/dS > 1 for positive selection. As above, we compared the two sets of models using LRTs and corrected for multiple-testing. We identified sites under positive selection using the Bayes Empirical Bayesian (BEB) in the PAML output (Yang et al. 2005). We considered sites to be under positive selection when the BEB posterior probability was >0.9. The site-specific models are particularly sensitive to alignment issues and poor alignment can result in false positives. We repeated the analysis using the aligner PRANK, which uses phylogenetic information to distinguish between gaps caused by deletions and those by insertions, producing good alignments for evolutionary purposes, resulting in a lower false positive rate than other commonly used aligners (Markova-Raina and Petrov 2011; Löytynoja 2014). We combined these PRANK results with our MAFFT-alignment results and considered a clade or site to be under positive selection only when identified by both analyses.

To determine where on the proteins these positively selected residues are located, we created protein models using ColabFold (Mirdita et al. 2021). ColabFold predicts protein structure based on AlphaFold2 (Jumper et al. 2021). We used the first sequence from each clade as a representative. We then visualized the resulting models using PyMOL (DeLano 2002).

We made figures using the packages ggplot2 and cowplot in R version 3.6.2 (R Core Team 2019; Wickham 2009, p. 2).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.