Properties of “Stable” Mosquito Cytochrome P450 Enzymes

1. Introduction

Mosquitos are infectious disease vectors having a profound impact on human health. Anopheles mosquitoes transmit the parasite Plasmodium falciparum, the causative agent of malaria, Aedes aegypti is the vector of transmission of yellow fever, and Culex mosquitoes transmit West Nile virus and Saint Louis encephalitis virus, as well as the nematode that causes lymphatic filariasis. In the absence of vaccines for these diseases, insecticides are widely applied as a means of mosquito control. This strategy has limitations in that insecticides are not only highly ecotoxic but are also becoming increasingly ineffective due to acquired resistance of the targeted insects. Resistance is mediated not only by mutations in the target insect receptors but also through insecticide detoxification. Cytochrome P450 enzymes have been shown to be essential for insect development and survival but, probably a larger number of P450s are involved in insecticide metabolism and resistance.

Cytochrome P450 enzymes, also known as CYPs, are haem-containing monooxygenases found in all kingdoms of life [1]. They first appeared in the scientific literature in 1958 [2].The Pfam database of protein families and domains classifies CYPs as a distinct family of enzymes (Pfam entry PF00067) involved in the oxidative degradation of various compounds [3]. In P450s, the haem iron forms a pentacoordinate system with the axial sulphur of a conserved cysteine (Cys) protein amino acid residue in the active site of the enzyme [4]. The conserved Cys residue exists as a thiolate anion giving a characteristic absorption band of 450 nm which is also typical of the spectral properties of other cytochromes such as the b-type haem-proteins haemoglobin and myoglobin. However, unlike these cytochromes, CYP monooxygenases are not electron-transfer proteins and function in the presence of redox partners, usually NAD(P)H-dependent ferredoxin or FAD- and FMN- containing CPR-type reductases by transferring molecular oxygen to -CH, -NH and -SH moieties of substrates with concomitant reduction of the oxygen atom to water [5]. Consequently, a more appropriate term for CYP enzymes would be P450 ‘haem-thiolate’ or P450 monooxygenase proteins rather than cytochrome P450s.

CYPs show an extraordinary diversity in their reaction chemistry. In mammals, they participate in oxidative, peroxidative and reductive metabolism of numerous endogenous compounds such as fatty acids, cholesterol, steroids, retinoids, vitamin derivatives, bile acids, porphyrins, thromboxane A2, prostacyclins, eicosanoids and other lipid mediators. In plants, they are involved in the secondary metabolism of, amongst others, phenolic compounds, alkaloids, gibberellins. The diversity of biosynthetic reactions catalysed by P450 is reviewed by Fujiyama, et al. [6]. Last but not least, CYPs metabolise an enormous range of xenobiotics and endobiotics including drugs, insecticides, environmental chemicals and pollutants, as well as natural plant products, and bacterial metabolites. P450s are involved in the metabolic detoxification of xenobiotics (phase I metabolism). This involves the addition or unmasking of polar groups, such as such hydroxy, amine or sulphydryl groups, in the xenobiotic substrate followed by hydrolysis, oxidation or reduction. The resultant reaction intermediates are then further metabolised by phase II enzymes, mainly, glutathione-S-transferases [7], and transported into the extracellular space through interactions with transmembrane proteins (phase IIII), mainly ATP-Binding Cassette (ABC) Transporters [8].

By 2018, more than 300,000 CYP sequences had been mined and collected in all areas of the tree of life [9] but the precise function of individual CYP proteins remains largely unknown. It is estimated that only 0.2% of the genes deposited in different databanks have been functionally characterised [10,11]. CYPs were first classified into distinct groups by Nerbert and Gonzalez in 1987. Following phylogenetic criteria, gene organisation and sequence similarity, CYPs are grouped into kingdom-specific clusters, named clans, each of which represent genes that diverged from a single common ancestor and can include one or multiple families [12,13]. Generally, groups of proteins having amino acid identity over 40% are assigned to the same family, whereas proteins having identity above 55% are allocated to the same subfamily. Gene families that repeatedly cluster in the same phylogenetic clade are grouped into the same clan. In this classification scheme a family-specific number is given after the root symbol CYP, followed by a letter and a number indicative of the subfamily and the gene, respectively. Clans can consist of one or multiple families and are given the name of the smallest family number present in the clade. For example, if a single clade in phylogenetic analysis comprises the CYP7, CYP8 and CYP39 families, these families become part of clan 7 [14]. Clan-based classification can be volatile and the addition of new, more distant sequences can lead to interleaved branches obscuring the boundaries between different classes [15].

Amongst animals, insects constitute the taxon with most P450 genes, which are grouped into six clans, namely, CYP2, CYP3, CYP4, CYP16, CYP20 and mitochondrial (mito). Two clans, CYP16 and CYP20, are restricted to certain species in Apterygota and Paleoptera [11]. Mosquitoes, second only to the deer tick (Ixodes scapularis), have a disproportionately large number of P450 genes [16,17]. CYP protein coding genes in the diptera D. melanogaster (Dmel), A. gambiae (Agam), A. aegypti (Aaeg) and C. quinquefasciatus (Cqui) form groups in Clans 2, 4, 6 and the mitochondrial clan. The number of CYPs populating each clan vary considerably (Table 1).

Gene assignment to clans and families may change with each new genome release and be further complicated by the fact that genomes can consist of different numbers of CYPs found in natural populations [20]. CYP assignment discrepancies may also arise from genome annotations differing in the coverage and quality of assembly of sequenced genomes, or due to homozygosity and polymorphisms in gene copy numbers [17,21].

Genes of the CYP2 and mitochondrial clans are relatively conserved and they form many families with few or even single members. Members of these clans have been shown to participate in core developmental and physiological functions and their evolution is well conserved with many families having few or even single members [15,22]. CYP clans 3 and 4 are by far the most populated due to lineage-specific gene amplification of paralogs or “phylogenetic blooms” [20,23]. These “blooms” are thought to result in response to environmental stimuli and, in particular, to the selection pressure exerted by polluted habitats and, in particular, due to the widespread application of chemical insecticides [18,24].

Phylogenetic studies have shown that the expansion of the P450 gene repertoire is shaped by gene duplication, gene birth and death, and gene neo-functionalisation [20]. These events are driven by selection acting on what is available at the time and can be described by a power law typical of simple birth and death models. An implication of this model is that CYPs may switch their function from physiology to detoxification and vice versa [16] or that some CYPs genes may be functionally redundant [25].

Yet some CYPs may deviate from this stochastic evolutionary model and a small number of CYPs have been shown to be “stable” in evolutionary terms, as for example in D. melanogaster, where 31 such “stable” genes have been identified [26]. Amongst these are the so-called “Halloween” genes, such as spook (spo), phantom (phm), disembodied (dib), shadow (sad), spookier (spok) and spookiest. These genes are essential for the biosynthesis of moulting hormones [27]. Functional investigation in different tissues of D. melanogaster, using RNAi screens, revealed that knockdown of 9 of these CYPs, resulted in lethality [28]. In the light of this evidence, it has been proposed that natural selection may also act on copy number polymorphisms within a species. Although this “selectionist” model has been contested, it is agreed that “conserved orthologous genes” are present in some phylogenetic branches of CYPs [11].

Regardless of evolutionary model, it is still not possible to distinguish phylogenetic branches related to essential physiological and developmental functions from those related to detoxification. Even in the case of D. melanogaster which serves as a model organism, well over 60% of the total number of P450s are “orphan” in terms of functional characterisation [29]. Likewise, phylogenetic branches of CYPs cannot be assigned ubiquitous xenobiotic detoxification functions as (i) different CYPs respond to insecticide selection in different mosquito populations and geographical regions [30], and (ii) CYPs are also known to play an indirect role in insecticide resistance [31,32]. As a result, susceptible and resistant strains with a common genetic background are usually not available to compare, and susceptible reference strains are biased by geographic variation or genetic drift. This makes attempts to identify one or a few metabolic genes associated with insecticide adaptation difficult or even futile [33,34].

“Gene essentiality” refers to those genes where a single gene–knockout results in lethality or severe loss of fitness. Earlier work has established that essential and non-essential genes can be recognised on the basis of their respective physical, chemical and evolutionary properties [35,36,37,38]. Based on this, bioinformatic methods were used to differentiate “stable” from “labile” mosquito P450s. Given that “stable” CYP genes may also be involved in auxiliary, non-essential, metabolic networks the analogy with “gene essentiality” may be tentative.

2. Materials, Methods and Datasets

CYPs identified in the sources specified in Table 1 were submitted to VectorBase [39] to source genomic and protein data and CYPs with ambiguous annotations were removed. The total number of CYPs in the dataset was 486 (SI1.xlsx). These were mapped onto 83 ortholog groups derived from OrthoMCL [40,41]. The latter is a genome-scale algorithm for the identification of orthologous protein sequences and provides not only groups shared by two or more species/genomes, but also groups representing species-specific gene expansion families [40,41]. The dataset was further split into two subsets, one of which included CYPs with zero or up to three paralogs and one that included CYPs with four or more paralogs (Supplementary Information 1). The division is based on earlier work with D. melanogaster CYPs in which genes with up to three duplications were shown to derive from evolutionarily stable clades, and genes with more than three duplications from “dynamic/unstable” clades. In this paper, the term “labile” is used in preference to “dynamic/unstable” [26].

Redundancy reduction. Each subset was purged by removing genes/proteins having sequence similarity > 70%. Redundancy reduction is necessary to avoid undesirable bias in statistical analysis arising from the presence of multiple copies of similar genes/proteins [42]. SkipRedundant of the EMBOSS suite of bioinformatics tools [43] was used to obtain non-redundant (culled) datasets, one from the “stable” and “labile” subsets. The numbers of CYPs in each of the resulting subsets is summarised in Table 2.

Gene and protein sequence-based features. The gene, protein, functional features and bioinformatic methods used to characterise “stable” and “labile” P450 are summarised in Table 3.

Protein sequence properties. The length of our “stable” and “labile” proteins was obtained from VectorBase. The EMBOSS program Pepstats (https://www.ebi.ac.uk/jdispatcher/seqstats/emboss_pepstats) was used to obtain statistics of a number of properties of FASTA formatted protein sequences. These attributes include: molecular weight, number of residues, charge, isoelectric point, and amino acid composition. Hydrophobicity was determined by multiplying the frequency of each amino acid in each P450 with its Kyte and Doolitttle hydropathicity index [47] and summing up the resulting values.

Pepstats groups amino acids into nine categories: Tiny (A, C, G, S and T); Small (A, B, C, D, G, N, P, S, T and V); Aliphatic (I, L and V); Aromatic (F, H, W and Y); Non-polar (A, C, F, G, I, L, M, P, V, W and Y); Polar (D, E, H, K, N, Q, R, S, T and Z); Charged (B, D, E, H, K, R and Z); Basic (H, K and R) and Acidic (B, D, E, Z). Pepstats was run with the default parameters setting. MATLAB scripts were used to extract features values from the output files and determine corresponding statistical values (MATLAB. (2024). Version R2024a. Natick, Massachusetts: The MathWorks Inc.).

Genomic properties. Genomic data (gene, CDS, intron/exon lengths) were obtained from VectorBase. %CG content data were obtained from the Ensembl Metazoa release 60 (http://metazoa.ensembl.org/index.html). For genes with multiple transcripts, the longest length transcript was used to determine the number of exons and total exon length. The intron length of a gene was calculated by subtracting the total exon length from the corresponding gene length.

Phosphorylation. Protein phosphorylation was predicted using PhosNet 3.0 [44]. The program predicts serine, threonine or tyrosine phosphorylation sites in eukaryotic proteins using neural networks. Kinase-specific predictions were used for serine and threonine, and generic predictions were used for tyrosine. In all cases, the confidence score used was > 0.7.

Signal peptides, transmembrane domains and protein localisation. Signal peptides, transmembrane domains and protein subcellular localisations were predicted using BUSCA [45]. The program runs as a web server which integrates several resources to predict sub-cellular localisation including protein feature predictors (DeepSig, TPPred3, PredGPI, BetAware and ENSEMBLE3.0) and protein localisation predictors (MemLoci, BaCelLo, SChloro). In addition, BUSCA annotates relevant protein features, such as signal/transit peptides, GPI anchors and transmembrane domains.

Gene ontology terms. GO term enrichment was obtained using the web-based application NetGO 3.0. The program predicts molecular function ontology from sequence using protein language models [46]. NetGO addresses shortcomings of GO functional annotations and is considered as one of the best methods at predicting function ontology [48,49]. Contextual visualisation of the GO terms was performed using GOATOOLS [50].

Statistical analysis. Statistical tests were carried out throughout using the statistics toolbox of MATLAB. The sequence properties did not show a normal distribution and the statistical significance of each property was determined using the two-tailed nonparametric Wilcoxon rank sum test. The test is equivalent to a Mann-Whitney U-test. Statistical significance was determined at the 0.05 level and the Bonferroni correction was applied to calculate corrected p-values. The Chi–squared (χ2) test was also carried out to check whether the frequencies of a particular feature in “stable” and “labile” genes differ from each other.

3. Results

3.1. Analysis of Genomic Features

Gene length, %CG content, transcript and exon numbers, as well as intron and exon lengths were analysed in all datasets. The results are shown in Table 4. Compared to “labile” P450s, “stable” genes were shown to be generally longer and have a greater number of exons. The opposite is true for intron lengths, where introns in the “labile” dataset are on the average considerably shorter (Table 4). The same trends are observed in the subsets of “culled” genes (Figure 1). The GC content of the “stable” CYPs is less than that of the “labile” CYPs but the difference may not be statistically significant.

3.2. Analysis of Protein Features

Protein average molecular weight, charge, isoelectric point and frequencies of different amino acid categories were analysed statistically to determine significant differences in the two groups of P450 (Table 5). “Labile” proteins were found to have greater proportions of acidic, basic, charged and polar amino acids. In the case of charged amino acids, the p-value in the culled dataset was marginally above the Bonferroni correction. “Stable” proteins had a higher median value for aliphatic amino acids and hydrophobicity. Differences between the protein features examined follow the same trend in the culled datasets (Figure 2).

Variations in the frequencies of amino acids found in “stable” and “labile” P450s are shown in Table 6. Statistically significant variations were observed in 11 amino acids (8 in the culled dataset). The relative proportion of Cys, Arg, Leu, and Trp is higher in “stable” P450s, and that of Glu, Lys, and Met higher in “labile” proteins. Variations in four of these (Ala, His) were not statistically significant in culled datasets of proteins.

3.3. Phosphorylation

Protein phosphorylation plays crucial roles in the regulation of cellular and metabolic processes such as cell differentiation and cell division. Cytochrome P450 have been shown to be subject to phosphorylation mediated by different protein kinases, cAMP-dependent protein kinase A being the most prominent one. In insects, phosphorylation is thought to play a role in the control of genotoxic metabolites [51]. P450 phosphorylation was predicted using NetPhos-3.1 [44] (SI2.xlsx). The median number of phosphorylated residues was 3 with the exception of threonine in the unculled “labile” dataset (median=4). The statistical differences between “labile” and “stable” P450s in the culled and unculled datasets were not statistically significant.

3.4. Signal Peptides, Transmembrane Domains, and Protein Location

Insect P450s are transmembrane proteins bound either to the endoplasmic reticulum (ER) or the inner mitochondrial membrane by means of the N-terminal transmembrane helix (TMH) [52]. In mitochondrial CYPs, the N-terminal anchor is missing. Instead, a topogenic sequence present at the nascent enzyme serves as a signal peptide to ensure the enzyme is transported and incorporated into the mitochondrial membrane [53]. Both transmembrane and mitochondrial P450s share the same fold and very similar tertiary structures [54]. The transmembrane domains usually adopt an α-helical structure while passing through the lipid bilayer once (single-pass proteins) or multiple times (multiple-pass proteins). In non-mitochondrial P450s, the catalytic domain of the enzyme is on the cytosolic side of ER and the N-terminus α-helix protrudes on the luminal side of ER. This posture of CYPs in relation to the membrane is stabilised by hydrophobic interactions between the residues of the N-terminal helix and the lipophilic ER environment. The N-terminal anchor has been associated with CYP trafficking into the ER or mitochondria, in interactions with different phospholipids, and as a mediator of CYP heteromer formation.

Signal peptides, subcellular localisation and number of transmembrane helices were predicted using BUSCA [45] (SI3.pdf). A summary of the results is given in Table 7. All CYPs were predicted to have at least two transmembrane domains. Approximately 80% of CYPs, both in the unculled and culled datasets, were predicted to lack the N-terminal helix. The differences between the datasets were not statistically significant.

3.5. Gene Enrichment Analysis

Gene Ontology (GO) [55]enables the classification of gene functions through the application of controlled vocabularies (ontology) to annotate the functional properties of gene and gene products across species. Each GO term is annotated with information which includes the type of gene product (e.g., protein, tRNA, etc.) and an evidence code describing the type of evidence (e.g., experimental, phylogenetic, text mining etc.). Annotations that are not curated manually are described as ‘IEA’ (inferred from electronic annotation). GO is structured as a graph comprising nodes representing each GO term. Edges between the nodes represent relationships between terms. The GO graph is hierarchical with ‘parent’ and ‘child’ terms, but unlike strict hierarchy, a given term may have more than one ‘parent’. Genes are annotated by (a) molecular function (b) cellular component and (c) biological process. One of the most common methods to perform GO term gene-set enrichment analysis is through the identification of common functions in a group of genes by means of over-representation analysis (ORA). The method compares the GO terms corresponding to a test list of genes to the functions of selected background genes (for example, a database with functionally annotated genes). If the number of functions in the test list is greater than the number of functions obtained by chance, then these functions are considered to be over-represented.

NetGO 3.0 [46] was used to predict enrichment of GO terms of the “stable” and “viable” P450s (SI4.xlsx). “Stable” and “labile” proteins were differentially enriched in the terms shown in Table 8. Although individual proteins in both datasets can show enrichment in the same GO term, the relative proportions of proteins in each set are substantially different. Comparing the numbers of “labile” and “stable” P450, considerably more “stable” P450s showed enrichment of GO terms related to biosynthetic and developmental processes. This trend is reversed in the case of deteoxification catabolic processes. GOATOOLS [50] was used to generate graphical views of enriched GO terms in biosynthetic, developmental and catabolic processes. These are shown in Figure 3 and Figure 4, respectively.

4. Discussion

The non-culled dataset was split into two subgroups which contained a total of 245 defined as phylogenetically “stable” and 241 phylogenetically “labile”. The genes were obtained from VectorBase database and cross-referred to Ensembl (Metazoa Genes 60 database). To avoid bias in statistical analyses, redundant proteins were removed from the datasets to generate non-redundant or culled datasets of “labile” and “stable” P450s. The culled and non-culled datasets were compared for a number of genomic and protein sequence properties. In general, sequence redundancy did not affect the overall picture and all sets follow similar statistical trends. Table 8 shows the gene and protein features that are strongly associated with “stable” P450 genes.

Analysis of genomic features of CYPs showed that “stable” genes are longer, have more exons, and longer intron lengths. This is consistent with recent findings showing that longer genes correlate with functions that are important in the development stages of an organism [56,57]. Additionally, “stable” CYPs the average intron length was shown to be almost twice the size of introns in “labile” P450s. This is in agreement with research showing intonic “burden” plays an important role in the evolutionary conservation of genes [58]. “Stable” CYPs were shown to have lower GC content than the “labile” ones and although the difference was not statistically significant it is noted that GC content varies inversely with intron and exon length [59].

Analysis of physicochemical and protein sequence features revealed several differences between “labile” and “stable” P450s. “Labile” proteins were found to have a greater proportion of polar and charged residues and be less hydrophobic than “stable” P450s which were relatively enriched in aliphatic amino acids. It is noted, however, that the method used to determine the relative hydrophobicity of P450s is a poor predictor as it does not take into account protein folding. It is known that the aggregation propensity of proteins is frequently linked to imperfect folding. In terms of amino acid composition, “stable” P450s were found to have Arg, Cys, Leu, Trp in greater proportion than “labile” P450s. Arg is a known chaotrope associated with the prevention of protein aggregates and increasing the solubility of hydrophobes [60]. Recent studies have shown that Trp/Tyr chains are common in enzymes utilising O₂ as a substrate and that redox active Trp/Tyr chains extend the functional lifetimes of P450s [61]. The relative enrichment of “stable” P450s in Cys residues may have a two fold importance. On one hand, Cys residues may impose structural constraints necessary for specific, rather than promiscuous, protein-ligand interactions as in the case of interactions that are developmentaly important [62]. On the other, Cys has a tendency to behave as a hydrophobic residue in folded proteins, despite possessing a polar sulfhydryl group [63].

Functional analysis of “stable” and “labile” P450s was studied by enrichment analysis of GO terms. “Stable” P450s were predicted to be enriched in GO terms relevant to the biosynthesis of lipids and hormones. These biosynthetic processes are also essential for instar larval or pupal morphogenesis and play essential roles in the growth, reproduction and defence systems of insects [64,65]. “Labile” P450s were predicted to be enriched in GO terms related to biological processes involving the detoxification of xenobiotics. As a cautionary note, it is pointed out that GO gene enrichment analysis is not without methodological challenges and pitfalls [66,67,68] and results can be unreliable due to inadequate coverage and resolution of the annotation databases [69,70].

The analysis of biological features, such as signal peptides and cellular localisation, was inconclusive and merits additional investigation. For example, proportion of “labile” CYPs predicted to be located in the mitochondrion is higher than that for “stable” CYP, albeit the fact that the absolute numbers are small to draw an definitve conclusions. This is contrary to phylogenetic studies showing that mitochondrial CYPs are relatively well conserved [22]. However, the name “mitochondrial” refering to the CYP clan should not be taken as definitive evidence of subcellular localisation and CYPs of this clan can be microsomal or have dual localisations [11]. Another important feature is the composition of tranmembrane domains. The reducing power of CYPs depends on the coupling of the endoplasmic reticulum N-terminal transmembrane helix with membrane-anchored NADPH-dependent cytochrome P450 oxidoreductase [71]. Bioinformatic analysis of of sequence conservation of transmembrane anchors and the characterisation of motif patterns of signal and transit peptides can provide useful insights on the localisation and function of CYPs [72].

5. Conclusions

“Stable” and “labile” P450s were shown to be significantly different in a number of genomic and protein features. The interdependency of these features implies that multiple aspects of biology unite to determine whether a gene is “stable” or “labile” in mosquitos. The range of features analysed and discussed in this paper can be broadened and expanded to other organisms thus gaining further insights into gene “essentiality”.Broadening the scope of the bioinformatic analysis presented in this paper may enable the functional classification of insect CYPs by training machine learning classifiers for the functional identification of insect P450s. Last but not least, the identification of “stable” genes associated with “lethality” may facilitate the identification of non-viable null mutants thus the identification of new insecticide targets.