2.1 ζ- crystallin sequences

The sequence of C. dromedarius ζ-crystallin fetched from European Nucleotide Archive (EMBL accession no: GADU01072870.1). The accession numbers of Homo sapiens, Cavia porcellus, and Hyla japonica ζ-crystallin were BAD96870.1, NP_001166407.1, and BAB41213.1, respectively.

2.2 Multiple sequence alignment and phylogenetic analysis relationships

The amino acid sequences were compared with those obtained from the protein database using PSI-BLAST algorithm. Next, homologous sequences of ζ- crystallin were aligned using the MAFFT sequence alignment technique (Katoh and Standley 2013). Next, the alignment output was color-coded using the Jalview workbench (Waterhouse et al., 2009). Finally, evolutionary divergence among the homologous sequences and the phylogenetic tree were constructed using the BLOcks SUbstitution Matrix (Henikoff and Henikoff 1992).

2.3 Sequence analysis

The physicochemical properties of C. dromedarius ζ-crystallin and its sequence homologs such as molecular weight, aliphatic index, theoretical pI, grand average hydropathy (GRAVY), and instability index (II) were calculated using the ProtParam tool (Gill and von Hippel 1989).

2.4 Homologous modeling

The three-dimensional structures of ζ-crystallin belonging to C. dromedarius, C. porcellus, and H. japonica were predicted using the Swiss-model server (Biasini et al., 2014). The crystal structure of H. sapiens ζ-crystallin (PDB id: 1YB5) at 1.85 Å resolution was used as a template. Homology modeling was carried out based on the following criteria: (1) Sequence identity > 60%; (2) 100% coverage; and (3) Template structure with a bound ligand. All the structures were analyzed and processed using PyMOL (PyMOL, 2006).

The three-dimensional structure models of C. dromedarius, C. porcellus, and H. japonica ζ-crystallin were validated using PROCHECK (Morris et al., 1992), ProSA-Web (Sippl 1995, Wiederstein and Sippl 2007) and QMEAN (Benkert et al., 2009) servers.

2.5 Molecular docking

The modeled 3D structures of ζ-crystallin belonging to all species in the trial were used for docking studies. NADPH (ZINC85551904) ligand was taken from the Zinc15 database and converted to a PDB format. AutoDock Vina (Trott and Olson 2010) was used to determine the binding site between ligands and proteins. Initially, potential binding sites were detected based on blind docking where the boxes were sufficiently large enough to cover the entire protein. PyMOL and Discovery Studio Visualizer (Biovia 2015) were used to visualize the structures.

2.6 Molecular dynamic simulations

The GROMACS 5.1.2 version was used to carry out the MD simulation for 50 ns at 300 K under a pressure of 1 bar using a GROMOS96 43a1 force-field. Computational tools such as PyMOL, Visual Molecular Dynamics (VMD) (Humphrey et al., 1996), and QtGrace were used to visualize, evaluate, and analyze MD trajectories. The structural coordinates of the ζ-crystallin of all species and NADPH were pre-processed in SPDBV. The topology and force-field parameters for NADPH were generated from the PRODRG server and merged into the parent.gro file of species’ ζ-crystallin to make the complex. Energy minimization was carried out using the steepest descent to remove bad connections to prevent structural hindrances in the solvated systems. The final MD ran at 50 ns for both systems, and the resulting trajectories were saved for further analysis using inbuilt utilities of GROMACS such as gmx rms, gmx rmsf, gmx gyrate, and gmx sasa (Naqvi et al., 2020).

3.1 Multiple sequence alignment of ζ-crystallin

The amino acid sequence of ζ-crystallin of C. dromedarius was 100% identical to the two-hump domesticated Bactrian camel (C. ferus) (accession no: XP_006178374.1). Sequence alignment of C. dromedarius ζ-crystallin with H. sapiens, C. porcellus, and H. japonica showed 87, 83, and 64% identity, respectively (Table 1). The sequence homology was 79–97% between C. dromedaries ζ-crystallin and other mammalian sources. Average sequence homology of the thirteen mammalian ζ-crystallins was 87%. Whereas the sequence homology of amphibian (H. japonica) ζ-crystallin was 64%. The multiple sequence alignment of ζ-crystallin from various sources is presented in Fig. 1. The output showed a high degree of conservancy.

Close

Fig. 1

Amino acid sequence alignment of 13 mammalian and one amphibian ζ- crystallins. The multiple sequence alignment was generated with the MAFFT Multiple Sequence Alignment on a Jalview workbench. According to conservancy, residues are color-coded. Ligand (NADP) binding residues are labeled by (*).

Export to PPT

3.2 Phylogenetic analysis of ζ-crystallin

The phylogenetic relationships among ζ-crystallins from various sources and evolutionary divergence are shown in Fig. 2. The C. dromedarius ζ-crystallin grouped with that of Lama guanicoe and Sus scrofa. H. sapiens and the other mammals showed evolutionary divergence. Due to low similarity, H. japonica (ζ-crystallin) was separated from other mammalian species.

Close

Fig. 2

The phylogenetic tree of camel ζ-crystallin and homologous genes. The phylogenetic tree was generated with the BLOSUM62 from MAFFT Multiple Sequence Alignment.

Export to PPT

3.3 Physicochemical properties of ζ-crystallin

The physicochemical properties of C. dromedarius ζ-crystallin and its sequence homologs are shown in Table 1. The calculated pI of C. dromedarius, C. porcellus, and H. japonica ζ- crystallin was close to neutral pH (7.1, 7.69, and 6.76, respectively). The ζ-crystallin was expressed at very high levels in the lens of C. dromedarius, C. porcellus, and H. japonica constituting approximately 10% of the total lens soluble proteins in these animals (Garland et al., 1991). Except for Bos taurus (bovine) ζ-crystallin, the pI of all other mammalian ζ- crystallin is in the alkaline region (pI > 8.0), suggesting that the solubility of the proteins is minimal at this pI. Interestingly, ζ-crystallin made up most of the soluble lens proteins in C. dromedarius, C. porcellus, and H. japonica, and their isoelectric points were close to the physiological pH. It is noteworthy that maintaining such a high degree of solubility is required throughout their lifespan. It may be possible that α-crystallin acts as a molecular chaperone to assist ζ-crystallin in maintaining its solubility. The pH inside the lenses of these animals is not known. It is also possible that the pH inside the lenses is different in these animals or some other factors may enhance their solubility.

The Grand Average of Hydropathy (GRAVY) value of a protein indicates its solubility (Khan et al., 2018). A positive value is indicative of a hydrophobic nature, while a negative value points to the hydrophilic nature of the protein (Kyte and Doolittle 1982). The ζ-crystallin of C. dromedaries, L. guanicoe, M. musculus, and H. japonica was relatively more hydrophilic than the other species (Table 1). C. porcellus ζ-crystallin was more hydrophobic than that of C. dromedarius, M. musculus, and H. japonica. In earlier studies, it was observed that when C. porcellus ζ-crystallin was expressed in E.coli, it was largely recovered as insoluble inclusion bodies. But, when C. porcellus ζ-crystallin was co-expressed with the GroEL/ES chaperone system, the solubility increased markedly (Goenka and Rao 2001). In contrast, M. musculus and H. japonica ζ-crystallin were expressed in E.coli without any assistance from a chaperone system (Fujii et al., 2001, Simpanya et al., 2010). Similarly, C. dromedarius recombinant ζ-crystallin was produced in E.coli without a chaperone system in our lab (Malik et al., 2016, Malik et al., 2017). Moreover, the yield of pure recombinant M. musculus ζ-crystallin was five-fold higher than recombinant C. porcellus ζ-crystallin (Simpanya et al., 2010), which could be due to the hydrophilic nature of the M. musculus ζ-crystallin.

The instability index provides an estimate of the in vitro stability of the test protein. An instability index value < 40 indicates that the protein is stable (Guruprasad et al., 1990). The instability index of all three ζ-crystallins examined were < 40. The instability index of the M. musculus ζ-crystallin (31.89) was the highest, while H. japonica ζ-crystallin was the lowest (22.28). The aliphatic index (AI) is used to indicate the thermostability of globular proteins (Ikai 1980). All (ζ-crystallins had a high aliphatic index, which varied between 88.05 and 104.33.

3.4 Structural studies of ζ-crystallin

The physiological role of ζ-crystallin in lenses is not well understood. Truncated mutants of ζ-crystallin cause cataracts in C. Porcellus. The structure of ζ-crystallin from C. dromedarius, C. porcellus, and H. japonica ζ-crystallin is not known, but a high-resolution structure of homologous H. sapiens ζ-crystallin (PDB: 1YB5) is available. A Swiss model server generated a 3D model structure of the (ζ-crystallin of the three species under study. Ligand (NADPH) bound human ζ-crystallin was superimposed on the ζ-crystallin of C. dromedarius, C. porcellus, and H. japonica (supplementary Fig. S1). The overall topology of ζ-crystallin from these animals was similar to human ζ-crystallin. The sequence identity of H. sapiens ζ-crystallin with C. dromedarius, C. porcellus, and H. japonica was 87, 83, and 64%, respectively.

Different parameters of the modeled structures were evaluated using PROCHECK, ProSA-Web, and QMEAN servers. The PROCHECK analyses evaluate the stereo-chemical quality of the protein structure. The PROCHECK analysis was generated using a Ramachandran plot. As shown in supplementary Fig. S2, this plot ensured the quality of the model structure of C. dromedarius, C. porcellus, and H. japonica ζ-crystallins with 92.5, 91.8 and 91.6% of the residues in the ‘most favored region’, respectively. According to the Ramachandran plot, 7.1, 7.8, and 8% of the residues were in the ‘additional allowed region’. Furthermore, only 0.4% of the residues were in the disallowed regions in all model structures.

The QMEAN Z-score evaluates the degree of nativeness in the model structure by comparing it to the template structure solved by X-ray crystallography or NMR spectroscopy. The QMEAN Z-scores range between 0 and 1. The QMEAN analysis of all three models is shown in Supplementary Fig S3. The three plots showed that the Z-score QMEAN of C. dromedarius, C. porcellus, and H. japonica were 0.48, 0.76, and 0.4, respectively. The red marker in the QMEAN graph of C. dromedarius and H. japonica ζ-crystallin structure is located within the dark zone while in the graph of C. porcellus it is located on the boundary of the dark zone, indicating an overall good quality of all modeled structures.

The ProSA analysis was used to analyze the overall quality of the three structure models of C. dromedarius, C. porcellus, and H. japonica ζ-crystallin by calculating Z-score based on the plot of experimentally determined protein structure (Supplementary Fig. S4). The Z-score for C. dromedarius, C. porcellus, and H. japonica ζ-crystallin structures were −8.71, −8.44, and −9.09, respectively. The Z-score in all three modeled structures displayed in the plot was in the range of all experimentally determined protein chains in the current PDB. If the ProSA Z-score is located outside a characteristic range of native conformation of proteins, then the model or experimental structures are erroneous (Sippl 1995).

3.5 Molecular docking

The docking results showed that NADPH has good docking scores with ζ-crystallin of the three species. The ligand, NADPH, showed strong binding affinities of −9.3, −10.8, and −9.8 kcal/mol towards C. dromedarius, C. porcellus, and H. japonica ζ-crystallin, respectively (Table 2). As shown in Fig. 3A, the C. dromedarius ζ-crystallin-NADPH complex interacted with the aid of seven hydrogen bonds within a distance range of 1.98–2.98 Å and six hydrophobic interactions within a distance range of 3.04–5.10 Å. The interaction between C. porcellus ζ-crystallin and NADPH) was relatively stronger than between the NADPH complexes with C. dromedarius and H. japonica ζ-crystallin. It occurred due to NADPH forming nine hydrogen bonds with distances between 1.93 and 3.39 Å and thirteen hydrophobic bonds at distances between 2.99 and 5.37 Å with the side chains of C. porcellus ζ-crystallin(Fig. 3B). The interaction between NADPH and H. japonica ζ-crystallin was relatively stronger than C. dromedarius ζ-crystallin as it formed eleven hydrogen bonds (with distances between 2.69 and 3.38 Å) and ten hydrophobic interactions (with distances between 3.33 and 5.34 Å) (Fig. 3C).

Close

Fig. 3

Interaction between ζ-crystallin and NADPH through molecular docking. (A) animated structure of C. dromedarius ζ-crystallin binding with NADPH in the central calyx as the best position of binding, (B) animated structure of C. porcellus ζ-crystallin binding with NADPH in the central calyx, (C) animated structure of H. japonica ζ-crystallin binding with NADPH in the central calyx. The H-bonding is shown in a dashed green line and hydrophobic interactions in a dashed pink line.

Export to PPT

3.6 Molecular dynamic study of ζ-crystallin and ligand

MD simulations have been used extensively to understand the conformational changes, stability, and dynamics of the binding mechanisms of ligands to a protein under explicit solvent conditions (Al-Shabib et al., 2019a,b, Khan et al., 2019). Here, the stability of the complexes (C. dromedarius, C. porcellus, and H. japonica ζ-crystallin-NADPH docked complexes) were evaluated by performing molecular dynamic simulations for 50 ns at 300 K. Binding of a small molecule in the binding pocket of a protein can lead to large conformational changes (Ali et al., 2014, Al-Shabib et al., 2020). Fig. 4 shows the root mean square deviations (RMSDs) of C. dromedarius, C. porcellus, and H. japonica ζ-crystallin with and without the ligand. The RMSD analysis provides valuable information on the conformational changes in structure during the simulation. We have determined the RMSD across Cα-atoms of all three ζ-crystallins in the absence and presence of NADPH. As shown in Fig. 4A, it is apparent that the RMSD of C. dromedarius ζ-crystallin without the ligand increased initially and stabilized within 2 ns of simulation. Conversely, the RMSD of the C. dromedarius ζ-crystallin-NADPH complex stabilized within 15 ns and remained constant during the 15–50 ns of simulation time, indicating a stable complex formation. The average RMSD values of C. dromedarius ζ-crystallin in the absence and presence of NADPH were estimated to be 0.34 and 0.49 nm, respectively, within the acceptable limit of 0.2 nm. The average RMSD values for C. porcellus ζ-crystallin in the absence and presence of NADPH were found 0.52 nm and 0.39 nm, respectively (Fig. 4B). Similarly, the average RMSD values for H. japonica ζ-crystallin alone or complexed with NADPH were found to be 0.44 nm and 0.42 nm, respectively (Fig. 4C). Here, in all three simulations, the RMSD values fluctuated by<0.2 nm, confirming the stability of the ζ-crystallin-NADPH complex during the 50 ns of simulation. After adding NADPH, the RMSD values of C. dromedarius, C. porcellus, and H. japonica ζ-crystallin complexes were stabilized in approximately 15, 2, and 5 ns, respectively. Stabilization of the ζ-crystallin-NADPH complexes depended on the complexes' binding energies, and results indicated that C. dromedarius ζ-crystallin-NADPH > C. porcellus ζ-crystallin-NADPH > H. japonica ζ-crystallin -NADPH (Table 2). Overall, the RMSD plots of C. dromedarius, C. porcellus, and H. japonica ζ-crystallins showed equilibration throughout the 50 ns MD simulation, suggesting stable NADPH complexes' formation (Fig. 4).

Close

Fig. 4

Characterization of ζ-crystallin-NADPH interaction parameters as a function of MD simulation time. (A) RMSD of C. dromedarius ζ-crystallin alone and complexed with NADPH, (B) RMSD of C. porcellus ζ-crystallin alone and complexed with NADPH, (C) RMSD of H. japonica ζ-crystallin alone and complexed with NADPH.

Export to PPT

Root mean square fluctuation (RMSF) measures the flexibility in the protein 3D structure due to side-chain movements. This study estimated the flexibility and rigidity of the ζ-crystalline side chains in the absence and presence of NADPH by measuring the RMSF during a 50 ns simulation. In Fig. S5A, it is observed that the average RMSF values of C. dromedarius ζ-crystallins alone and complexed with NADPH were 0.16 and 0.19 nm, respectively. The average RMSF of C. porcellus ζ-crystallins in the absence or presence of NADPH was 0.18 and 0.19 nm, respectively (Fig. S5B). Similarly, the average RMSF of H. japonica ζ-crystallins with and without NADPH was also 0.18 and 0.19 nm, respectively (Fig. S5C). The fluctuating regions on the ζ-crystallins of C. dromedarius, C. porcellus, and H. japonica spanned across the protein molecule regardless of the presence of NADPH. These residual fluctuations were found to be minimized upon binding of the NADPH throughout the simulation. However, several random residual fluctuations were also observed when ζ-crystallin bound to NADPH (Fig. S5).

The radius of gyration (Rg) measures the overall compactness of a protein–ligand complex during MD simulation. Moreover, Rg analyses also provide information on the protein folding pattern, stability, and conformational changes in the absence and presence of ligands. Here, we have assessed the Rg of C. dromedarius, C. porcellus, and H. japonica ζ-crystallins alone and as complexes with NADPH as a function of simulation time (Fig. S6). The average Rg values of C. dromedarius ζ-crystallins alone or complexed with NADPH were determined to be 3.01 and 2.78 nm, respectively. A 0.23 nm decrease in C. dromedarius ζ-crystallin-NADPH complex indicates that the binding of NADPH increased the compactness of C. dromedarius ζ-crystallins (supplementary Fig. S6A). NADPH binding with C. porcellus and H. japonica ζ-crystallins lead to a slight increase in the Rg of both ζ-crystallins indicating that the binding of NADPH resulted in a partial loss of compactness in both ζ-crystallins (Table 3 and Fig. S6B & C).

Measurement of the solvent-accessible surface area (SASA) of a protein shows the changes in solvent exposure of amino acid residues after ligand binding (Khan et al., 2019, Al-Shabib et al., 2020). This study assessed how the SASA was affected by the ligand binding of ζ-crystallins obtained from the three species under study at different simulation times (Fig. S7). We observed that the mean SASA of C. dromedarius ζ-crystallins alone and in its complex form was 281.98 and 259.14 nm², respectively (Fig. S7A). A decrease in the SASA of the C. dromedarius ζ-crystallin-NADPH complex indicated that the solvent surface exposure area decreased with NADPH binding. The mean SASA values for C. porcellus alone and when complexed with NADPH were found to be 250.73 and 266.15 nm², respectively (Table 3). An increment in C. porcellus (-crystallins) SASA was observed with NADPH binding (Fig. S7B). This increment in the SASA can be presumed to result from the exposure of some of the internal residues in C. porcellus ζ-crystallins to solvents due to conformational changes arising from complexing with NADPH. The solvent surface exposure area of H. japonica ζ-crystallins increased when bound to NADPH, indicating that minor conformational changes occurred after NADPH binding (Table 3 and Fig. S7C). After the simulation, the volume and density of the (-crystallins alone and when complexed with NADPH were also calculated (Table 3).