Structural basis of βTrCP1-associated GLI3 processing

Controlled ubiquitin-mediated protein degradation is essential for various cellular processes. GLI family regulates the transcriptional events of the sonic hedgehog pathway genes that are implicated in almost one fourth of human tumors. GLI3 phosphorylation by Ser/Thr kinases is a primary factor for their transcriptional activity that incurs the formation of both GLI3 repressor and activator forms. GLI3 processing is triggered in an ubiquitin-dependent manner via SCFβTrCP1 complex; however, structural characterization, mode of action based on sequence of phosphorylation signatures and induced conformational readjustments remain elusive. Here, through structural analysis and molecular dynamics simulation assays, we explored comparative binding pattern of GLI3 phosphopeptides against βTrCP1. A comprehensive and thorough analysis demarcated GLI3 presence in the binding cleft shared by inter-bladed binding grooves of β-propeller. Our results revealed the involvement of all seven WD40 repeats of βTrCP1 in GLI3 interaction. Conversely, GLI3 phosphorylation pattern at primary protein kinase A (PKA) sites and secondary casein kinase 1 (CK1) or glycogen synthase kinase 3 (GSK3) sites was carefully evaluated. Our results indicated that GLI3 processing depends on the 19 phosphorylation sites (849, 852, 855, 856, 860, 861, 864, 865, 868, 872, 873, 876, 877, 880, 899, 903, 906, 907 and 910 positions) by a cascade of PKA, GSK3β and CSKI kinases. The presence of a sequential phosphorylation in the binding induction of GLI3 and βTrCP1 may be a hallmark to authenticate GLI3 processing. We speculate that mechanistic information of the individual residual contributions through structure-guided approaches may be pivotal for the rational design of specific and more potent inhibitors against activated GLI3 with a special emphasis on the anticancer activity.

. Schematic illustration of GLI3 translocation, processing, and degradation via SCF (SKP1, Cullin, F-box protein complex). In the absence of hedgehog (right panel), Ptch1 constitutively inhibits Smo, preventing its ciliary localization. In this state, GLI proteins are retained in a complex with Sufu, where PKA, CK1ε, and GSK3 phosphorylate them. Phosphorylated GLI3 binds the SCF complex that is partially ubiquitinated and processed by the proteasome into GLI3-R, which translocates to the nucleus and represses transcription of Hh target genes, including Akt, Gli1, Ptch1 etc., prior to degradation by an unknown E3 ligase complex. In the presence of Hh (Left panel), PKA phosphorylation is inhibited. The activated form of GLI3 not only prevents it from processing but also permits its subsequent transport to the nucleus to allow activation of transcription of Hh target genes. Following this, GLI3-A is ubiquitinated by the SPOP/CUL3 complex and degraded by the proteasome. P, phosphate group; Ptch1, Patched 1; SUFU, Suppressor of Fused; SMO, Smoothened; PKA, Protein kinase A, GSK3, Glycogen synthase kinase 3; CSKI, Casein kinase 1; CUL1, Cullin 1; CUL3, Cullin 3; RBX1, RING-box protein 1; SKP1, S-Phase Kinase-Associated Protein 1; SPOP, Speckle Type BTB/POZ Protein; GLI3-R, GLI3 repressor; GLI3-A, GLI3 activator.

Molecular docking analysis.
In order to explore the binding pattern of βTrCP1 and GLI3 peptides (GLI3-un, GLI3-β1, GLI3-β2, GLI3-β3, GLI3-β4 and GLI3-β1-4), molecular docking was accomplished by HADDOCK 26 . Prior to docking analysis, 3D structures of βTrCP1 and GLI3 peptides were submitted to CPORT server to predict the interface residues 27 . CPORT provides a list of active and passive residues of the interaction interfaces that are further employed for the preparation of ambiguous interaction restraints (AIR) by HADDOCK 28,29 . Docking simulations were carried out with default parameters, among them HIS protonation state, random removal of restraints, number of structures to dock and refine, electrostatic scaling constant, restraint energy constants, scoring parameters, temperature and time steps of refinement processes were defined automatically by the interface of HADDOCK server 30 . HADDOCK produces 1000 docked structures through experimental data to drive docking. HADDOCK scores each model using Equation 1 and retains the top 200 solutions for consequent flexible refinement, where E AIR , E elec , E vdW and E desolv are the restraints violation , electrostatic, van der Waals and desolvation energies, respectively. BSA is the buried surface area and E data indicates the energy of other restricted data. HADDOCK score is weighted as the sum of the following four terms: electrostatic energy (weight: 0.2), Van der Waals energy (weight: 1.0), desolvation energy (weight: 1.0) and restraint violation energy (weight: 0.1). Furthermore, selected models are subjected to a semiflexible refinement followed by water refinement step in torsion angle space and explicit water shell, respectively. These parameters are scored using Equations 2 and 3, respectively. The solutions are clustered using a 7.5 Å limit created on their pairwise ligand interface. Root mean square deviation (RMSD) values and cluster ranks are rendered to the average score of the top four structures in each cluster. The best docked complexes of top ranked clusters were selected and visualized by UCSF Chimera 1.11.2 to analyze their conformational readjustments. Residual interactions such as hydrogen bonds, hydrophobic, electrostatic interactions and bond lengths were evaluated by employing DIMplot embedded in LigPlus 31 . Furthermore, these complexes were subjected to molecular dynamic simulations for detailed analysis.
Molecular dynamics simulation analysis. In order to gain a deep insight into the mechanism of GLI3 phosphodegron recognition by apo-βTrCP1 and its bound forms with GLI3-β1, GLI3-β2, GLI3-β3, GLI3-β4 and GLI3-β1-4, these complexes were subjected to molecular dynamics (MD) simulation assays to evaluate the stability, folding, conformational changes and dynamic behavior of interacting proteins. Groningen Machine for Chemicals Simulations (GROMACS) 5.1.4 package was used to perform MD simulation assay. All MD simulations were performed by GROMOS96 43a1 force field 32 to acquire the equilibrated system. All systems were www.nature.com/scientificreports www.nature.com/scientificreports/ solvated using SPC water model 33 in a periodic box, followed by the addition of Na + and Cl − counter ions to neutralize the system. Before MD simulations, systems were subjected to energy minimization to remove initial steric clashes using 1000 steps of steepest descent algorithm via a tolerance of 10 KJ/mol/nm −1 . Systems were equilibrated for 1000 ps at 300 K and 1 bar pressure in NVT 34 and NPT 35 ensembles, respectively and their equilibration profiles were evaluated (Fig. S1). The hydrogen bond length constraints were employed at a time step of 2 fs for numerical integration with Verlet leap-frog algorithm 36 . PME (Particle Mesh Ewald) algorithm 37 was employed to evaluate the long-range (LR) electrostatic interactions. Finally, MD simulations were run for 40 ns time scale. Stability and time-dependent behavior of each system were investigated using system conformations extracted every 5 ns from the MD trajectories and analyzed using UCSF Chimera and GROMACS tools. GROMACS modules such as g_rms, g_rmsf, g_energy and g_hbond were utilized to analyze the stability and behavior of each system. The secondary structure database (DSSP) was utilized to analyze the time-dependent secondary structure fluctuations in the bound complexes 38 .
Binding free energy calculation. Poisson-Boltzmann or generalized Born and surface area continuum solvation (MM/PBSA) method 39 was employed to calculate the binding free energy of the system. This method provides inclusive energy composition and improves docking energy via incorporating protein flexibility. The binding energy of ligand-protein complex was calculated using the following equation: binding c omplex protein l igand G complex is the total free energy of the protein-ligand complex; G protein and G ligand are total energy of separated protein and ligand in solvent, respectively. The free energy for each individual G complex, G protein and G ligand were estimated by: x is the protein-ligand complex. E MM is the molecular mechanics energy and G solvation is free energy of solvation. The molecular mechanics potential energy was calculated in vacuum as following: E bonded is bonded interaction including bond, angle, dihedral and improper interactions and E non-bonded is non-bonded interactions that consist of van der Waals (E vdw ) and electrostatic (E elec ) interactions.The solvation free energy (G solvation ) was estimated as the sum of electrostatic solvation free energy (G polar ) and apolar solvation free energy (G non-polar ): solvation p olar non polar G polar was computed using the Poisson-Boltzmann (PB) equation 40 and G non-polar was calculated using a solvent accessible surface area (SASA) as follows: non polar γ is a coefficient related to surface tension of the solvent and b is fitting parameter.

Results
Structural evaluation of GLI3 peptides. The predicted structures of GLI3 peptides (GLI3-β1, GLI3-β2, GLI3-β3, GLI3-β4 and GLI3-β1-4) were evaluated by Ramachandran plots (Fig. S2), where blue colour indicated favorable region (sterically allowed regions), while no outliers were observed. Approximately, 92-95% residues were resided in the blue region. Additionally, parameters like peptide bond planarity, non-bonded interactions, Cα-tetrahedral distortion, main chain H-bond energy values and overall G-factors for the predicted models were lying in the favourable ranges. GLI3 peptide structures optimized through GROMACS tool were further evaluated by RAMPAGE 41 .

Phosphopeptide binding and conformational transitions.
In order to evaluate mechanism of substrate recognition by βTrCP1, GLI3 phosphopeptides were subjected to molecular docking analysis. Given a max- .0%, 65.0%, 57.5%, 46.5%, 58.5% and 57.0% of water-refined models, respectively. The statistics of top 10 clusters (ranked on the basis of lowest overall energy and Z-score values) were shown by HADDOCK, out of which scores of the optimal clusters for each βTrCP1-peptide complexes are illustrated in Table 1. The more negative HADDOCK and Z-scores indicate a reliable interaction. Z-score is the quantitative measure of cluster standard from the average score. All βTrCP1-peptide complexes were carefully characterized to access their binding patterns. In case of βTrCP1-GLI3 PKA complex, phosphopeptide exhibited binding with the 1st, 2nd and 7th WD40 repeats of βTrCP1 having a score of −17.6 ( Fig. 2A). In contrast, GLI3 GSK3β and GLI3 CKIϵ peptides did not exhibit binding with βTrCP1 ( Fig. 2B,C). In βTrCP1 and GLI3-β1-4 complex, phosphopeptide binding was observed at the upper (2019) 9:6865 | https://doi.org/10.1038/s41598-019-43392-3 www.nature.com/scientificreports www.nature.com/scientificreports/ interface of β-propeller (Fig. 2D). Thus GLI3 phosphorylation by all three enzymes (PKA, GSK3β and CKIε) resulted in accurate binding with βTrCP1 substrate binding site.
Furthermore, PDB files were characterized to measure the conformational switches in GLI3 phosphopeptides upon binding to βTrCP1. All phosphopeptides exhibited quite stable binding patterns at 25 ns. Particularly, upon binding to βTrCP1, both helical regions (Ile854-Ser864 and Thr900-Glu908) were shortened in GLI3-β1-4 to accommodate it in the cavity formed by β-propellers (Fig. 7E). A profound conformational change was observed in Thr900-Glu908 helical region (Fig. 7D), as upon binding to βTrCP1, this helix was completely missing. This trend was observed throughout MD simulation run as evident from the analysis of time-dependent secondary structure fluctuations via DSSP module (Fig. S3). Another notable secondary structural amendment was witnessed in the loop region of GLI3, where Sep875-Glu878 region of GLI3-β2 adopted a α-helical conformation upon binding to βTrCP1 (Figs S3B and  7B). In βTrCP1-bound GLI3-β1, GLI3-β3, GLI3-β4 and GLI3-β1-4 peptides, this region remained structurally preserved (Fig. S3). Subsequent analysis of RMSF indicated residual flexibility of phosphorylated residues upon GLI3 binding to βTrCP1. In case of GLI3-β1 and GLI3-β3 binding, major fluctuations up to 10 Å and 4.5 Å were perceived in all phosphorylated residues (Fig. 7F). Correspondingly, GLI3-β2 and GLI3-β4 peptides exhibited minor rate (up to 2.8 Å) of fluctuations as compared to other simulated systems. In case of GLI3-β1-4, significant fluctuations were detected in Sep899, Sep903, Sep906, Sep907 and Sep910 residues (4-11 Å) to assist in binding, while phosphorylated residues involved in binding (Sep852, Sep855, Sep872, Sep873, Sep876, Sep877 and Sep880) were quite stable (Fig. 7E). These results specified that Sep899-Sep910 of GLI3-β1-4 exhibited more fluctuations thus suggesting that Sep899-Sep910 region of GLI3 may be crucial for βTrCP1 binding.  Table 2. Binding residues of βTrCP1 and GLI3 phosphopeptides. Residues involved in hydrogen bonding and hydrophobic associations are indicated in bold and normal forms, respectively.

Discussion
Ubiquitination plays a crucial role in Hh signaling activity via GLI proteins 2 that act as zinc finger transcriptional effectors and regulate vertebrate development 18 . In the absence of Hh, phosphorylation-mediated ubiquitination keeps GLI3 in the repressor form. Despite the fact that mechanistic role of phosphorylation-mediated GLI3 degradation is well-established, the structural-functional paradigm is largely unknown. Here, we explored the potential role of multisite phosphorylation in the selective binding of βTrCP1 and GLI3 phosphopeptides via in silico approaches. Evidently, key substrate binding residues (Arg285, Arg307, Ser327, Ile344, Cys347, Lys365, Arg367, Asn394, Ser448, Arg474 and Arg521) of βTrCP1 were consistent with the earlier studies 19,42 . RMSD analysis demonstrated stability (2-4 Å) in all systems at 12 ns. Further analysis elucidated multiple conformational changes that invoked specificities in the β-propeller upon phosphopeptide binding. The overall topology  www.nature.com/scientificreports www.nature.com/scientificreports/ of β-strands remained preserved in the βTrCP1 structure (Table 3). A predominant transformation of β-strand (Thr381-Leu386) into loop conformation facilitated the binding via flexibility. Other prominent positional readjustments observed in the β-strands were localized in WD40 repeat-1 (Arg301-Leu303 and Leu313-Tyr315) and Table 3. Secondary structure changes during MD simulations in phosphopeptide-bound βTrCP1 states with reference to apo-βTrCP1. Figure 6. Structural details of βTrCP1 and GLI3 phosphopeptide binding. βTrCP1 is represented by light gray ribbon, while pale yellow ribbons represent phosphopeptide GLI3-β1-4 with interacting residues indicated by coral ball and stick mode. Illustration of four sequence motifs (β1 to β4) related to the βTrCP1 binding site are underlined that are phosphorylated by a putative cascade of PKA, GSK3β and CK1. PKA phosphorylated serines (phosphoserine) in the sequence motifs are colored in red. GSK3β phosphorylates serines (green) four residues N-terminal to a phosphoserine, while CK1 phosphorylates serines (blue) three residues C-terminal to a phosphoserine; both can chronologically multiphosphorylate GLI3 after priming. Middle panel shows the conservation pattern of βTrCP1 binding residues upon phosphopeptide binding. X-axis indicates the binding residues of βTrCP1 and Y-axis indicates the GLI3 phosphopeptides (GLI3-β1, GLI3-β2, GLI3-β3, GLI3-β4 and GLI3-β1-4). Dot represents the contribution of respective residue in binding to phosphopeptide. www.nature.com/scientificreports www.nature.com/scientificreports/ repeat-6 (Ile492-Trp495 and Ile532-Ser534), leading to GLI3-β1-4 binding. In RMSF analysis, βTrCP1 binding region (Arg285-Arg521) attained more stability upon binding to GLI3-β-4 (Fig. 4D).
In agreement to the previous observations 17 , where crucial role of GLI3 motif-4 has been reported in βTrCP1 binding and GLI3 processing, our findings indicate that motif-4-specific Ser899 phosphorylation invokes other phosphorylated serines to impart active role in binding to βTrCP1 (Fig. 8). The interaction of βTrCP1 and GLI3 was significantly influenced by the positional readjustments of residues lying in two helices (Ile854-Sep865 and Thr900-Sep907) due to phosphorylation of paired neighboring residues that induced flexibility differences through helix-loop inter-conversion. Generally, introduction of phosphate group targets loop conformation by rearranging the hydrogen bonding network of side chains lying at the vicinity of loop region 48 . These transitions in the surrounding regions render helical shifting into loop that acts as a conformational switch for the binding cleft geometry 49 . The presence of diverse hydrogen bonding pattern and conformational switching due to phosphorylation is crucial for the recognition of GLI3 by βTrCP1. Any change in this pattern may impair their binding affinity due to imbalanced phosphorylation level. It is however unclear at the moment how energy barrier overcomes the phosphorylation or other post-translational modification-induced conformational space.
GLI3 contains multiple binding sites for βTrCP1, where approximately, two-third of GLI3 contacts involve phosphorylated PKA sites and secondary CK1/GSK3 sites. The potential involvements of GLI3-specific primary (Ser852, Ser873, Ser877, Ser880 and Ser903) and secondary phosphorylated (Ser855, Ser872, Ser876, Ser906) residues in βTrCP1 binding indicate that both primary and secondary phosphorylations are required for βTrCP1 binding. Study of interdependent phosphorylation status through structural knowledge may expand the www.nature.com/scientificreports www.nature.com/scientificreports/ repertoire of GLI3 processing. Indeed, any mutation at the PKA-specific sites may significantly reduce the binding of βTrCP1 to GLI3 7 . Taken together, our results are in good agreement with the experimental outcomes 7,17 . This study may uncover the spectrum of structural linkages in association with the kinase-mediated phosphorylation paradigm to illustrate the molecular basis of GLI3 processing in Hh signaling. Further studies will be needed to elaborate the effect of putative phosphorylation site mutations at structural level.

Data Availability
All data generated or analyzed during this study are included in this published article (and its supplementary information files). Table 4. Free energy (kJ/mol) calculation for βTrCP1 in complex with GLI3 phosphopeptides.