Structural, functional and molecular dynamics analysis of cathepsin B gene SNPs associated with tropical calcific pancreatitis, a rare disease of tropics

Data curation

The single nucleotide polymorphisms (SNPs) in CTSB gene associated with TCP was extracted from an article published in 2006 (Mahurkar et al., 2006). In this article, researchers have done direct exome sequencing of CTSB gene, taking samples from 25 controls and 51 TCP patients, and further replicating the sequencing in 130 controls and 89 TCP patients from the same cohort, in order to ensure their results. In the current study, we have mapped the SNPs extracted from the literature on current human genome assembly, GRCh38.7. The mRNA accession number, NM_147782.2, and protein accession number, NP_680092.1 of gene Cathepsin B (CTSB), was used in our computational analysis. The current data about these SNVs were retrieved from human SNP database, dbSNPbuild150. The workflow for the computational analysis performed to decipher the significance of SNPs is depicted in Fig. 1.

Sequence retrieval and alignment

The sequence of Cathepsin B (CTSB) was retrieved from UniProt database: P07858 (CATB_HUMAN). The non-synonymous variants (L26V and S53G) were manually inserted in the wild-type protein sequence for further analysis.

Non-synonymous SNP analysis

The functional effect of mutations was predicted using the following algorithms: SIFT (Sorting Intolerant from Tolerant) (Ng & Henikoff, 2003), PolyPhen-2 (Polymorphism Phenotyping v2) (Adzhubei et al., 2010) and PANTHER (Mi et al., 2013). SIFT predicts whether the non-synonymous coding mutation affects protein function or not, based on sequence homology and physical properties of amino acids. SIFT calculates median conservation value for each amino acid position and thus measures the diversity of sequence. Finally, it gives the score which is the normalized probability of an amino acid change. Score of less than 0.05 are deleterious substitutions. PolyPhen-2 predicts the impact of substitution on structure and function of a protein, after annotating the substitutions and finally building conservation profiles. The prediction algorithm of PolyPhen-2 calculates Naïve Bayes posterior probability about the damaging effect of mutation and gives prediction sensitivity scores also. PolyPhen-2 annotates the substitution as “Possibly damaging”, “Probably damaging” or “Benign”, based on their scores. PANTHER predicts the functional effect of coding mutation based on “evolutionary preservation” metric of a given substitution and calculates preservation time- position-specific evolutionary preservation (PESP). Longer the PESP time, the more likely that substitution will have a deleterious effect. All these softwares: SIFT, PolyPhen-2, and PANTHER, were based on evolutionary conservation-based algorithms. Another tool, ProtParam (Bendtsen et al., 2004) was used to calculate the hydropathicity or the GRAVY (grand average of hydropathicity) score (Kyte & Doolittle, 1982) of mutated procathepsin B sequences. Hydrogen bond length and the rotational angles of main chain hydrogen bonds is a significant descriptor to study the conformation and dynamics of a protein. Therefore, to calculate the altered hydrogen bonding patterns in the mutated three-dimensional procathepsin B structures, Vadar v1.8 program (Willard et al., 2003) was used. This program calculates the H-bond distances in the main chain, side chain and the bond angle between main chain residues.

Further, the stability of the mutated models was calculated by SDM (Site Directed Mutator) server (Worth, Preissner & Blundell, 2011), I-Mutant suite (Capriotti, Fariselli & Casadio, 2005) and mCSM (Mutation Cutoff Scanning Matrix Calculation) (Pires, Ascher & Blundell, 2014) webservers. These algorithms calculate the difference in change in Gibbs free energy (ΔΔG). SDM server is used to calculate the difference in thermal stability of wild-type protein structure and mutated protein structure, using constrained environment specific substitution tables (ESSTs). I-Mutant suite is a support vector machine-based algorithm which predicts the protein stability upon mutation, by taking datasets from ProTherm (Bava et al., 2004) database. This server also calculates the change in Gibbs free energy (ΔΔG) between wild-type and mutant protein structures. mCSM server is used to predict the effect of mutations in proteins using graph-based signatures. It predicts the effect of single point mutation on protein stability by extracting data sets about various thermodynamic parameters from the ProTherm database and calculates change in Gibbs free energy (ΔΔG) between wild-type and mutant protein structures. Altogether, these three algorithms calculate ΔΔG for a protein on mutation.

The change in Gibbs free energy (ΔΔG) is as follows: (1)${\displaystyle \Delta G=\Delta H-T\Delta S}$ (2)${\displaystyle \Delta \Delta G=\Delta G_W-\Delta G_M}$ ΔG = Change in Gibbs free energy of a system from unfavourable to favourable condition ΔH = Change in enthalpy of the system ΔS = Change in entropy of the system, T = Temperature of the system ΔΔG = Value of free energy stability change of a protein upon mutation ΔG_W = Change in Gibbs free energy of the wild-type protein from unfavourable to favourable conditions ΔG_m = Change in Gibbs free energy of mutant from unfavourable to favourable conditions. ΔΔG > 0 = Increase protein stability upon mutation ΔΔG < 0 = Decrease protein stability upon mutation.

Additionally, we have also done multiple sequence alignment (MSA) of procathepsin B protein using sequences from 8 model organisms along with Homo sapiens (NP_680092.1): Mus musculus (NP_031824.1), Sus scrofa (NP_001090927.1), Macaca mulatta (NP_001181828.1), Rattus norvegicus (NP_072119.2), Ovis aries (NP_001295516.1), Danio rerio (NP_998501.1), Bos taurus (NP_776456.1). The MSA was performed using Clustal Omega program (Sievers et al., 2011).

Homology modeling

The two non-synonymous SNPs (L26V and S53G) were modelled to analyze the structural effect of variants on protein. The position specific iterated blast program (PSI-BLAST) with protein databank database (PDB) and default advanced settings, was used to find the template for homology modelling. MODELLER 9.15 (Šali & Blundell, 1993) was used to build mutated models of procathepsin B. The best predicted models according to the lowest value of DOPE score (Discrete Optimized Protein Energy), was used for further evaluation and analysis. The predicted 3D mutated models, L26V and S53G, were evaluated for their quality by using Verify-3D (Eisenberg, Lüthy & Bowie, 1997) and ProSA (Protein Structure Analysis) servers (Wiederstein & Sippl, 2007). Verify-3D examines the correctness of 3D-structures by comparing the 3D structure to the 1D structure. If 3D–1D score of each amino acid is ≥ 1, then the model is correct. ProSA analysis the correctness of theoretical models by calculating the Z-score of the input structure. ProSA considers C-alpha atoms of protein structure and calculates Z-scores based on the similarity of crystal and NMR structures of the same size. If the Z-score for a model is negative, then it is a model with minimum or no errors. After these models pass the respective thresholds, the Ramachandran plot was evaluated (Lovell et al., 2002). TM-score (Zhang & Skolnick, 2005) and root mean square deviation (RMSD) values of mutant structures (L26V and S53G), were calculated with respect to wild-type by using TM-Align web server (Berendsen, Van der Spoel & van Drunen, 1995).

Molecular dynamics simulation

The molecular dynamic simulation was performed with the Gromacs-5.0 package (Kutzner et al., 2015) on the native (PDB ID: 3PBH) and mutant structures (S53G and L26V). This computational investigation was done with a viewpoint to examine if these single nucleotide variants might lead to changes in surface properties or distort the protein orientation. The protein molecule was solvated in a dodecahedron box with SPC216 water molecules at 1.5 Åmarginal radiuses. The system was made neutral by adding 7 Na+(Sodium ions) because the initial charge of the system is −7. Subsequently, the molecular system was subjected to steepest distance energy minimization until reaching the criterion of 1,000 kJ/mol (the minimization is converged when maximum force is less than 1,000 kJ/mol) with OPLS-all atom force field (Kutzner et al., 2015). Berendson temperature coupling method (Berendsen et al., 1984) was used to regulate the temperature inside the box at 300 k. Isotropic pressure coupling was performed using Parinello-Rahman method (Martoňák, Laio & Parrinello, 2003), and the pressure of the system was maintained at 1 bar. LINCS algorithm (Hess et al., 1997) was used to treat bond lengths including H-bonds. Van der Waals and Coulomb interactions were truncated at 1nm, and Particle Mesh Ewald method (Darden, York & Pedersen, 1993) was used to compute electrostatic interaction. Finally, the simulation was performed for 35ns. The structural deviations between native and mutated structures were subjected to comparative analysis by computing RMSD (Root mean square deviation) and RMSF (Root mean square fluctuation). The trajectories were analyzed, and finally, protein compactness was studied by calculating the radius of gyration (Rg). The secondary structure analysis of wild-type and simulated mutant structures was also done using the do-dssp program of Gromacs.

Analysis of SNPs in UTR region

5′UTR and 3′UTR regions in a gene play a crucial role in regulating gene expression at the post-translational level. UTRs regulate the exit of mRNAs from the nucleus, translation efficiency, sub-cellular localization and mRNA stability (Mignone et al., 2002). The effect of SNPs in UTR regions was analyzed using databases like (1). miRdSNP (Bruno et al., 2012), (2). PolymiRTS (Bhattacharya, Ziebarth & Cui, 2014) and (3). miRNASNP (Gong et al., 2012).

Data curation

The SNPs associated with TCP were extracted from literature and is tabulated in Table 1. The SNPs were further categorized according to their type. There were total of 23 SNPs in CTSB gene found to be associated with TCP. The non-coding region included 20 SNPs (one deletion, six in 5′UTR, two in 3′UTR and 11 in introns). Coding region had 2 missense variants and 1 synonymous variant (Table 1).

Sequence retrieval and alignment

The protein sequence (NP_680092.1) was retrieved from the UniProt database, and the desired variants were manually inserted in the sequence Fig. 2. The UniProt ID for the procathepsin B sequence is P07858 (CATB_HUMAN). Mutant 1 (L26V) where residue L (Leucine) is substituted with residue V (Valine). Mutant 2 (S53G) where residue S (Serine) is substituted with residue G (Glycine).

Homology modeling

The 3D-structure of mutant (L26V and S53G) proteins was predicted after template searching by PSI-BLAST. The protein sequence (NP_680092.1) of procathepsin B was used as a query, and the resulting templates were then filtered. Finally, the X-ray crystal structure of human procathepsin B (PDB ID: 3PBH) with 2.5 Å resolution (Sequence identity: 100% and Query coverage: 93%) was used as a template for homology modeling. The DOPE scores, TM-scores, and RMSD of the predicted best models by Modeller 9.15 are shown in Table 2. The stereochemical properties of the mutated procathepsin B structures were evaluated using the Ramachandran plot from RAMPAGE. The plot defines the amino acids in favoured, allowed and outlier regions in the mutated structures as well as in the wild-type Cathepsin B structures (Table 3). Verify-3D did structure validation of predicted models, and it was observed that 99.05% amino acids had average 3D–1D protein score in a 21 residue sliding window ≤0.2 for L26V mutated model and 97.16% amino acids had average 3D–1D protein score in a 21 residue sliding window ≤0.2 in S53G mutated model. Additionally, the ProSA web server was also used to evaluate the quality of predicted 3D mutated models. The Z-score (by the ProSA webserver) of L26V model was −7.32 and of S53G model was −7.47, which were within the acceptable range of X-ray and NMR studies. The interaction energy analyzed by ProSA tool was negative for maximum residues in L26V, and S53G predicted models, in a sliding window of 10 and 40 respectively. Since the mutations were present in the propeptide region of procathepsin B, only the mutated propeptide structures are shown in Fig. 3 and Fig. 4.

Non-synonymous SNP analysis

The functional effect of the mutations predicted by using algorithms described in the “Methods” section are tabulated in Table 4. Analysis by SIFT and PANTHER suggests that S53G and L26V mutations can have damaging effects. The GRAVY scores of wild-type WT was −0.470, and the mutants (L26V and S53G) was −0.469 for both. Thus, it could be concluded that there is no significant effect of mutations on hydropathicity of the protein. The comparative analysis of hydrogen bond lengths and the rotational angles between WT and mutants were calculated using Vadar v1.5 server, at 10 different regions of protein, playing imperative role in the functioning of procathepsin B. Remarkable differences in the H-bond lengths and Bond angles between WT and mutants (S53G and L26V) were observed as represented in Fig. 5. Hence, it could be interpreted that these mutations alter the binding between residues in mutated structures.,

The change in Gibbs free energy (Δ ΔG) of the mutated structures calculated by SDM server, I-Mutant 2.0, and mCSM webservers indicatedestabilization of mutated proteins (Table 5). Additionally, the results of MSA revealed that leucine at the 26th position and Serine at the 53rd position are conserved residues. Extrapolating these results indicate that any alteration in conserved residues will affect the structure and function of the protein (Fig. 6). Altogether, L26V and S53G mutations were predicted to have a deleterious effect on the structure and function of the protein, through non-synonymous SNP analysis algorithms.

Analysis of SNPs in UTR region

UTRs play an essential role in mRNA processing during post-translational mechanism. Hence, the SNPs in the UTR region can significantly affect the functionality of UTRs, provided they affect the miRNA binding sites. The 3′UTR region is essential for microRNA (miRNA) binding which can lead to degradation or transcriptional suppression of mRNA and thus, can further affect the downstream processing. The databases miRdbSNP, PolymiRTS, miRNASNP, were used to predict the significance of SNPs in the 3′UTR region (Table 6). The two SNPs present in 3′UTR region, rs709821 and rs8898, were predicted to be present in miRNA binding sites and therefore are significant. The SNP, rs8898 was predicted to create a new miRNA site, while rs709821 disrupt a non-conserved miRNA site as predicted by PolymiRTS database. The two miRNAs targeting CTSB gene having both non-coding variants, hsa-miR-96 and hsa-miR-1271, are pancreas specific miRNAs, deciphered from the miRNet database (Fan et al., 2016).

Molecular dynamics simulations

The comparative analysis of trajectories by calculating RMSD, RMSF and Radius of gyration after MD simulation of 35 ns, for both native and mutants, was performed. Interestingly, it was observed from RMSD of backbone residues, that both mutated structures (S53G and L26V) were conspicuously deviated from the native structure (PDB ID: 3PBH). To infer effect of mutations on the dynamic behaviour of each residue, RMSF of C α-atoms was calculated. There was fluctuation observed in mutated protein structures as compared to the native structure. The protein compactness was determined by the radius of gyration (Rg). It was observed that Rg of mutated structures were distinctly fluctuated as compared to native structure, throughout the simulation (Fig. 7). The secondary structure analysis of wild-type and simulated mutant structures, by do-dssp program, implicated that mutations had caused deviation to the protein (Fig. 8).