Discovery of Klotho peptide antagonists against Wnt3 and Wnt3a target proteins using combination of protein engineering, protein–protein docking, peptide docking and molecular dynamics simulations

Abstract

The Klotho is known as lifespan enhancing protein involved in antagonizing the effect of Wnt proteins. Wnt proteins are stem cell regulators, and uninterrupted exposure of Wnt proteins to the cell can cause stem and progenitor cell senescence, which may lead to aging. Keeping in mind the importance of Klotho in Wnt signaling, in silico approaches have been applied to study the important interactions between Klotho and Wnt3 and Wnt3a (wingless-type mouse mammary tumor virus (MMTV) integration site family members 3 and 3a). The main aim of the study is to identify important residues of the Klotho that help in designing peptides which can act as Wnt antagonists. For this aim, a protein engineering study is performed for Klotho, Wnt3 and Wnt3a. During the theoretical analysis of homology models, unexpected role of number of disulfide bonds and secondary structure elements has been witnessed in case of Wnt3 and Wnt3a proteins. Different in silico experiments were carried out to observe the effect of correct number of disulfide bonds on 3D protein models. For this aim, total of 10 molecular dynamics (MD) simulations were carried out for each system. Based on the protein–protein docking simulations of selected protein models of Klotho with Wnt3 and Wnt3a, different peptides derived from Klotho have been designed. Wnt3 and Wnt3a proteins have three important domains: Index finger, N-terminal domain and a patch of ∼10 residues on the solvent exposed surface of palm domain. Protein–peptide docking of designed peptides of Klotho against three important domains of palmitoylated Wnt3 and Wnt3a yields encouraging results and leads better understanding of the Wnt protein inhibition by proposed Klotho peptides. Further in vitro studies can be carried out to verify effects of novel designed peptides as Wnt antagonists.

Keywords:

• Homology modeling
• Klotho
• molecular dynamics (MD) simulations
• peptide design
• protein engineering
• protein–protein docking
• Wnt3
• Wnt3a

Introduction

The aging comes with various health complications caused by degenerative processes such as osteoporosis, decreased fertility, arteriosclerosis and skin atrophy1^,2. It is often conjectured that the stem and progenitor cell function abrogation and cell senescence contribute to aging. Absence of Klotho is found to be related to stem and progenitor cell function abrogation and ultimately to cell death3. Klotho, which includes 1012 amino acid residues, is a trans-membrane (TM) protein with only 10 residues belonging to intracellular (IC) and 925 residues located to extracellular (EC) domains. Interestingly, along with being membrane bound, the Klotho also acts as secreted protein in circulation4. EC part of the Klotho consists of two domains (KL1 and KL2)5. These domains share similarity to Family I glycosidases. The secreted Klotho is believed to play an important role in longevity of the protein6.

The stem cell regulators (i.e. Wnt proteins) are postulated to have functional interaction with the Klotho7. The Wnt3 and Wnt3a (Wingless-type Mouse Mammary Tumor Virus (MMTV) integration sites 3 and 3a) proteins are 40 kDa in size Wnt family proteins. Wnt signaling has been recognized for its importance in embryonic development8. The embryonic processes orchestrated by Wnt signaling include cell proliferation, cell migration, cell fate specification and tissue regeneration in adult bone marrow9^,10. The Wnt signaling is often involved in stem cell control as a self-renewal signal11. Aberrant Wnt genes or Wnt pathway components lead to various developmental defects. Increased uninterrupted Wnt signaling is being evident in Klotho deficient animals which causes stem cell and progenitor cell senescence12. The aberrant expression of Wnt proteins is found to be antagonized by Klotho. The single EC domain KL1 (amino acids 1–285) interacts with Wnt3 and Wnt3a, and mediates the response that antagonizes the aberrant expression of these proteins. It is believed that Klotho is the secreted antagonist of Wnt3 and Wnt3a7.

Due to the highly hydrophobic nature of Wnt proteins only two proteins, so far only xWnt8 and WntD, have been crystalized. The structural studies of secreted lipid modified Wnt proteins revealed some highly unusual structures in C-terminal domain (CTD) and N-terminal domain (NTD), forming folds which mimic the index finger and thumb as well as palm like features, respectively. Another structural and functional distinctions in Wnt protein are the Cys–Cys disulfide bonds. The members of Wnt family proteins have 22–25 Cys amino acid residues distributed throughout the structure, mostly involved in disulfide linkage13^,14. These Cys–Cys disulfide bonds function to maintain the natural folds of the protein. The Cysteine residues distribution in Wnt3 protein is shown in Figure 1. The protein sequence alignment of important domains of xWnt8 with Wnt3 and Wnt3a is given in Table 1.

Figure 1. Distribution of 24 Cysteine residues over different domains of Wnt3 protein structure. The finger like domain, including Cysteine knot and thumb domain are Cysteine hotspots. The position of Cysteine residues in the protein structure gives the clue about disulfide bonds importance in the stability of critical Wnt protein domains (NTD: N-terminal domain; CTD: C-terminal domain).

Table 1. Frizzled receptor/co-receptors interaction sites. xWnt8 annotation is based on its experimental values. Wnt3 and Wnt3a annotation is based on sequence similarity studies considering xWnt8 as template.

Important interaction sites	xWnt8	Wnt3	Wnt3a
Thumb domain residues	181–188	206–213	203–210
Finger domain residues	261–338	332–342	329–339
A patch of ∼10 residues on palm domain	216–219	241–244	238–241
	249–252	273–276	270–273
	256–259	280–283	277–280
Palmitoylation site	Ser187	Ser212	Ser209

Keeping in mind the importance of Klotho in Wnt signaling, in this study combined in silico approaches have been applied to better understand the important interactions between Klotho and Wnt3/Wnt3a. The main aim of the work is to identify important residues of the Klotho (in form of peptide) that could potentially act as Wnt antagonist using peptide docking techniques with three important domains of Wnt3 and Wnt3a proteins (i.e. index finger domain, thumb domain and a patch of ∼10 residues at solvent exposed palm domain)15–17.

Methods

Homology model building

Homology models of Wnt3, Wnt3a and Klotho proteins were prepared using SWISS-MODEL, I-TASSER and Prime module of Schrodinger’s Maestro molecular modeling Suit18–20. xWnt8 (PDB ID: 4F0A) was used as template protein for the 3D models of Wnt3 and Wnt3a target structures using 34% sequence identity. Model evaluation was done by Ramachandran plot determined by Maestro21. Quantitative comparison between the structure of studied proteins and the native state of its template was carried out by PROCHECK and ERRAT 2.022^,23. Along with the model evaluating software and tools, major part was played by secondary structure analysis of the proteins in selecting potential model and number of formed disulfide bonds within the protein structure.

Klotho protein modeling

The Wnt binding domain is located within the amino-terminal portion of Klotho’s KL1 domain (amino acids 1–285). Therefore, the KL1 domain of Klotho protein was considered for generating model and further experimentation. The homology model was generated using a template Klotho-related protein (Cytosolic neutral β-glycosylceramidase (PDB ID: 2E9M)). Sequence identity between template and target sequences was found as 46%. The generated model was further directed to energy minimization and molecular dynamics (MD) simulations.

Molecular dynamics (MD) simulations

Molecular dynamics simulations played an important role in the evaluation of potential model of Wnt and Klotho proteins by studying the dynamical behavior of proteins. MD simulations were performed by Gromacs 5.1 package24^,25. To completely immerse the protein systems in water, initial structures of the proposed systems were solvated with simple point charge (SPC) model water molecules in a cubic box of 10 Å following periodic boundary conditions (PBC). GROMOS96 43a1 force field26 was used for all MD simulations. The constructed system was then neutralized by “genion tool” plugin of Gromacs package. The energy minimization for the whole systems were carried out by steepest descent (SD) algorithm for 50 000 iterations. Particle Mesh Ewald (PME) method was used to calculate the long-range electrostatic interactions. A cutoff radius of 14.0 Å was applied for van der Waals and Coulomb interactions. The studied systems were simulated under isothermal-isobaric ensemble (NPT) with temperature of 300 K and pressure of 1 atm using velocity rescaling (modified Berendsen) temperature coupling and Parrinello–Rahman pressure coupling methods27^,28. The equilibrations (NVT and NPT) were performed for 50000 iterations. After implementation of minimization and equilibration strategies, total 10 systems of Wnt3, Wnt3a and Klotho were subjected to 20 ns each (totaling 200 ns) production run of MD simulations using 2 fs as time step. Coordinates of all atoms were recorded in the form of .trr Gromacs trajectory file for every 10 ps (i.e. 2000 frames each simulation). Post-processing MD simulations analyses including root mean square deviation (RMSD), root-mean square fluctuation (RMSF) were done by visual molecular dynamics (VMD) program29.

Protein–protein docking

In order to understand the important residues more viable of making interaction with Wnt proteins, the docking of Klotho KL1 domain with the Wnt3 and Wnt3a was carried out. The ClusPro and HADDOCK docking programs were used for protein–protein docking simulations with default parameters30^,31. The resulting docked complexes from both docking programs were analyzed to understand interactions of Klotho amino acid residues with Wnt proteins. Furthermore, three new peptides in Klotho’s KL1 domain were predicted by Data Mining for Enzymes Search Utility (DME) which were expected having glycosidase activity32.

Peptide preparation

The marked strings of (adjacent) amino acid sequences making good interactions with Wnt3 and Wnt3a proteins were manually cleaved from the generated model PDB file of KL1 domain and saved into separate coordinate pdb file as peptide. These peptides were prepared for peptide docking using Maestro molecular modeling package33. Appropriate charged ends were added to each peptide followed by peptide docking.

Protein preparation

All the proteins and peptides were prepared for the in silico experimentation using protein preparation module of Schrodinger’s Maestro Molecular modeling Suit. Hydrogen atoms were added followed by energy minimization and optimization using OPLS2005 force field. Protonation states were located at physiological pH 7.4 using PROPKA34^,35.

Peptide docking

Peptide docking was carried out by Glide SP-Peptide docking and Induced Fit Docking (IFD) modules of Maestro molecular modeling package36–38. The grid files for the receptors Wnt3 and Wnt3a were generated by selecting appropriate parameters for different set of fields for better representation of shape and properties of the receptor in Glide SP-Peptide docking39. The grid covers residues of interest in Wnt3 and Wnt3a proteins, palmitoylated Ser212 and Ser209, respectively, for thumb domain. Likewise, the grid files were generated for important residues at the active sites of finger and palm domains (Table 1). The Glide protocol is further used for SP-Peptide docking into the defined binding site of protein given in grid file. IFD provide the selected approach using Glide (XP) and Prime tools giving extra conformational sampling for peptide and protein binding site. The principle approach used by IFD involves the docking of peptides into target binding sites using Glide SP. The target binding site and peptide van der Waals scaling is set to 0.50 Å by default. The resulting peptide–protein complexes with high docking scores are then selected. The residues residing within the 5 Å of docked peptides in the target site are further refined by Prime module of Maestro and side chains were optimized. Eventually, the peptides are re-docked into the refined binding site of protein using Glide (XP) docking protocol.

Protein

Palmitoylation

Since the Wnt3 and Wnt3a proteins are active, and secreted after palmitoylation, further experimentations on the lipid modified protein models were also performed. Palmitoleic acid is the chemical ingredient which makes the Wnt3 and Wnt3a mature when attached at conserved Serine residues (Ser212 and Ser209, respectively)40^,41. For palmitoylation, the Palmitoleic acid was docked to the specific residues position of the Wnt protein. Glide/SP, Glide/XP and Glide/IFD modules of Maestro were used for this aim. After getting the reasonable coordinates near the mentioned residues of the Wnt proteins, covalent bond was created manually between the palmitoleic acid and conserved Serine residue. The covalent bond formation was further proceeded by energy minimization using Macromodel module of Maestro. The proteins mimicking its mature form were prepared for docking studies of Klotho derived peptides.

Results

Model evaluation

Various important points were considered for the best protein model evaluation including secondary structure elements, and appropriate number of disulfide bonds in the protein structure.

Importance of secondary structure elements

A 20 ns-long trajectories was derived on pre-selected Wnt3 (Wnt3-M0, Wnt3-Model-0) and Wnt3a (Wnt3a-M0, Wnt3a-Model-0) protein models to determine the importance of secondary structure in the protein model selection. The selected homology model that is generated by Prime module of Schrodinger’s Maestro molecular modeling Suit, comprises of finger domain (β-strand hairpin loop) with few patches of β-sheets, as well as an α-helix at its base and six disulfide bonds in each models. The MD simulations analysis shows 180° movement of the finger domain in both proteins as well as abrupt movement of thumb domain (see Figure S1, Supplementary Materials, as well as Animation-1). The 20 ns MD simulations was also performed for the Wnt template protein xWnt8 (PDB entry 4F0A) to compare the results (control). The core of template protein xWnt8 remained very stable during the MD simulations span. The comparison and in-depth analysis of the protein structure such as RMSD, RMSF and visual intimation (Supplementary Material, Animation-1) indicate the importance of secondary structure elements (β-sheets, α helices) in the Wnt protein family CTD and NTD42. Keeping in mind the results we gain in form of deformed Wnt protein models (Wnt3-M0 and Wnt3a-M0) at the final step of MD simulations, we are unable to use these models for further experiments. The new homology models (Wnt3-M1 and Wnt3a-M1) have been created from the scratch. This time, models have been selected by comparing the sites where secondary structure elements are important (as the structural flaws showed by bent Wnt models), with xWnt8 protein. These models have the substantial number of β-sheets and α-helices in overall protein structure, specially, β-strand hairpin loop structure of finger domain and Cysteine knot. These two models (Wnt3-M1 and Wnt3a-M1) were subjected to 20 ns MD simulations to observe their structural and dynamical behaviors. Unlike previous models (Wnt3-M0 and Wnt3a-M0), the new models this time did not deform. This is the first step toward selecting appropriate protein model for Wnt proteins.

Importance of disulfide bonds in protein

The disulfide bonds play important role in the protein structure stability. The previously described in silico experiment with Wnt3-M0 and Wnt3a-M0, resulting in bent posture of β-strand hairpin loop structure of finger domain indicates apparently the importance of α-helix and β-sheets in the protein structure. The next step was to understand the importance of number of disulfide bonds in the Wnt protein family. The substantial number of Cysteine residues in the Wnt protein family enhance the significance of disulfide bonds in Wnt3 and Wnt3a protein structures. The new favorable models selected, keeping in mind the indications availed from models regarding secondary structures, were used for further experimentation to introduce disulfide bonds. The homology models generated by SWISS-MODEL, I-TASSER and Prime were taken into account and assessed one by one for protein models with higher number of secondary structure elements (β-sheets, α helices) in the CTD and NTD of Wnt3 and Wnt3a. The model selection procedure was carried out with the comparison of xWnt8 (4F0A) as well. Its observed that the selected Wnt3 homology model (Wnt3-M1) included two disulfide bonds which are generated automatically by the SWISS-MODEL, I-TASSER and Prime. Likewise, selected homology model for Wnt3a (Wnt3a-M1) involved five disulfide bonds. Indeed, there is uncertainty about the number of disulfide bonds present in the structure of the Wnt proteins. The template protein xWnt8 has 11 disulfide bonds in pdb structure, however, it is reported to have 12 disulfide bonds14. The one missing disulfide bond is due to the missing residues at the terminal point of the protein, which is difficult to crystallize. In some studies, Wnt3a is reported to have 12 disulfide bonds while UniProt shows only nine disulfide bonds (UniProtKB Wnt3 ID: P56703, Wnt3a ID: P56704)43^,44. There is no exact information available about the number of disulfide bonds within Wnt3 protein. Therefore, we conducted three different experiments each on Wnt3 and Wnt3a proteins to understand the importance of disulfide bonds. In the first experiment, we used the selected homology models including default disulfide bonds formed in the processes of homology modeling (two disulfide bonds in model Wnt3: Wnt3-M1 (Wnt3-Model-1) and five disulfide bonds in model Wnt3a: Wnt3a-M1 (Wnt3a-Model-1)). For the second experiment, UniProt data was used as reference, and nine disulfide bonds between the Cysteine residues were used as mentioned for both proteins resulting in models Wnt3-M2 (Wnt3-Model-2) and Wnt3a-M2 (Wnt3a-Model-2). For the third experiment, we have manually asserted 12 disulfide bonds in Wnt3 and Wnt3a inferred from their sequence similarity with xWnt8 and resulting models named as Wnt3-M3 (Wnt3-Model-3) and Wnt3a-M3 (Wnt3a-Model-3). It must be noted that a total of 12 disulfide bonds has been experimentally reported for the proper functioning of the Wnt3a45. The Wnt3 and Wnt3a proteins share 84.2% sequence identity. The reason behind devising experiments on Wnt3a model with Wnt3 for number of disulfide bonds is to understand structural and dynamical behaviors of highly identical proteins (Wnt3 and Wnt3a) under same in silico experimental conditions. This behavioral analysis with Wnt3a anticipated to help to understand the appropriate number of disulfide bonds in Wnt3. All these six models (i.e. three models for Wnt3 and three models for Wnt3a) were then directed toward 20 ns MD simulations using Gromacs.

RMSD and RMSF analysis of xWnt8

The RMSD and RMSF graphs clearly indicate that protein is stabilized during simulations. The RMSF shows most fluctuating residues within CTD i.e. residues from 300 to 338 and residues from 150 to 200, and 100 to 150 in NTD (Figure 2). xWnt8 is the template protein which is used for the 3D models of Wnt3 and Wnt3a target structures. The comparative analysis of homology models of Wnt protein was carried out considering MD simulations analysis of xWnt8.

Figure 2. RMSD and RMSF analysis of xWnt8 protein MD simulations. Average RMSD for xWnt8 is 4.2 Å.

RMSD and RMSF analysis of Wnt3 models

Model Wnt3-M1: The two disulfide bonds in Wnt3-M1 give the protein final structure after MD simulations with finger and thumb domain shifts away from the palm domain as compared to its initial stage structure before MD simulations, as visualized by VMD, Maestro and Chimera46 (Figure 3). The average RMSD remains 4.7 Å for 20 ns MD simulations for this model (Figure 4). Since the CTD has a β-strand hairpin comprising of 40 amino acids, it needs stability which is provided by extensive network of disulfide bonds. During initial 4 ns of MD simulations, the sudden adjustment and fluctuation of Cysteine knot have been observed. That could be the one reason behind massive fluctuations in finger domain residues as seen in RMSF graph. According to the RMSF graph, the four major domains show high fluctuations, and these domains comprise of amino acid residues ranging from 28 to 65 with two Cysteine residues, 140 to 165 with two Cysteine residues, 205 to 220 with four and 306 to 354 with 10 Cysteine residues (Figure 4).

Figure 3. Wnt3 comparison of change in domains in all three experiments from initial structure. Alignment of these structures have been done by Needleman–Wunsch Alignment algorithm and BLOSUM-62 matrix in Chimera molecular visualization tool by keeping initial stage structure as reference. Right side 90° angle view of finger domain has been presented. Cysteine knot shows huge difference among all three models like other domains. The thumb domain appears comparatively stable. (Yellow and red colored cartoons show before (yellow) and after (red) MD simulations of Wnt3 with 12 disulfide bonded model (Wnt3-M3). Blue and green-colored cartoons show after MD simulations of model with nine disulfide bonds (Wnt3-M2) and model with three disulfide bonds (Wnt3-M1), respectively).

Figure 4. RMSD and RMSF analysis of three Wnt3 protein homology models (Wnt3-M1: Wnt3 homology model with two disulfide bonds; Wnt3-M2: Wnt3 homology model with nine disulfide bonds; Wnt3-M3: Wnt3 homology model with 12 disulfide bonds). Average RMSDs for Wnt3-M1, Wnt3-M2 and Wnt3-M3 are calculated as 4.7 Å, 7.6 Å and 4.3 Å, respectively.

Model Wnt3-M2: Increasing the number of disulfide bonds to nine in the Wnt3-M2 model gives interesting movements in the thumb and prominently finger domain. The thumb domain moves abruptly throughout the simulation as seen in previous Wnt3-M1 experiment with two disulfide bonded model, as well. The finger domain bends ∼90° sidewise from the initial position for Wnt3-M2. The RMSD graph shows high level of fluctuation as average of 7.6 Å throughout the simulation span (Figure 4). The RMSF graph shows that the highly fluctuating residues in this experiment are the same as in the previous experiment with two disulfide bonded model (Wnt3-M1). However, the intensity of fluctuation is relatively high, where domain spanning 306–354 residues have fluctuations of 12.6 Å (Figure 4).

Model Wnt3-M3: Interestingly, the structures with 12 disulfide bonds show stable system compared to other two, as shown in RMSD and RMSF graphs (Figure 4). However, unexpectedly, the finger domain, irrespective to the number of disulfide bonds in the structure, bends inward to the palm domain in the middle of MD simulations span, with frequent movements in whole domain (Figure 3; red cartoon). Similarly, thumb domain has shown some arbitrary movements during the MD simulations. The finger domain tried to get stable position while bending towards the palm domain which is different configuration as compared to the template (xWnt8) throughout the simulations (Supplementary Material, Animation-2). The RMSD analysis gives slightly lower values as compared to the two disulfide bonded model Wnt3-M1 with average RMSD of 4.3 Å. Furthermore, the RMSF analysis shows the most fluctuating amino acids are from 300 to 354 with RMSF of 10 Å which is lower than the experiment with nine disulfide bonded model Wnt3-M2 and slightly higher than the Wnt3 experiment with two disulfide bonded model Wnt3-M1 (Figure 4).

RMSD and RMSF analysis of Wnt3a models

Model Wnt3a-M1: Structural instability is evident during the MD simulations span when the Wnt3a protein with five disulfide bonds is experimented (Wnt3a-M1). An unexpected fold has been observed during MD simulations. Both the finger and thumb domains bend inward making a hole as seen in protein structure (Figure 5; Green cartoon). The RMSD value reaches to 8.7 Å during 17–19 ns of MD simulations. The average calculated RMSD was 6.8 Å for this model (Figure 6). The most fluctuating domains can be explained by RMSF values of individual residues, as four major domains exhibit maximum fluctuations. The domains comprise of residues 20–60, 128–160, 180–230 and 306–351 with number of Cysteine residues in each domain, 2, 4, 4 and 10, respectively.

Figure 5. Wnt3a comparison of change in domains in all three experiments from initial structure. Alignment of these structures have been done by Needleman–Wunsch Alignment algorithm and BLOSUM-62 matrix in Chimera molecular visualization tool by keeping initial stage structure as reference. Right side 90° angle view of finger domain has been presented. Cysteine knot shows huge difference among all three models like other domains. The thumb domain appears comparatively stable. (Yellow and red colored cartoons show before (yellow) and after (red) MD simulations of Wnt3a with 12 disulfide bonded model (Wnt3a-M3). Blue and green-colored cartoons show after MD simulations of model with nine disulfide bonds (Wnt3a-M2) and model with three disulfide bonds (Wnt3a-M1), respectively).

Figure 6. RMSD and RMSF analysis of three Wnt3a protein homology models (Wnt3a-M1: Wnt3a homology model with five disulfide bonds; Wnt3a-M2: Wnt3a homology model with nine disulfide bonds; Wnt3a-M3: Wnt3a homology model with 12 disulfide bonds). Average RMSDs for Wnt3a-M1, Wnt3a-M2 and Wnt3a-M3 are 6.8, 5.4 and 5.5 Å, respectively.

Model Wnt3a-M2: The in silico experiment with nine disulfide bonds shows different behavior of Wnt3a protein as compared to that of Wnt3. Unlike Wnt3-M2, the finger domain in Wnt3a-M2 does not bend ∼90° sidewise (Figure 5; blue cartoon). Both finger and thumb domains keep on fluctuating throughout the MD simulations span as seen from the RMSD and RMSF analyses (Figure 6). The MD trajectory analysis shows the average RMSD of 5.4 Å throughout the simulation span (Figure 6). The RMSF graph indicates same four distinct domains as of Wnt3a experiment with five disulfide bonds (Wnt3a-M1) showing fluctuations maximum at CTD with 10 Å RMSF. The comparison of RMSD graphs of Wnt3a models with five and nine disulfide bonds show that the average RMSD value is slightly lower for the nine-disulfide bonded model (5.4 Å).

Model Wnt3a-M3: The model with 12 disulfide bonds also showed fluctuations during the MD simulations (Supplementary Materials, Animation-3). The RMSD and RMSF values show the fluctuation and deviations in the protein structure during MD simulations. The comparative analysis among all three models indicates that the RMSD values for Wnt3a with 12 disulfide bonds (Wnt3a-M3) remain lower at the second half (11 ns to 20 ns) of the MD simulations span, compared to the other two models. The average RMSD for this experiment remains 5.5 Å throughout the MD simulations which is greater than the average RMSD 4.2 Å of xWnt8 (Figure 6). The RMSF analysis of 12 disulfide bonded model shows the similar fluctuations in same four domains as seen in previous experiments of Wnt3a with and nine disulfide bonded model (Wnt3a-M2). However, the residues starting from 30 to 50 amino acids show slightly abrupt pattern of fluctuation compared to other models (Figure 6).

The analysis of MD trajectories of xWnt8, Wnt3 and Wnt3a clearly demonstrates that the CTD is always found to be less stable domain compared to the other parts of the protein in all models irrespective of number of disulfide bonds. The MD simulations of other family members of Wnt proteins (Wnt4) have also shown the most fluctuating residues residing in CTD15.

Addition of palmitoleic acid

As per our theoretical analysis of the structural aspects of the Wnt3 and Wnt3a proteins, the model with 12 disulfide bonds were selected for lipid modification (Wnt3-M3 and Wnt3a-M3). The target protein files of selected models for docking are representative structures of trajectory frames which has the lowest RMSD to the average structure within 18–20 ns of MD simulations span. The average structure is calculated by RMSD trajectory tool of VMD. Palmitoleic acid is covalently attached to the conserved Serine residues, i.e. Ser212 and Ser209 in Wnt3 and Wnt3a, respectively.

Protein–protein docking simulations

The resulting docked complexes from both docking programs, Cluspro and HADDOCK were analyzed to understand interactions of KL1 domain of Klotho residues with Wnt proteins. The strings of adjacent residues of Klotho making interactions with the Wnt proteins were considered for catalytic domain of Klotho and these domains were used in peptide preparation (Figure S2, Supplementary Materials). Accumulatively, eight strings of residues (Peptide 4–11) of KL1 domain of Klotho were marked as making good interactions with Wnt3 and Wnt3a. Furthermore, three new peptides (Peptide 1–3) were also predicted having glycosidase activity using DME32. Thus, in total 11 peptides were taken for docking simulations with important domains of Wnt3 and Wnt3a (Table 2).

Table 2. Klotho derived peptides docking (IFD) against three important binding sites of Wnt3 and Wnt3a. Docking scores are in kcal/mol.

Peptide ID	Peptide residue sequence in Klotho’s KL1 domain	Peptidesequence	Docking analysis of Wnt3			Docking analysis of Wnt3a
			Thumb domain/PAM	Finger domain	Co-receptorbinding site	Thumb domain/PAM	Fingerdomain	Co-receptorbinding site
			1	147–153	YRFSISW	−1.9	−10.7	−11.4	−7.5	−4.1	−9.1
2	204–211	YGGWANRA	−5.1	−12.9	−10.4	−6.6	−6.2	−9.9
3	231–237	VKYWITI	−6.4	−11.2	−8.3	−5.2	−4.3	−8.9
4	140–147	RELGVTHY	−6.3	−5.3	−7.2	−6.1	−5.5	−9.8
5	70–75	GSAAYQ	−7.5	−9.1	−11.5	−7.0	−8.6	−10.2
6	60–66	TFPDGFL	−8.9	−8.1	−7.4	−6.3	−7.1	−9.0
7	233–240	YWITIDNP	−6.3	−6.9	−9.8	−6.6	−7.8	−8.8
8	233–239	YWITIDN	−4.9	−8.0	−11.9	−8.0	−7.0	−9.7
9	283–287	NTSFR	−8.3	−9.0	−8.3	−5.6	−7.9	−9.3
10	265–269	YLVAH	−5.4	−7.1	−10.1	−5.3	−7.0	−9.1
11	265–272	YLVAHNLL	−5.4	−7.7	−8.9	−5.4	−9.1	−8.2

Peptide docking

The docking scores of docked complexes of Klotho derived peptides with Wnt3 and Wnt3a are represented in Table 2. All peptides show good docking scores (IFD scores) and interactions with all three important domains including palmitoylated Serine and surrounding residues in thumb domain, important residues in finger domains and co-receptor binding sites present in palm domain. For instance, the peptide-1 gives docking score of −10.7 kcal/mol and −11.4 kcal/mol against finger and palm domains of Wnt3, respectively. Similarly, the peptide-5 shows docking sores of −7 kcal/mol, −8.6 kcal/mol and −10.2 kcal/mol against the thumb, finger and palm domains of Wnt3a, respectively (Table 2). Interacting residues of all three domains of Wnt3 and Wnt3a with peptides are given in Table 3. The docking studies of all 11 Klotho-derived peptides come up with very good docking scores and interactions with the important residues at least one of the three binding sites of Wnt protein (Figure 7). The docking poses of peptides-4 and 11 have been shown for Wnt3 and Wnt3a, respectively (Figure 7). These results highlight the importance of potential peptides which can be used instead of full length Klotho protein as Wnt protein antagonist to cease its aberrant expression.

Figure 7. Interaction of peptide with all three considered domains of Wnt3 and Wnt3a using induced fit docking (IFD) docking approach. (a) Docked complex of Wnt3 (green) with peptide 4 (magenta) against three binding pockets (cyan). (b) Docked complex of Wnt3a (cyan) with peptide 11 (magenta) against three binding pockets (yellow). Blue-colored lines show palmitoleic acid attached to Wnt proteins.

Table 3. Interacting residues of Wnt3 and Wnt3a with peptides derived by IFD docking approach.

Peptide ID	Peptide sequence	Docking analysis Wnt3			Docking analysis Wnt3a
		Thumb domain/PAM	Finger domain	Co-receptor binding site	Thumb domain/PAM	Finger domain	Co-receptorbinding site
		1	YRFSISW	Glu216 His209	Arg129 Glu160 Glu133 Lys205 Trp221 His335	Arg88 Arg110 Ser240 Asp273 Leu265 Glu243 Tyr276 Asn279 Ser280	Leu208 Glu213 Ser209 Lys204	Arg328, His332, Lys326	Arg247, Glu240, Met274, Ala238, Tyr262
2	YGGWANRA	Ser212, Glu216, Leu211, Cys208	Ser28, Tyr23, Cys338, Gly32, Glu71, Glu344, Gln33, Val340, Lys329	Arg88, Arg110, Ala241, Glu243, Tyr276, Ser242, Asp273	Lys204, Cys205, Gly210, Cys212, Glu213, Ser211	Arg344, Lys326, Glu325, Gln340, Ser338, Cys339, Glu341	Lys234, Asp236, Asp233, Asp270, Ser277, Tyr273, Glu240, Ser237
3	VKYWITI	Lys154, His209, Leu211, Gly152, Glu216	Tyr23, Ser28, Glu71, Arg85, His335, Glu133	Arg88, Arg110, Glu243, Pro268, Asp273, Asp239, Tyr276, Asn282, Asn279, Phe283,	Ser209, Glu213	His332, Phe331, Trp333, Cys329, Lys326	Arg84, Arg85, Glu275, Glu240, Tyr273, Ser277
4	RELGVTHY	Lys207, Gly210, Glu216, Ser214, Ser212	Cys332, Phe334, Trp336, Gln343, Ser341	Arg88, Arg110, Glu111, Asn279, Ser280, Glu243, Ser242, Pro268	Cys212, Leu208, Glu213, Ser209	Arg328, Arg344, Gln340, Cys339, Lys326, Glu341	Arg84, Arg85, Lys232, Asp236, Ser239, Met241, Glu240, Tyr273, Arg107, Tyr235
5	GSAAYQ	Lys154, His209, Leu211, Glu216	Glu71, Trp26, Trp27, Leu29, His335tyr339	Arg88, Arg110, Ser280, Asn279, Glu243, Ala241, Ser240, Asp273	Cys212, Ser211, Glu213, Leu208	Val330, Tyr336, Phe331, Glu325, Arg344, Arg328, Cys329	Glu240, Arg247, Glu275, Glu108, Arg107, Ala238, Tyr273
6	TFPDGFL	Lys154, His209, Leu211, Glu216, Lys207, Ser212, Gly213	Glu71, Arg129, Glu133, Gln224, Lys205	Arg110, Glu243, Asn279, Asn282	Ser209, Glu213, Cys212, Leu208	Tyr336, Arg328, Phe331	Glu240, Ala238, Tyr273, Arg107, Tyr274, Glu275
7	YWITIDNP	Lys154, Lys207, Cys208	Gln33, Phe334, His335, Trp336, Gln343, Ser341	Arg88, Arg110, Glu243, Tyr275, Asp273, Ser242	Lys204, Gly210, Cys212, Glu213	Phe331, Cys329, Arg328, Cys339, Tyr336, Lys326, Val330	Glu240, Ser239, Tyr273, Ala238, Arg85, Ser277, Arg107, Arg84, Asp236
8	YWITIDN	Glu216, Cys215	Gly32, Gln33, Phe334, Trp336, Gln343, Lys329, Glu344, His331	Arg8, Arg110, Tyr238, Ser280, Phe283, Glu243, Ser242, Leu265, Asn279	Lys204, Cys205, Gly207, His206, Ser209, Ser211, Glu213, Leu208	Cys329, Lys326, Gln340, Tyr336	Ala238, Glu240, Glu108, Arg85, Arg84, Arg107, Tyr273
9	NTSFR	Glu216, Leu211, Lys154, Gly150, Trp153	Trp27, Glu133, Asp226, Arg228	Arg88, Arg110, Glu243, Asn279, Ser280, Asp239, Ile307	Ser209, Glu213, Ser211, Cys212	Tyr336, Ser338, Val337, Arg328, Cys329	Ser239, Met241, Arg259, Arg247, Glu240
10	YLVAH	Lys154, His209, Leu211, Glu216	Ser28, Glu133, Arg228	Arg110, Tur276, Ser280, Glu243, Ser242, Asn279, Leu265	Glu213, His206, Gly207, Leu208	Phe331, Cys329, Tyr336, Cys339, Lys326	Arg85, Glu108, Arg107, Tyr273, Glu240, Ser277, Phe280
11	YLVAHNLL	Lys154, Lys207, His209, Ser212	Glu71, Arg129, His335, Trp336	Arg88, Arg110, Tyr276, Glu243, Glu247, Val245	Glu213, Ser209, Leu208, Lys204, Cys205	Cys335, Tyr336, Val337, Cys339, Val330	Ala238, Ser237, Lys234, Tyr273, Arg84, Arg107, Phe280, Arg85

The peptide-4 make main interactions with Lys207, Gly210, Ser212, Ser214 and Glu216 of thumb domain, Cys332, Phe334, Trp336, Ser341 and Gln343 of finger domain and Ser242, Glu243, Pro268, Asn279 and Ser280 of palm domain of Wnt3 (Figure 8) (see also Figure S3, Supplementary Materials). Similarly, peptide-11 forms main interactions with Lys204, Cys205, Leu208, Ser209 and Glu213 of thumb domain, Val330, Cys335, Tyr336, Val337 and Cys339 of finger domain and Lys234, Ser237, Ala238, Tyr273 and Phe280 of palm domain of Wnt3a (Figure 9) (see also Figure S4, Supplementary Materials).

Figure 8. 2D-ligand interactions diagram of Wnt3 with peptide-4 derived by IFD docking. (a) Top-docking pose against thumb domain residues, (b) top-docking pose against finger domain and (c) top-docking pose against a patch of palm domain.

Figure 9. 2D-ligand interactions diagram of Wnt3a with peptide-11 derived by IFD docking (a) top-docking pose against thumb domain residues, (b) top-docking pose against finger domain and (c) top-docking pose against a patch of palm domain.

Discussion

Klotho is an important protein known for its role in enhancing cognition, suppressing age like symptoms and extending life span47. It is seen that Klotho deficiency causes chronic renal failure, which is one of the factors underlying degenerative processes such as osteoporosis, skin atrophy, arteriosclerosis aging and bone loss48. Elevating the amount of Klotho in mice has enhanced its cognition and life span. The experiments performed on Klotho deficient mice revealed a decrease in number of stem cell and progenitor cell death7^,49^,50. The Wnt proteins are stem cell regulators known for their role in embryonic development including body axis patterning, bone formation, cell proliferation, cell migration and cell fate specification51^–53. Klotho KL1 domain is found to act as antagonist of Wnt protein. The KL1 domain is 285 amino acids long which interact with Wnt proteins. Analyzing the specific interaction points between Wnt and Klotho proteins can provide significant information about the importance of related residues, and it can enable us to develop peptide inhibitor drugs derived from Klotho as Wnt antagonists.

Due to the nature of Wnt proteins and unavailability of resources, it is difficult to solve 3D structures of whole proteins. Therefore, computational based studies can be applied to address this problem and analyze the structural features. For the comparison studies of our predicted model protein structures with the template structure (xWnt8: 4F0A), we have run 20 ns MD simulations on xWnt8 protein. The secondary structure elements are important for the stability of protein and its domains54. Contrary to the β-sheets, the loop structure is very flexible and makes the CTD including finger domain (β-strand hairpin loop) turn 180° to achieve new stable conformations, which are very different from the initial stage conformations (before MD simulations) as seen in model Wnt3-M0 and Wnt3a-M0 experiment (Figure S1, Supplementary Materials, and Animation-1). Interestingly, the finger domain has taken 180° shift from upper part of its Cysteine knot which also has α-helix. This further supports the importance of secondary structure elements in protein function as well as it gives the indications about some other factors that play important role in the structure stability. The Wnt protein family members have 24–25 Cysteine residues in each protein. In conjunction with the importance of secondary structure elements, disulfide bonds between the Cysteine residues could play important role in stability of a protein structure. Although, the Cysteine residues are scattered over whole protein, the Wnt family shares the same pattern of Cysteine residues hotspots throughout the structure. These hotspots include index finger like domain and Cysteine knot, from CTD, thumb domain and palm domain as well as amino terminal from NTD (Figure 1). Cysteine residues present in the structure help to stabilize the protein, and mutation in any of Cysteine residues results in ectopic intermolecular disulfide bonding which abolishes activity45. To get the reasonable protein model with correct number of disulfide bonds we ran in silico experiments on the homology models with 2 and 5 disulfide bonded models of Wnt3 and Wnt3a, respectively, generated automatically by the SWISS-MODEL, I-TASSER and Prime. The protein structures of Wnt3 and Wnt3a behave entirely different under the same used MD protocols. In Wnt3, the whole protein fluctuates tremendously giving a final structure where the thumb and finger domains get farther from the palm domain. Contrary to Wnt3, the index finger in CTD and thumb domains of Wnt3a comes closer to each other forming a Vitarka mudrā like gesture. This conformation, however, does not correlate with the reference structure of xWnt8. Therefore, another experiment with 9 disulfide bonded models based on the sequence in UniProt was run for both proteins (Wnt3-M2 and Wnt3a-M2). The Cysteine knot domain moved swiftly, and as a consequence, the finger domain moved sideways away from palm domain, leading a ∼90° shift from its initial position in Wnt3-M2. Interestingly, the Wnt3a CTD follows the same pattern of Cysteine knot movement as seen in Wnt3 with nine disulfide bonds, likely due to the unfavorable disulfide bond formation. However, in contrast to Wnt3-M2, the finger domain in Wnt3a-M3 does not show ∼90° movement sideways. The UniProt assignment of nine disulfide bonds in Wnt3 and Wnt3a is actually based on manual assertion inferred from their sequence similarity with WNT8_XENLA, which is also known as xWnt8. The similarity-based approach creates a disulfide bond between Cysteine residues numbered, 300 and 315 in Wnt3 and 297 and 312 in Wnt3a as represented in UniProt. However, MacDonald et al.45 have reported disulfide bonds between 297 and 307 as well as 281 and 312 in the case of Wnt3a. It must be noted that the sequence and number of disulfide bonds reported by UniProt have missing disulfide bond in the base of CTD (Cysteine knot) of Wnt3 and Wnt3a proteins as per our analysis. Furthermore, in contrast to the UniProt annotation regarding disulfide bonds, our study suggests disulfide bonds between 284 and 315, 300 and 310 in Wnt3, similarly, 297 and 307 as well as 281 and 312 in case of Wnt3a. The unstable Cysteine knot could be the reason behind ectopic movement of the CTD. Therefore, the results from the experiment with nine disulfide bonds were again not very reasonable as concluded from the RMSD and RMSF analysis as well (Figures 4 and 6). Ultimately, we carried out another in silico experiment with 12 disulfide bonds following the disulfide bonding sequence between the Cysteine residues as present in xWnt845. In this case, the CTD of Wnt3 (Wnt3-M3) moved towards palm domain with lower RMSD values (average 4.3 Å) as compared to the other two models (Wnt3-M1 and Wnt3-M2) (Figure 3, red cartoon). The model Wnt3a-M3 creates a slight curve in CTD resembling xWnt8 giving precise donut shaped pose of overall protein structure. The protein structure alignment of xWnt8, Wnt3 and Wnt3a having 12 disulfide bonds have also been performed (Figure S5, Supplementary Materials).

The MD simulations for Wnt3 and Wnt3a model with 12 disulfide bonds have also been extended to 30 ns. The Wnt3 protein shows the stable RMSD for MD simulation span from 20 to 30 ns (Figure S6, Supplementary Materials). For the Wnt3a protein model, the graph shows that the RMSD from 20 to 23 ns reaches as high as 9 Å. However, the MD simulation span starting 23–30 ns shows decline and stable in RMSD (Figure S7, Supplementary Materials).

The homology model of Klotho KL1 domain has been generated using cytosolic β-glucosidase as template (PDB ID: 2E9M with 46% identity). The generated model was further directed to energy minimization and MD simulations (Figure S8, Supplementary Materials). Peptide docking has given promising results for the selected peptides (Table 2). For the validation purpose, we have run additional in silico experiments by mutating the peptides and important binding site residues of Wnt protein (Alanin mutagenesis studies) using Glide Standard Precision (SP) docking approach. The residues in all peptides are mutated by adding Alanine one by one in each peptide, totaling 71 individual mutated peptides. These mutated peptides were directed to peptide docking against all three important binding sites of Wnt proteins (Table S1, Supplementary Materials). The peptide docking results showed that the mutation in the peptides, in general, does not give better docking scores compared to the wild type peptides, as expected. For instance, the peptide-1 has the docking score of −8.1 kcal/mol against Wnt3 thumb domain. Mutating from first residue (Tyrosine to Alanine) in peptide-1 was lowering the docking scores to −6.0 kcal/mol with the difference of 2.1 kcal/mol. Similarly, peptide-2 is having the docking score of −8.1 kcal/mol against Wnt3 palm domain. Mutating the first residue (Tyrosine to Alanine) decreases the results to −5.1 kcal/mol with the difference of 3.0 kcal/mol. There are more examples (peptide 3, 7 and 10) that have decrease in docking scores after mutation in all the domains of Wnt3 and Wnt3a, as seen in Table S1, Supplementary Materials. However, some of the peptide mutations resulted with improved docking scores compared to the wild type peptides, for instance, peptide-2 against finger domain Wnt3a. Most prominently, six different combinations of peptide-2 with Alanine mutations played more actively in peptide docking against Wnt3a finger domain. The wild type peptide-2 gives −4.4 kcal/mol docking score against the same binding pocket. Mutating the third residue Glycine to Alanine leads docking score to −8.1 kcal/mol and improved the binding scores by 3.7 kcal/mol. Similarly, mutated peptide-2 combination 2.5 and 2.6 improved the results by 2.9 and 3.5 kcal/mol against same binding pocket (Table S1, Supplementary Materials). Peptide-8 mutations also have shown slightly better results against finger domain binding site of Wnt3a compared to the wild type peptide. For example, the peptide 8.6 with Alanine mutation in place of Aspartic acid (residue number 6) gives −7.8 kcal/mol docking score compared to the −4.0 kcal/mol as wild type peptide docking score. Although, the mutated peptide-8 combinations are active against two domains (finger and thumb domains), yet these combinations are not found to be as active as wild type peptides against other domains of Wnt3 and Wnt3a. Except some of the mentioned mutated peptides in the detail, wild type peptides are overall potentially active against Wnt3 and Wnt3a, even the single residue mutation is affecting and decreasing predicted binding affinity of the peptides against binding sites of Wnt3 and Wnt3a.

In previous section, Alanine scanning mutagenesis studies have been performed for the Klotho-derived peptides. In the next in silico experiment, we have mutated one important binding site residue to Alanine in all three important binding sites of Wnt3 and Wnt3a proteins. Selection of critically important residues for Wnt3 and Wnt3a is based on sequence similarity to the Lys172 and Trp319 of xWnt8 for thumb and finger domain, respectively14. In case of Wnt3, we have mutated following amino acids Lys207Ala, Trp336Ala and Glu243Ala. Likewise, the residues Lys204, Trp333 and Glu240 have been mutated to Ala in Wnt3a (Table S2, Supplementary Materials). The Glu243 and Glu240 from Wnt3 and Wnt3a have been selected based on its active participation in making interaction with peptides in our main experiment. The noteworthy decline in the peptide docking score in this study after mutating the critical amino acid in binding site indicates the strong interaction and specificity of wild type peptides against the critical amino acids in the binding site of Wnt3 and Wnt3a. As Table 2 shows, we can see that the wild type peptides 4, 6, 7, and 11 are interacting with the Trp336 and its surrounding residues. However, most prominently, the peptide-3 does not show specificity against the Trp336 and gave improved docking score of −8.2 kcal/mol upon mutation.

Moreover, in order to test the peptide screening procedure, random peptides with similar length were derived. For this aim, we used our in house developed peptide generator script (Figure S9, Supplementary Materials). Randomly and automatically picked peptides were then docked at the Wnt3 and Wnt3a target sites and their docking results are compared with the proposed peptides. Results showed that proposed peptides with similar length usually have better docking scores than randomly picked peptides (Table S3, Supplementary Materials).

The peptide docking of wild type peptides against protein xWnt8, which is used as template for modeling of Wnt3 and Wnt3a, has also been carried out (Table S4, Supplementary Materials). The docking results have shown the effects of peptides against all three important binding sites of xWnt8 nearly as good as in case of Wnt3 and Wnt3a. Although the peptides 4, 8 and 11 could not provide any important docking pose against co-receptor binding site at palm domain, the peptides 5, 7, 9 and 10 give prominently good docking scores in all the three binding sites (Table S3, Supplementary Materials).

Conclusions

All carried in silico experiments regarding the identification of appropriate conformation of secondary structure elements and the correct number of disulfide bonds in Wnt3 and Wnt3a provided a detailed insight into the versatility of Wnt family proteins. Despite the high sequence identity between Wnt3 and Wnt3a proteins (i.e. 84.2%), these proteins have significantly different structure and dynamics under the same conditions (i.e. same MD protocols) indicating that the 15.8% difference in sequence identity has a considerable impact on protein structure as well as dynamics55–60. The behavior of protein domains towards the number of disulfide bonds is an artifact as seen in the experiments performed in our study. In case of Wnt3, lesser the number of disulfide bonds, afar the finger and thumb domains move. In contrast, lesser the number of disulfide bonds in Wnt3a structure, closer the distance between the finger and thumb domains.

These secreted Cysteine rich lipoproteins (Wnt proteins) bind to Frizzled (FZD) receptor and LDL receptor-related protein 6 (LRP6). As described in solved structure of xWnt8, the index finger domain/CTD and thumb domain/NTD engage FZD. All the Wnt proteins undergo characteristic glycation, acylation, and other modifications leading to insolubility. One of such example is the attachment of mono-unsaturated fatty acid (palmitoleic acid) at conserved Serine, which interacts to a pocket in the Frizzled Cysteine rich domain, showing conformation that emphasizes the importance of the lipid for signaling. The lipid modification on the Wnts is required for active signaling, and perhaps important for Wnt secretion. Predicted location of palmitoylation after multiple sequence alignment of Wnt3 and Wnt3a is Ser 212 and Ser 209, respectively. Addition of palmitoleic acid at NTD leads to an important domain for peptide binding.

The protein–protein docking analysis of KL1 domain and Wnt proteins provided insight into the important residues that could form good interactions with the important domains of Wnt3 and Wnt3a. These important residues were chopped out from the protein as a string of adjacent amino acids, and used to prepare peptides employing in silico based techniques. Out of 11 peptides, three were specific peptides that were extracted by DME, which can perform glycosidases activity, as predicted by DME. These peptides were further validated by putting them into peptide docking experiments against the three domains of Wnt3 and Wnt3a (Finger domain, thumb domain with palmitoleic acid attached, a patch of ∼10 residues at the solvent exposed site of palm domain) and results showed excellent docking scores against most of these important domains as listed in Table 2. Validation of these peptides has been carried out in silico by applying mutations in peptides and Wnt protein binding sites. Result of the in silico validation experiments gives satisfactory results about the predicted activities of the peptides. These Klotho-derived peptides can serve as potential antagonist against Wnt proteins, however, this still needs to be validated by wet-lab experiments.

Disclosure statement

The authors report no declarations of interest.

Funding

This study is supported by The Scientific and Technological Research Council of Turkey (TÜBİTAK) under TUBITAK-2216 program (Research Fellowship Program for International Researcher). PhD student SBM is supported by this program.