Screening of high-risk deleterious missense variations in the CYP1B1 gene implicated in the pathogenesis of primary congenital glaucoma: A comprehensive in silico approach

Data retrieval

The list of missense variations in the CYP1B1 gene of Homo sapiens was retrieved from the dbSNP database of the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/projects/SNP/). The encoded protein sequence with reference ID: NP_000095.2 was also downloaded from the same database.

Functional annotation of variants

The experiment is carried out computationally and is called the “in silico method.” These methods use databases, network analysis, and machine learning approaches that are faster, more critical, and more cost-effective than in vitro and in vivo studies (Azfaralariff et al., 2022). Therefore, in the present study, to identify the high-risk deleterious variations, the functional annotation of missense variations in the CYP1B1 gene was carried out by employing multiple bioinformatics tools. Seven different tools, which were SIFT, PolyPhen-2, PROVEAN, SNAP2, PANTHER, PhD-SNP, and Predict-SNP, were used in this screening assay. The first screening tool employed in this assay to investigate the deleterious variants from the data retrieved from the dbSNP NCBI database was called SIFT (sorting intolerant from tolerant) (https://sift.bii.a-star.edu.sg/). It predicts the functional importance of the altered amino acid substitution using evolutionary conservation via sequence homology alignment and the physiochemical characteristics of a protein. The SIFT algorithm predicts substitutions with a probability score ≥ 0.05 as tolerated substitutions, while those with a probability score <0.05 are considered deleterious (Ng & Henikoff, 2003; Shahid et al., 2020). Polyphen 2.0 (polymorphism phenotyping v2) (http://genetics.bwh.harvard.edu/pph2/) was utilized, acknowledging its structural and functional impact prediction of a single amino acid alteration in a protein. This tool evaluates an amino acid substitution using comparative homology protein models and generates a position specific independent count (PSIC) score. It annotates a substitution as probably damaging, possibly damaging, or benign if the probability score is >0.85, >0.15, and <0.15, respectively (Adzhubei et al., 2010). Further, the PROVEAN (Protein Variation Effect Analyzer) (http://provean.jcvi.org/index.php) server was accessed in this functional annotation assay. This server predicts the functional impact of an amino acid alteration in a protein and categorizes them as neutral or deleterious, if the generated score is greater than −2.5 or less than −2.5, respectively (Choi et al., 2012). The next screening was performed by the PANTHER (Protein analysis through evolutionary relationship-position-specific evolutionary preservation) (http://www.pantherdb.org/) server. This tool calculates the score of an amino acid variant based on evolutionary preservation (time in millions of years) using the HMM algorithm. It annotates an amino acid alteration as benign: time <200 my, possibly damaging: 450 my to 200 my, or probably damaging: 450 my, respectively (Tang & Thomas, 2016). SNAP2 (Screening of non-acceptable Polymorphism) (https://rostlab.org/services/snap2web/) tool was employed to determine the pathogenicity of an amino acid alteration in a protein. It is an ensemble of neural network models that predicts both the clinical and molecular effects of a missense variant in the form of strongly effective (score: +100) or strongly neutral (score: −100) (Bromberg & Rost, 2007). Moreover, the PhD-SNP (Predictor of human deleterious single nucleotide polymorphism) tool, accessed from https://snps.biofold.org/phd-snp/phd-snp.html, was also utilized in this functional assay. It is a Support Vector Machine (SVM) tool that predicts the pathogenicity of an amino acid substitution in a protein based on evolutionary information and a hybrid predictive model. This tool classifies the missense variants as neutral or pathogenic and provides a reliability index (RI) ranging from 0-10, with 10 being the highest reliability (Capriotti, Calabrese & Casadio, 2006). Additionally, the Predict-SNP (https://loschmidt.chemi.muni.cz/predictsnp/) tool was also approached to predict the disease related variations of the CYP1B1 gene. It uses consensus-based methods to generate its output (Bendl et al., 2014). The allelic frequency of the variants was determined using the gnomAD (Genome Aggregation Database) version 2.1.1 (https://gnomad.broadinstitute.org/). The gnomAD is the most widely used database that harbors the largest collection of sequencing data for population variations (Gudmundsson et al., 2022).

Effect of variations on CYP1B1 protein stability

To evaluate the stability of the protein upon any amino acid alteration, two machine learning tools, namely, I-Mutant 2.0 and MUpro, were utilized. I-Mutant 2.0 (https://folding.biofold.org/i-mutant/i-mutant2.0.html) predicts protein stability using an SVM-based algorithm and calculates the amino acid change in the form of Gibbs free energy (DDG) at a given pH and temperature. It categorizes the protein stability as increased or decreased and provides a reliability index (RI) ranging from 0 to 10, with 10 being the highest reliability (Capriotti, Fariselli & Casadio, 2005). The input conditions were set to pH 7 and temperature 37 °C respectively. Additionally, the MUpro server (http://mupro.proteomics.ics.uci.edu/) was approached to validate the I-Mutant 2.0 outcomes. This server ensembles both SVM and neural network algorithms to predict the effect of amino acid alteration in a protein and classifies the protein stability as increased or decreased if the confidence scores are <0 or >0, respectively (Cheng et al., 2006).

Conservancy analysis of high-risk variations and secondary structure of CY1B1

The evolutionary conservation of the high risk deleterious amino acid variations was evaluated by utilizing the ConSurf server (https://consurf.tau.ac.il/). This server works on the empirical Bayesian method and multiple sequence alignment approach to estimate the phylogenetic conservation score, which ranges from 1–9. Residues of conservation scores of 7–9 were considered highly conserved, 4–6 as intermediate, and 1–3 were represented as variable and less conserved. In addition, this tool also predicted the characteristics of residues based on their positions in the protein’s structure: either buried (b), exposed (e), highly conserved and exposed (functional, f), or highly conserved and buried (s), which further embellished the structural and functional significance of amino acid residues in a protein (Ashkenazy et al., 2010; Shahid et al., 2022). For the secondary structure prediction of CYP1B1 protein, the SOPMA (Self-optimized Prediction Method with Alignment) bioinformatic tool was utilized by using its default query parameters (window width, 17; similarity threshold, 8; and number of states, 4). This server ensembles the self-optimized prediction method and a neural network method (PHD) to generate the secondary structure of a protein (Geourjon & Deleage, 1995).

Structural examination of native and mutant types of CYP1B1

The three-dimensional (3D) structure of CYP1B1 protein was downloaded from the Protein Data Bank (PDB) (https://www.rcsb.org/). To generate the complete structure of the CYP1B1 protein, the I-TASSER (Iterative Threading Assembly Refinement) (https://zhanggroup.org/I-TASSER/) most advanced biological server was approached (Zhang, 2008). This is a tool that constructs the 3D structures of proteins through hierarchical approaches by first searching the templates from the PDB database and then constructing the models through iterative template fragment assembly simulation. Moreover, 19 mutant structures of CYP1B1 protein were also generated through the I-TASSER tool by incorporating the amino acid alteration of each high-risk deleterious variation into the wild-type FASTA protein sequence. The quality (degree of appropriateness) of all the native and mutant structures was verified through the ERRAT and VERIFY3D, which are in-built tools of the SAVES meta-server version 6.0 (https://saves.mbi.ucla.edu/). Subsequently, the resultant structures were further subjected to the TM-Align server (https://zhanglab.dcmb.med.umich.edu/TM-align/) which calculates the RMSD (Root Mean Square Deviation) to analyze the deviation of mutant structures from the wild-type (Zhang & Skolnick, 2005). Graphical inspection and superimposition of wild-type and mutant structures were performed by UCSF Chimera software version 1.14 (Pettersen et al., 2004).

Molecular dynamics simulation (MDS) analysis

The molecular dynamic simulation (MDS) for the native and mutant CYP1B1 protein structures was carried out using GROMACS simulation software accessed from https://simlab.uams.edu (Bekker et al., 1993). For the MDS run, the wild-type and mutant protein structures were prepared using the GROMOS96 43a1 force field solvated by the Single Point Charge (SPC) water model and fixed in a periodic cubic solvated box. Other requisite parameters in molecular dynamic simulation, such as neutralization by adding the salt (NaCl) concentration (0.15 M), energy minimization for 5,000 steps of the prepared system using the steepest descent method, equilibration type consisting of NVT/NPT at 300 K and 1 bar of pressure using the leap-frog method. The final MD simulation step was processed for 50 ns to analyze the deviation of mutant structures from the wild-type protein structure.

Investigation of CYP1B1-interacting candidates

The co-expression analysis of the CYP1B1 gene was evaluated through the GeneMANIA server (https://genemania.org/). It is a publically available bioinformatics tool that harbors a large database for gene functional analysis, including co-expression, pathway interaction, physical interaction, co-localization, genetic interaction, and shared protein domains (Shahid et al., 2021; Warde-Farley et al., 2010). The submitted query was run by setting Homo sapiens as the species organism and the maximum number of resultant genes and attributes as 10. On top of that, the protein-protein interaction (PPI) of CYP1B1 was also determined by utilizing the STRING protein database (https://string-db.org/). The STRING database query terms were cut-off at 0.7 and the maximum additional interaction was 10 (Azfaralariff et al., 2022; Szklarczyk et al., 2021). Cytoscape software version 3.8.2 (https://cytoscape.org/) was employed to construct the highly co-expressed genes and PPI networks in the present study (Shannon et al., 2003).

Gene ontology and pathway enrichment analysis of CYP1B1

The functional annotation of three gene ontology (GO) terms, which were molecular functions (MF), cellular components (CC), and biological processes (BP), and the KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway of CYP1B1 were analysed with Enrichr database (https://maayanlab.cloud/Enrichr/), a publically available online tool for gene expression profiling, transcriptional factor analysis, PPI networks, pathway analysis, toxicogenomics and pharmacogenomics studies (Chen et al., 2013). The P-value <0.05 was set as a cut-off criterion for the submitted query. Cytoscape software was utilized to construct all the networks.

Data retrieval and selection

The dbSNP is a unique database that holds all the necessary information, including the variation data of many of the human genes. This database with CYP1B1 reported a total of 4788 SNPs, of which 610 were missense, 276 were synonymous, 1359 were in the intronic region, and the rest belonged to other categories (Accessed: December 21, 2021). In this study, the missense variations were selected in order to investigate their damaging effects and how they influence the phenotypic characteristics of the respective encoded protein. These missense variations are responsible for various complex diseases by affecting the protein’s functional activities and structural stability or folding (Khalid & Almaghrabi, 2020; Shahid et al., 2022). A list of all 610 missense variants is presented in Table S1.

Functional annotation of missense variants

The missense variations of a gene can produce amino acid allelic variants that may influence the functions of that particular gene product by disrupting the structural conformation. Hence, missense variants of the CYP1B1 gene were chosen for further investigation. Variations of 610 CYP1B1 missense retrieved from the dbSNP NCBI database were first subjected to the SIFT server, which predicted 92 variants to be either deleterious or tolerated (Table 1; Table S1 ). These 92 variants were further evaluated by six other computational tools, namely PolyPhen-2, PhD-SNP, SNAP2, PANTHER, PROVEAN, and Predict-SNP, which used different algorithms to scrutinize a variant. The predictions and categories of these tools are provided in Table 1. From conducted analysis, a total of 19 missense variants (P52L, G61E, G90R, P118L, E173K, D291G, Y349D, G365W, G365R, R368H, R368C, D374N, N423Y, D430E, P442A, R444Q, F445L, R469W, and C470Y) were found to be common in all the functional analysis tools by manual concordance, and all were deleterious, probable damaging, and effect variants. Therefore, these 19 variants were considered the high-risk deleterious variants of the CYP1B1 gene and are promising candidates for further research. The detailed results of 19 high-risk variants are presented in Table 2 and their allele frequencies information are in Table 3.

Effect of variations on CYP1B1 protein stability

Two biological tools, namely I-Mutant 2.0 and MUpro, were independently used to investigate the effects of all the 19 missense variations on the CYP1B1 protein stability. The I-Mutant 2.0 server predicted the effect of amino acid change in the form of protein stability Gibb’s free energy change (DDG) and calculated the RI of the variations, while MUpro predicted the effect of protein stability for single site mutation in the form of Gibb’s free energy change (DDG). Of the 19 high-risk missense variations identified, I-Mutant 2.0 revealed that all amino acid variations profoundly affected the CYP1B1 protein by decreasing its stability, except for four variants: G61E, P118L, Y349D, and D430E, which did not influence the protein stability. Apart from I-Mutant 2.0, MUpro identified only P52L had no effect on protein stability and the remaining exhibited a reduction in protein stability. The details of I-Mutant 2.0 and MUpro analysis are provided in Table 3 and Table S2.

Conservancy analysis of high-risk variations and secondary structure of CY1B1

The evolutionary conservation characteristic of all the common 19 missense variants was accomplished by the Consurf server. It is essential to determine the magnitude of an amino acid substitution to alter the protein’s structure and function. Results revealed that all amino acid variations occurred in the highly conserved region of the CYP1B1 protein except for one substitution (Y349D) that was observed in the variable region (Fig. 1). Additionally, the solvent accessibility feature predicted that 13 variants (P52L, G90R, P118L, E173K, D291G, G365W, G365R, R368H, R368C, D374N, D430E, P442A, and R444Q) were found to be highly conserved and exposed residues, four variants (G61E, N423Y, F445L, and C470Y) were observed to be highly conserved and buried residues and two variants (Y349D and R469W) were identified as buried and exposed residues. The evolutionary conservation and solvent accessibility feature of all the 19 substitutions are demonstrated in (Fig. 1). Generally, it has been considered that exposed residues have an impact on protein function and, thereby, highly conserved and exposed residues are recognized to have functional importance. On the other hand, buried residues are incorporated into the protein structure maintenance, and thus highly conserved and buried residues reveal their structural importance (Zhang et al., 2017). Based on this analysis, it can be concluded that changes in 19 regions of the CYP1B1 protein, which play an important role in the structural and functional behavior of protein integrals, confirm their deleterious nature.

Furthermore, of the total 543 amino acids, 255 amino acids (46.96%) were found in the alpha helix, 207 amino acids (38.12%) in the random coils region, 61 (11.23%) in the extended strand, and 20 (3.68%) in the beta turns of the CYP1B1 protein as predicted by SOPMA (Fig. S1).

Structural examination of native and mutants of CYP1B1

To predict the structural impact of these 19 high-risk variations, first the 3D structure of the CYP1B1 protein was obtained from the PDB database (PDB ID: 3PM0). The structure had some missing residues; therefore, to generate the complete structure of the CYP1B1 protein, the I-TASSER server was utilized. From the generated models of I-TASSER, the best model with a C-score of 0.16, an estimated TM-score of 0.73 ± 0.11, and an estimated RMSD of 7.0 ± 4.1 A was selected and subjected to further evaluation. The quality of the best generated model structure was validated through the ERRAT and VERIFY3D tools of the SAVES server. The ERRAT tool provided the overall quality score, which was 92.57. The VERIFY3D tool estimates the compatibility of a protein model based on local statistical potentials and predicts that the CYP1B1 protein model passed the criteria with an average score of 89.94% and is worth further investigation. Similarly, the 3D structures of all 19 mutants of CYP1B1 were generated, and their graphical illustration and superimposition with the wild-type structure were performed by UCSF Chimera software version 1.14. The 3D wild-type and some mutant structures are depicted in Fig. 2A. In addition, the RMSD of all 19 mutant structures was calculated using the TM-Align server, and their results are shown in Table 3. The results of the TM-Align server elucidated that all 19 mutant structures have profound impacts on the CYP1B1 protein. This suggests that when these variations occur in the CYP1B1 protein, they may interfere with its functional behavior.

Molecular dynamics simulation (MDS) analysis

The molecular dynamic simulation (MDS) was executed to evaluate the protein stability following a missense variation in the CYP1B1 protein. The MDS was run up to 50 ns and the RMSD values were noted for the wild-type and for all the 19 high-risk missense variations. Results demonstrated that the CYP1B1 wild-type structure showed an RMDS in the range of 0.36 nm and structural stability throughout the simulation period (Fig. 3). These results indicate that the native structure attained a stable conformation in the MDS period. Overall, the RMSD value was observed to be higher in all the 19 mutants. Six variants (P52L, D291G, R368H, R368C, D374N, and C470Y) demonstrated an RMSD ∼0.38 nm. The highest value of RMSD was observed in three variants (Y349D, D430E, and F445L) that had ∼0.52, 0.46, and 0.44 nm, respectively, compared to the wild-type protein structure. The remaining variants’ RMSD was observed to be greater than 0.40 nm. Figure 3 depicts the RMSD results of all 19 variants after the MDS run. Generally, the higher the RMSD value indicates, the higher the deviation of the mutant with respect to the native protein structure (Agrahari et al., 2019). Therefore, the results of MDS demonstrated that these variants were predicted to destabilize the native protein structure, which could interfere with the normal functioning of the CYP1B1 protein.

Investigation of CYP1B1-interacting candidates

Co-expression data for the CYP1B1 gene was obtained from the GeneMANIA database and then imported into the Cytoscape software to build a network based on the co-expression score (Fig. 4A). Results elucidate that CYP1B1 was highly co-expressed with CYP1A1, followed by PDCD6-AHRR and VCAN genes, respectively. The least expression of CYP1B1 was observed with the LUM, PLA2G7, and F13A1 genes. Furthermore, the protein-protein interaction (PPI) network was also generated via cytoscape software (Fig. 4B). Then, the cytoHubba plugin was employed to identify the highly interacting proteins with CYP1B1. Through cytoHubba, it was discovered that five candidates, which were: CYP1A1, HSD17B1, HSD17B2, AKR1C3, and HSD17B7, demonstrated close interaction with CYP1B1 protein (Fig. 4C).

Gene ontology (GO) and pathway enrichment analysis of CYP1B1

In this study, the Enricher server was utilized for Gene Ontology (GO) and pathway analyses. In total, 59 biological processes (BPs) were discovered for CYP1B1 as documented in Table S3. The p-value for all 59 BPs was less than 0.05, indicating that CYP1B1 plays an important biological role in these processes. Fig. 5A illustrates the top ten BPs. Subsequently, only four molecular functions were discovered, which were estrogen 16-alpha-hydroxylase activity, oxidoreductase activity, steroid hydroxylase activity, and heme binding, respectively (Fig. 5B). Likewise, the cellular component (CC) was widely observed in the endoplasmic reticulum membrane and intracellular membrane-bounded organelles. The detailed results of GO are tabulated in Table S3. Moreover, six KEGG pathways were discovered for CYP1B1, which were: Tryptophan metabolism, ovarian steroidogenesis, steroid hormone biosynthesis, metabolism of xenobiotics by cytochrome P450, chemical carcinogenesis, and microRNAs in cancer Fig. 5C Table S4).