Exploring protein phosphorylation by combining computational approaches and biochemical methods.

Post-translational modifications of proteins expand their functional diversity, regulating the response of cells to a variety of stimuli. Among these modifications, phosphorylation is the most ubiquitous and plays a prominent role in cell signaling. The addition of a phosphate often affects the function of a protein by altering its structure and dynamics. However, these alterations are often difficult to study and the functional and structural implications remain unresolved. New approaches are emerging to overcome common obstacles related to the production and manipulation of these samples. Here, we summarize the available methods for phosphoprotein purification and phosphomimetic engineering, highlighting the advantages and disadvantages of each. We propose a general workflow for protein phosphorylation analysis combining computational and biochemical approaches, building on recent advances that enable user-friendly and easy-to-access Molecular Dynamics simulations. We hope this innovative workflow will inform the best experimental approach to explore such post-translational modifications. We have applied this workflow to two different human protein models: the hemeprotein cytochrome c and the RNA binding protein HuR. Our results illustrate the usefulness of Molecular Dynamics as a decision-making tool to design the most appropriate phosphomimetic variant.

Introduction

The complexity of the human proteomes goes far beyond the ca. 20,000 ORFs (open reading frames) contained in the genome. A myriad of processes give rise to the structural and functional variety of proteins without increasing the genetic background. Alternative promoters and splicing, as well as mRNA editing yield a five-fold increase in the transcriptome. Post-translational modifications (PTMs) represent another tool by which cells can expand the functional diversity of proteins. These modifications may occur at any step in the “life cycle” of a protein and alter either its function or stability. We can sort PTMs into five groups: 1) addition of complex molecules, 2) addition of proteins and peptides, 3) amino acid modifications, 4) addition of chemical groups and 5) cleavage of part of the protein sequence [1]. This work will focus on phosphorylation, which is the main PTM involved in cell regulatory processes [2], [3], [4], [5], [6].

Metabolic regulation by protein phosphorylation is notable for its rapid reversibility, an essential feature in cell signaling [7]. This PTM consists of the addition of a phosphoryl group either to a hydroxyl-containing residue through an ester bond, or to a nitrogen-containing one to from an amidate, this process is generally mediated by kinases. The energy of this phosphate ester linkage is similar to that of the ATP β–γ anhydride bond. In fact, under certain experimental conditions, it is possible for kinases to act in the opposite direction, thereby generating ATP. Thus, high ATP:ADP ratios forestall dephosphorylation by protein kinases [6], [8], [9]. The residues most prone to phosphorylation are Ser (ca. 60%), followed by Thr (ca. 25%) and Tyr (ca. 15%) [6], [10]. Phosphorylation of other residues—such as histidine and arginine—is less stable and infrequent in eukaryotes [11], [12], and is thus excluded from this work.

Phosphorylation sites can be found within intrinsically disordered regions (IDRs) or structured domains. Analyses suggest that serine and threonine phosphorylation takes place predominantly in IDRs, while tyrosine phosphorylation occurs both in IDRs and structured regions [13]. Phosphorylation induces both disorder to order and order to disorder transitions. An example of order to disorder transition is nucleophosmin (NPM1), sequential phosphorylation of solvent exposed Ser and Thr regulates the thermodynamic stability of the NPM1 oligomerization domain. This phosphorylation event promotes the exposure of other sites buried within the structure for subsequent phosphorylation and drives the switch from structured pentamer to disordered monomer [14]. K-homology splicing regulator protein (KSRP) contains four K-homology domains (KH1-4), in which KH1 phosphorylation acts as a conformational switch between folded (non-phosphorylated protein) and essentially unfolded (phosphorylated protein) states [15]. Unfolding of KH1 domain creates a 14-3-3ζ binding site which prevents the recruitment of the exosome [16]. On the other hand, phosphorylation of IDR eukaryotic translation initiation factor 4E (eIF4E)-binding protein 2 (4E-BP2) protein constitutes an example of disorder to order transition by inducing folding into a four-stranded β-domain, which blocks the accessibility of eIF4E [17].

Overview of phosphoprotein detection and purification methods

Phosphoproteomic investigations usually require purification and often enrichment of endogenous proteins from cells or tissues. The substitution of the hydroxyl group of Ser, Thr and/or Tyr by a phosphate group affects the physicochemical and immunogenic properties of any protein, and is consequently used for its isolation and analysis. Purification methods for phosphorylated proteins usually take advantage of 1) changes in the isoelectric point (pI), 2) metal cation chelating capacity of phosphate groups, and 3) unique immunogenicity of isoforms [18].

The main purification techniques used for isolation of phosphorylated proteins are chromatofocusing (CF), ion exchange chromatography (IEX), immobilized metal ion affinity chromatography (IMAC), metal oxide affinity chromatography (MOAC) and immunoaffinity chromatography (IAC). CF and IEX are based on the downshift in protein pI caused by phosphorylation [19]. The magnitude of this shift depends strongly on amino acid composition, secondary structure and number of concurrently-phosphorylated residues [18], [19]. There are several pI prediction algorithms which calculate the influence of different PTMs, including phosphorylation, but such theoretical approaches are not exempt from errors and their accuracy for phosphorylated proteins depends on neighboring phosphoresidues and protein folding [20], [21].

Many phosphoprotein purification strategies utilize the tendency of phosphate groups to chelate metal ions and oxides. IMAC consists of phosphoproteins binding to metal-functionalized resins, usually containing Fe3+, and its elution by either a pH gradient or competitive displacement by chelating ligands, typically imidazole [22], [23]. To improve the selectivity and efficacy of IMAC, Kinoshita et al. [24] designed a phosphate chelating group, the so-called Phos-tag, which is an alkoxide-bridged binuclear metal (Zn2+ or Mn2+) complex. Phos-tag and its conjugated derivatives bind to phosphorylated residues in an amino acid-independent manner [25]. This technology has been extensively used for phosphoprotein detection by SDS-PAGE (sodium dodecyl sulfate–polyacrylamide gel electrophoresis) and western blot analyses [25], [26], [27], and has not yet been surpassed by other types of IMAC developed over the last few years [22], [23], [28]. As an alternative to this technique, MOAC employs metal oxides, most commonly TiO2, to selectively isolate phosphoproteins [29]. As with IMAC, phosphoproteins are complexed with metal oxides under acidic conditions; however, complete recovery of phosphorylated proteins in MOAC requires an alkaline pH gradient from 8.5 to 11.5 for elution of multiple phosphorylated isoforms [30].

When purifying phosphoproteins, IAC is more selective than IMAC or MOAC, though it still suffers from low specificity. IAC exploits the unique immunogenic properties of phosphorylated proteins, utilizing antibodies that detect specific phosphoresidues in a sequence-independent manner. Antibodies raised against kinase consensus sequences can also be used [31]. Currently, only anti-phosphorylated Tyr (pTyr) antibodies are widely used in phosphoproteomics, despite the lower abundance of phosphorylated Tyr residues, compared with Ser and Thr [32]. The relatively low specificity of anti-phosphorylated Ser (pSer) and anti-phosphorylated Thr (pThr) antibodies favors the use of more cost-efficient approaches such as IMAC/MOAC [31]. In fact, the best asset of IAC, i.e. the high affinity of antigen–antibody interactions, is offset by the high cost of antibody production [18]. Recent papers report combined purification strategies using IEX, TiO2-MOAC and/or IAC, either sequentially or as complementary techniques, to optimize protein enrichment into phosphotyrosine [33], [34], [35].

Unlike anti-pTyr antibodies, the Src homology 2 (SH2) domain has recently emerged as an affinity reagent in the field of phosphoproteomics. SH2 modular domains are present in many signaling proteins such as growth factor receptor-bound and signal transducers and activators of transcription (STAT) protein family members. Single SH2 domain expression and covalent coupling to chromatographic resins allows the purification of phosphotyrosine proteins [36], [37]. For selectivity improvement, Kaneko et al. [38] designed high-affinity SH2 domains by directed evolution methods. They discovered that the substitution of three key residues in the SH2 binding site of the tyrosine kinase Fyn notably increases its affinity for Tyr-phosphorylated proteins.

Obtaining phosphoproteins from natural sources has limitations in terms of sample quantity and possible downstream applications. Many researchers dealing with endogenous phosphoproteins conduct quantitative and temporal phosphoproteomic profiling of the components of cellular signaling networks under certain stress conditions, including disease [39], [40], [41]. The methodology used in these studies includes chromatographic techniques for phosphoprotein isolation and enrichment so far described, as well as typical mass spectrometry analyses [37], [42]. Nonetheless, natural sources of endogenous phosphoproteins are not optimal since large quantities are needed for structural and/or biophysical assays.

The simplest and most common strategy to produce high amounts of any protein is the use of overexpression vectors for bacterial (usually E. coli) or yeast systems. They allow the expression of the construct of interest, normally a gene encoding a fusion protein with a tag sequence—e.g. polyhistidine, glutathione S-transferase (GSH), or maltose-binding protein—to facilitate purification by affinity chromatography (Ni2+-IMAC, GSH-sepharose or amilose resin). Assuming that the expression and purification of the recombinant protein provides good yields, both in terms of quality and quantity of the sample, it is still necessary to phosphorylate in vitro the target residues using the appropriate kinase enzyme. However, there are very few examples in the literature of preparative methods for in vitro phosphorylation of proteins used in binding assays, i.e. to study the effect of protein phosphorylation on an interaction with a third protein [43], [44]. In most papers, the in vitro kinase reactions are intended to characterize the phosphorylation event itself – to determine parameters such as kinetics, kinase substrate specificity and target-site identification [45], [46], [47], [48].

Baudin and Müller [45], [46], [47], [48] took a different approach and in bacteria co-expressed both a kinase domain and a STAT protein for in vivo phosphorylation. Using this approach, they purified Tyr-phosphorylated STAT proteins, which were used to assess the binding of STAT proteins to specific DNA sequences by electrophoretic mobility shift assay (EMSA). This successful strategy could potentially be applied to other kinase-phosphoprotein pairs. Although, if any given kinase phosphorylates more than one residue, it will be necessary to resort to point (e.g. Tyr to Phe) and/or deletion mutants of the target protein to study the effect of each phosphoresidue separately.

The lack of information on specific kinases for the vast majority of the identified phosphorylation targets demands the adoption of other methodological approaches. Thus, the use of phosphomimetic mutants, either classic variants (e.g. Ser/Thr/Tyr to Asp or Glu) or those incorporating non-canonical amino acids (ncAAs), is a good alternative for studying the structure and interactions of phosphorylated proteins.

Classical phosphomimetic variants using standard amino acids

As discussed above, the use of endogenous phosphorylated protein is frequently not feasible because of its low yield when purified from cell extracts or the absence of information on specific kinases. In this case, the easiest approach is the use of phosphomimetic proteins. Since phosphorylation is a post-translational event, cells lack native tRNAs to carry phospho-residues and to transfer them to nascent polypeptides. Classically, Ser- and Thr-phosphorylation have been mimicked by mutations to aspartate and glutamate. These substitutions provide the negative charge of each phosphoresidue and are of a similar volume. Therefore, they constitute the simplest methodology to introduce canonical amino acids [49], [50]. Glu is a more suitable phosphomimetic replacement as it matches the phosphate oxygens of pSer and pThr, while Asp is less suitable as its chain is shorter. However, a major drawback of these mutations is they fail to feature the hydration layer and formal charges of phosphate [7].

Despite their utility, residue substitutions are not always possible. In some cases, it is necessary to introduce a pair of acidic residues to generate the double charge and mimic the behavior of the natural phosphoprotein [51], [52], [53]. However, this strategy may lead to local stress on protein structure due to side chain repulsion.

Glutamate is often used to mimic pTyr [54], [55]. However, the glutamate residue is unable to mimic either the charge or the volume of pTyr. Indeed, there are examples in the literature wherein Tyr phosphorylation could not be emulated by this mutation [56], [57]. Alternatively, several ncAAs can be used to mimic the abovementioned properties of the phosphate group (Table 1).

Table 1

Non-canonical amino acids mimicking phosphoresidues.

Phosphorylated residues	Unnatural amino acid	Citation
Phosphotyrosine	4′-Carboxymethyloxy-3′-phosphonophenylalanine	[85]
	N-[4-(2-{2-[3-(2-Bromo-acetylamino)-propionylamino]-3-hydroxy-propionylamino}-ethyl)-phenyl]-oxalamic acid	[86]
	3-(4-{2-[2-(2-Bromo-acetylamino)-ethyldisulfanyl]-ethylcarbamoyl}-cyclohexylcarbamoyl)-pyrazine-2-carboxylic acid	[86]
	Photocaged phosphotyrosine	[73]
	Sulfotyrosine	[87]
	p-Carboxymethyl-l-phenylalanine	[88]
	p-(Phosphonoamino)-phenylalanine	[89]
	O-phosphotyrosine	[90]
	Benzyl 2-benzyloxycarbonylamino-3-[4-bis(dimethylamino)phosphoryloxyphenyl] propanoate	[91]
	4-Phosphomethyl-L-phenylalanine	[77]
Phosphothreonine	Photocaged phosphothreonine	[73]
	O-phospho-l-threonine	[92]
Phosphoserine	Photocaged phosphoserine	[73]
	Phosphonomethylene alanine	[161]

Incorporation of unnatural amino acid to obtain phosphomimetic proteins

The use of ncAAs in recombinant protein expression enables modification of the chemical properties of one or more residues [58]. Moreover, the use of technologies to selectively and co-translationally incorporate phosphomimetic amino acids into proteins can increase sample preparation yields to enable biophysical analysis. There are two current options for introducing ncAAs into polypeptides: genetic code expansion and reprogramming. Both strategies rely on codon reassignment. The expansion approach targets stop or quadruplet codons, whereas reprogramming involves reassigning sense codons to ncAAs [59]. The most popular approach utilizes the evolved tRNA and aminoacyl-tRNA synthetase (aa-RS) pair system, which consists of a suppressor tRNA, normally targeting the low-usage amber (TAG) triplet, and an orthogonal aa-RS recognizing the tRNA bearing a specific ncAA. Thus, recombinant expression of aa-RS/tRNA pairs and the protein construct of interest, together with the addition of the ncAA to the culture medium, constitute the basis for phosphomimetic protein biosynthesis [60], [61], [62]. The aa-RS/tRNA pair system has successfully been used with a wide variety of cells, both prokaryotic and eukaryotic, including yeast, plants and animals [58]. However, the effectiveness of ncAA incorporation into proteins is variable. It depends on multiple factors, such as expression levels and performance of a given orthogonal aa-RS/tRNA pair, sequence context of repurposed codons within a gene, number and/or diversity of ncAAs present in the same protein, and cellular bioavailability of ncAAs [63]. Nehring et al. [64] examined the incorporation efficiency of several Tyr analogs by cognate aa-RS/tRNA pairs, at a specific position of a target protein in yeast. Using western blot analysis followed by mass spectrometry characterization, the authors reported significant differences in the performance of the aa-RSs tested [64]. Pott et al. [65] studied the influence of nucleotide context in amber codon suppression efficiency. They generated a library of sequences around the TAG triplet and performed expression assays, identifying hit sequence contexts with relative protein yields of 70–110% (compared to the levels of control proteins without ncAAs), three to five fold higher than the average yield obtained for random sequences [65]. Schwark et al. [63] assessed the role of the release factor-1 (RF-1) in the efficiency of amber codon reassignment. The results of an expression screening with ncAA-containing green fluorescent protein (GFP) mutants, both in wild type and RF-1 deleted E. coli cells, suggested that the contribution of RF-1 in translation efficiency is highly variable in single amber codon constructs, but has a consistent and significant negative effect on multisite mutants [63].

These few examples illustrate the convenience, even necessity, of optimizing the orthogonal translation systems. This can be achieved with high-yield protein expression vectors such as pEVOL, or selectivity-enhanced aa-RS, tRNA and/or bacterial elongation factor Tu mutants for stronger interactions [58], [66], [67], [68], [69]. It is also possible to use a knock-out of the tRNA competitor RF-1 for improved binding to amber/quadruplet codons, or engineered ribosomes for higher decoding efficiency [70], [71], [72]. Moreover, the ncAA bioavailability can also be optimized through protection against phosphatases or gene deletion of these enzymes for elevated intracellular levels of phosphoresidues [73], [74].

A common issue encountered when introducing phosphorylated residues, is the lack of permeability of cell membranes to highly-charged ncAAs [62], [75]. Two main strategies have been developed to increase cellular uptake of phosphoresidues and ncAAs: 1) site-directed mutagenesis of periplasmic transporters for ncAA-specific internalization [76], and 2) addition of ncAAs to the culture medium as dipeptides with any natural amino acid, which are imported by the ATP-dependent DppA transporter and eventually hydrolyzed, releasing the ncAAs inside the cell [77]. Alternatively, the use of cell-free protein synthesis (CFPS) for ncAA incorporation has gained popularity in recent years. This approach seems to be particularly suitable for solving the intrinsic cellular bioavailability constraints of phosphoresidues, among other issues [78]. Protein expression with this method occurs in an open reaction environment to which crude cell extracts containing an overexpressed tRNA/aa-RS pair are added along with the ncAA, which is perpetually accessible to the translation machinery. There are several reports of increased production yield due to the implementation of CFPS protocols, as corroborated by western blot and mass spectrometry analyses [78], [79]. Oza et al. [80] developed a CFPS-based method to enhance the synthesis efficiency of site-specific Ser-phosphorylated constructs. Notably, the application of this system allowed them to increase the relative yield (i.e. percentage of wild type protein expression) of doubly phosphorylated human MEK1 kinase by ca. 80% compared to previously reported data for the same construct using an in vivo approach [81]. Wakamori et al. [82] took advantage of the CFPS platform to boost the efficiency of the simultaneous incorporation of acetyl-lysine into four specific positions of histone H4, achieving virtually 100% expression of the tetra-acetylated species. Yanagisawa et al. [83] also resorted to CFPS to improve the ratio of Lys-methylated Histone H3 synthesized to ncAA consumed by 10- to 20-fold compared to in vivo expression, which is particularly relevant when working with expensive and/or difficult to synthesize ncAAs [78].

Despite the evident strengths of CFPS, a great variability in reaction performance has also been documented, the causes of which can be multiple (such as ncAA sequence location in the protein construct) and unmanageable [84]. Some authors have developed strategies to obtain pSer, pThr, and pTyr by adding a protective “hood” against the action of phosphatases [73]. This protective “hood” is a photocleavable o-nitrophenylethyl group, which uncaps the phosphorylated residue upon irradiation with 350 nm light.

There is a great variety of ncAAs developed to mimic both the charge and the volume of phosphoresidues (Table 1). pTyr residue is the amino acid for which the greatest variety of ncAAs are available, in an attempt to mimic the volume and properties of its aromatic ring [73], [77], [85], [86], [87], [88], [89], [90], [91]. However, most authors choose the ncAA p-carboxy-methyl-L-phenylalanine (pCMF) because of the relatively low cost of its synthesis—it is even commercially available—and the efficient development of the technique employed to introduce it in a desired position (see below). Recently, further strategies based on the pEVOL technique have been designed to allow the synthesis and incorporation of pTyr and pThr residues in vivo, thus preventing the action of phosphatases [90], [92] (Table 1).

Case study

Many researchers tend to picture how proteins or biomolecules function by using static atomic-level structures while neglecting the dynamics of the molecule. Nonetheless, Molecular Dynamics (MD) simulations have been used increasingly over the past few years to study macromolecular systems due to the constant progress in computing power. Recently, Hollingsworth and Dror (2018) published an extensive review on MD simulations highlighting their potential application in molecular biology and drug discovery [93]. For example, MD simulations can be used to figure out how a biological system will respond to certain perturbations, such as amino acid mutation [94] or post-translational modification [95].

We herein provide a methodological article to offer a guide for readers who want to study phosphorylated proteins by combining computational and biochemical approaches. To illustrate this, two human proteins—for which phosphorylation in several residues has been studied—have been chosen as models: cytochrome c (Cc) and human antigen R (HuR) (Fig. 1). These two proteins have been selected as their number of phosphorylation sites is known and their relatively small size reduces the computational time required for simulations.

Fig. 1

Phosphorylation sites of human cytochrome c and HuR. Richardson representations of human Cc (A; PDB ID: 2N9I; [142], [113]) and HuR1-101 (B; PDB ID: 5SZW; [118], [114]). Heme group and axial ligands (His18 and Met80) of Cc are colored in green. Phosphorylation sites described in the literature for mammalian Cc and HuR are highlighted in red. Thr at position 58 is the most conserved amino acid from bacteria to mammals, but it is not present in human Cc. N-terminal and C-terminal of the polypeptide chains are marked with N and C, respectively. Secondary structure elements (α-helix and β -sheet) are depicted. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Cc is an essential metalloprotein for mitochondrial metabolism and homeostasis, playing a dual-purpose role in cell life and death. Under homeostasis, Cc acts as an electron carrier within the respiratory chain in the mitochondria. However, stress stimuli induce the release of Cc to the cytosol and nucleus where it interacts with several targets [96], [97], [98], [99], [100], [101], [102]. In the cytosol, Cc triggers programmed cell death (PCD) by activating the caspase cascade [103], [104], [105]. In the nucleus, Cc is involved in DNA damage response through its interaction with histone chaperones [99], [100]. All of these Cc functions are modulated by post-translational modifications [106], [107], [108]. In mammals, phosphorylation of residues 28, 47, 48, 58 and 97 has been reported [4], [109], [110], [111], [112], [113] (Fig. 1A). Tyr48 is the folding nucleus of Cc, and its phosphorylation alters the dynamics of the hemeprotein, thus regulating oxidative phosphorylation and PCD-related functions [114].

HuR is an RNA-binding protein that modulates the processing of mRNAs involved in cell survival, proliferation and apoptosis. HuR is composed of three RNA recognition motifs (RRMs) connected by flexible regions. The modular structure of HuR means the function of the RRM domains can be studied independently and/or in tandem [115], [116], [117], [118]. For this reason, computational and biochemical analyses of single RRMs can provide information of molecular changes caused by PTMs. Its function is highly regulated by the phosphorylation catalyzed by several kinases [119]. For example, Janus kinase 3 (JAK3) can phosphorylate Tyr63, Tyr68 and Tyr200 [120]. Tyr63 and Tyr68 are located at the first motif (RRM1), which comprises residues 1–101 (HuR1-101) (Fig. 1B). On the other hand, Tyr200 is located at the linker region between RRM2 and RRM3. Functional data suggests that Tyr200 phosphorylation impairs HuR binding to target transcripts by either shifting its localization or altering the conformation of RRMs [120]. However, the effect of phosphorylation at RRM1 remains unveiled. Understanding conformational or dynamic changes induced by this PTM may provide a rationale for changes in functionality.

Methodology

In silico approach to study protein phosphorylation

Prediction of phosphorylation sites

Tyr-phosphorylation is the most frequent PTM found in cell signaling, due to the highly conserved pTyr binding domains and the low basal phosphorylation levels of kinase protein targets, essential for a signaling mechanism [3], [121]. However, there may be several residues in a protein that can be phosphorylated. Hence, as a first approach, it is useful to check the information on phosphorylation sites provided by database servers, such as PhosphoSite ([122]), Phospho.ELM ([123]) or NetworKIN ([124]) (Table 2). Several machine learning methods have been used to generate prediction tools for phosphorylation sites. These prediction tools are based on artificial neural networks (NetPhos [125], NetPhosBac [126], NetPhosYeast[127]), convolutional neural networks (NetPhospan[128]), position specific scoring matrices (GPS [129], ScanSite [130]), hidden Markov models (KinasePhos 1.0 [131], Predikin [132]), support vector machines (PHOSIDA [133], KinasePhos 2.0 [134]), logistic regression (DISPHOS [13]) or structural analysis (pKaPS [135]) (Table 2).

Table 2

Databases and prediction tools to study protein phosphorylation.

	Databases
Phosphosite	https://www.phosphosite.org	[122]
Phospho.ELM	http://phospho.elm.eu.org	[123]
NetworKIN	http://www.networkin.info	[124]
	Prediction tools
DISPHOS	http://www.dabi.temple.edu/disphos/	[13]
GPS 5.0	http://gps.biocuckoo.cn	[129]
Scansite 4.0	http://scansite.mit.edu	[130]
pKaPS	http://mendel.imp.ac.at/sat/pkaPS	[135]
KinasePhos 2.0	http://KinasePhos2.mbc.nctu.edu.tw	[134]
NetPhos	https://services.healthtech.dtu.dk/service.php?NetPhos-3.1	[125]
NetPhosBac	https://services.healthtech.dtu.dk/service.php?NetPhosBac-1.0	[126]
NetPhosYeast	https://services.healthtech.dtu.dk/service.php?NetPhosYeast-1.0	[127]
NetPhospan	https://services.healthtech.dtu.dk/service.php?NetPhospan-1.0	[128]
PHOSIDA	http://phosida.de	[133]
Predikin	http://predikin.biosci.uq.edu.au	[132]

Analysis of phosphorylation-induced perturbations on protein structure and dynamics

Comparison between MD simulations on wild-type (WT) and phosphorylated protein species enables us to anticipate possible perturbations in the target structure and dynamics induced by phosphorylation at certain positions. MD simulations allows the evaluation of the motions of residues and atoms within a protein structure. This in silico information may be useful in assessing the behavior of phosphomimetic variants (Fig. 2). The chosen phosphomimetic variant should present similar dynamic properties and conformational changes as the phosphorylated protein.

Fig. 2

Proposed experimental approach for protein phosphorylation analysis. The decision-making at each step is the result of assessing in silico approaches with wet-lab techniques.

Due to the extensive development of MD simulations, there is a broad range of software packages that can be used to perform the simulation. Some examples include GROMACS, AMBER, CHARMM, NAMD and OpenMM [136], [137], [138], [139], [140]. Latest versions of MD software packages usually work “out of the box” with detailed manuals and tools to enable first-time use of the software. All MD simulations performed in this work were carried out on a small workstation equipped with a i7-6700 Intel CPU processor and NVIDIA 1050-Ti GPU. With this setup, we achieved a simulation speed of 50 ns/day for a 30,000 atoms system, which could be hugely increased with a high-end workstation and high-performance computing clusters. In our case, we performed MD simulations using the OpenMM [140] software in an OpenCL platform with AMBER 14SB force field [141] and particle mesh Ewald (PME) electrostatics with an Ewald summation cutoff of 10 Å.

The initial structures were based on the nuclear magnetic resonance (NMR) structures of reduced human Cc (PDB ID: 2N9I; [142]), reduced human Y48pCMF Cc mutant (PDB ID: 2N3Y; [114]) and HuR1-101 (PDB ID: 5SZW; [118]). All mutants and phosphorylated species structural models were built using the WT protein as a template. The force field parameters for phosphorylated tyrosine [143], unnatural amino acid pCMF [144] and heme group [145] were used to generate the topology and coordinate files required for the simulation, with the TLEAP module of AMBER [139].

For each simulation, the system was neutralized with sodium and chlorine counter-ions according to the total charge of the protein, and solvated with optimal 3-charge, 4-point rigid water model (OPC) water molecules [146]. The whole system was subjected to 2500 steps of energy minimization at 298 K. Langevin thermostat was used to control the temperature with a friction coefficient of 1 ps−1 and a step size of 0.002 ps. Each system was subjected to 500 ns MD simulation. To study the stability of the phosphomimetic variants, similar MD simulations were carried out at different temperatures (400 K and 500 K) and the system was subjected to 10 ns simulation.

There are several software analysis frameworks available such as MDTRAJ [147], CPPTRAJ [148] or GROMACS [138]. We opted to use the CPPTRAJ module of AMBER for trajectory analysis [148]. For each simulation, root mean square deviation of both atomic positions (RSMD, rms command on CPPTRAJ), atomic fluctuations (RMSF, atomicfluct command on CPPTRAJ), structure content (secstruct command on CPPTRAJ), hydrogen bonds (hbond command on CPPTRAJ), surface accessible to solvent (surf command on CPPTRAJ) and native contacts (nativecontact command on CPPTRAJ) were analyzed to evaluate the perturbation induced by the addition of a phosphate group or a phosphomimetic variant at a specific position of the protein. Statistics of simulated proteins throughout the simulations are listed in Table 3.

Table 3

Statistics of simulated proteins along MD trajectories.

	Average RMSD (Å)a	Average RG (Å)a	Average SAS (Å2)b	Average RMSF (Å)c
Cc WT	0.85 ± 0.10	12.96 ± 0.04	39.20 ± 10.66	0.43 ± 0.07
Cc Y48E	0.73 ± 0.07	12.87 ± 0.04	61.33 ± 13.42	0.45 ± 0.07
Cc Y48pCMF	1.93 ± 0.08	13.24 ± 0.05	128.10 ± 18.21	0.65 ± 0.14
Cc pY48	2.13 ± 0.11	13.20 ± 0.05	123.61 ± 22.67	0.66 ± 0.12
HuR1-101 WT	1.59 ± 0.18	12.14 ± 0.07	−17.25 ± 3.34	0.47 ± 0.05
HuR1-101 Y68E	1.73 ± 0.19	12.21 ± 0.07	34.96 ± 12.93	0.83 ± 0.13
HuR1-101 Y68pCMF	2.23 ± 0.25	12.33 ± 0.07	−4.71 ± 8.46	1.26 ± 0.43
HuR1-101 pY68	2.35 ± 0.29	12.40 ± 0.08	102.76 ± 28.16	1.57 ± 0.53
HuR1-101 WT	1.59 ± 0.18	12.14 ± 0.07	18.84 ± 10.54	0.87 ± 0.52
HuR1-101 Y63E	1.56 ± 0.10	12.20 ± 0.07	66.36 ± 10.54	1.10 ± 0.70
HuR1-101 Y63pCMF	1.66 ± 0.12	12.18 ± 0.06	140.11 ± 10.29	0.67 ± 0.32
HuR1-101 pY63	1.92 ± 0.10	12.23 ± 0.06	102.27 ± 21.22	0.65 ± 0.30

Average RMSD/RG values were calculated for the plateau (last 100 ns). For HuR1-101 trajectories, values were calculated for the structured domain (residues 22–97). For all simulations, RMSD drift was smaller than 0.001 Å/ns.

Average SAS values were calculated for the plateau (last 100 ns). Values were calculated for the heme group (Cc trajectories), position 68 of HuR1-101 (Tyr68 mutants’ trajectories) and position 63 of HuR1-101 (Tyr63 mutants’ trajectories).

Average RMSF values of Cc were calculated for the loop 80–91, for the stretch 42–49 of HuR1-101 Tyr68 mutants, and for the stretch 52–69 of HuR1-101 Tyr63 mutants.

Overexpression of phosphomimetic variants

Site-directed mutagenesis

Incorporation of the unnatural amino acid pCMF was accomplished by co-expressing the unnatural tRNA/aminoacyl‐tRNA synthetase pair specific for pCMF (encoded by pEVOL-pCMF plasmid) along with the protein of interest, the DNA sequence of which was mutated to a TAG amber codon at the desired incorporation site. The DNA encoding HuR1-101 was cloned into the pGEX expression vector with an N-terminal 6xHis-tag and was used as a template for the mutagenic polymerase chain reaction (PCR). AMBER mutation was performed on a pBTR1 plasmid comprising the CYCS gene encoding human Cc, along with the CYC3 gene of the yeast Cc heme lyase. CYC3 is required for proper cytoplasmic maturation of human Cc, as it covalently attaches the heme group to the Cc apoprotein. Primers used for mutagenic PCR are detailed in Table S1.

Protein expression

The expression of WT Cc and its phosphomimetic variants (Y48E and Y48pCMF) were performed as previously described [144], [149] in E. coli BL21 (DE3) strain. Briefly, ncAA-containing protein was expressed as follows: pre-culture was incubated at 37 °C at 150 rpm in minimal medium M9 supplemented with ampicillin (100 μg /mL) and chloramphenicol (20 μg/mL). Pre-induction of the pre-culture with 0.02% arabinose is necessary to ensure the presence of the unnatural tRNA/aaRS pair. Protein expression was induced at an optical density (OD600) of 0.6 with 0.02% arabinose and 1 mM isopropyl-β-d-thiogalactopyranoside (IPTG); pCMF and δ-aminolevulinic acid were added after induction at a final concentration of 1 mM each. Cells were incubated for 20 h under agitation (150 rpm) at 30 °C before harvesting by centrifugation.

WT HuR1-101 and its phosphomimetic variants were expressed in the E. coli BL21 (DE3) strain. For WT HuR1-101 and non-containing pCMF mutants, the cells were cultured at 37 °C under agitation (150 rpm) in LB medium supplemented with 100 μg/mL ampicillin. When OD600 reached 0.6, protein expression was induced with 1 mM IPTG for 20 h under agitation (150 rpm) at 30 °C. For the HuR1-101 mutants Y63pCMF and Y68pCMF, cells were cultured at 37 °C at 150 rpm in minimal medium M9, supplemented with 100 μg /mL ampicillin and 20 μg/mL chloramphenicol. Protein expression was induced at OD600 of 0.6 with 0.02% arabinose and 1 mM IPTG, as described previously [144]. pCMF was added after induction at a final concentration of 1 mM. Cells were harvested by centrifugation, resuspended in the lysis buffer 20 mM Tris-HCl (pH 8.0), containing 800 mM NaCl, 10 mM imidazole, 1 mM phenylmethanesulfonyl fluoride, 0.2 mg/mL lysozyme and 0.02 mg/mL DNase, and ruptured by sonication. Afterwards, cellular debris was separated by centrifugation at 14,000 × g for 30 min at 4 °C.

Protein expression was monitored by SDS-PAGE and Western-Blot. Proteins were separated by SDS-PAGE and transferred onto polyvinylidene difluoride (PVDF) membranes. Membranes were blocked with 5% fat-free dry milk in phosphate buffered saline (PBS)-0.1% Tween-20. Immunoblotting was performed with specific mouse primary antibody against HuR (catalog number 32160702, Sigma) or 6xHis-tag (catalog number 11922416001, Roche) and horseradish peroxidase (HRP)-conjugated secondary antibody against mouse IgG (catalog number A4416, Sigma). The immunoreactive bands were detected using Immobilon Western Chemiluminescent HRP Substrate (Millipore).

Results and discussion

General workflow

A general workflow for protein phosphorylation studies is proposed in Fig. 2. The first decision to be made before carrying out any functional and/or structural analysis is choosing the source of phosphorylated protein. If the phosphorylated protein is obtained with a high yield, either by direct purification from tissues/cells extracts or further purification after a kinase in vitro reaction, the desired physicochemical analysis can be performed.

However, as pointed out in Section 1.1, obtaining large amounts of phosphorylated protein is not always feasible. For this reason, it is necessary to overexpress a phosphomimetic variant. At this stage, an in silico approach can help to indicate which variant is the most suitable for your system.

To perform the MD simulation screening, the 3D atomic coordinates for the desired protein (preferably from the crystallographic or NMR structure) must be available at the PDB. Homology modelling software such as SWISS-MODEL (https://swissmodel.expasy.org; [150]), MODELLER (https://salilab.org/modeller; [151]) or ROBETTA (http://robetta.bakerlab.org; [152]) can be used to generate the protein coordinates if they are not available. If X-ray diffraction data from homologous proteins is available crystallographic water molecules can be included in the model. Once structural models have been generated, a comparison of phosphorylation-induced perturbations should be performed by analyzing the MD simulations of the WT protein versus the phosphorylated species (see Section 2.1.2). Subsequent analysis of multiple MD simulations of candidate phosphomimetic variants a helps to identify which molecule best mirrors the conformational changes induced by phosphorylation. Computational-based interaction network analysis may also be helpful to evaluate phosphomimetic constructs [153].

As soon as the selection of the phosphomimetic variant has been completed based on the MD data (decision-making tool), overexpression and purification of the chosen protein may be performed. This stage requires optimization of the expression conditions as mutations may affect protein folding and stability. Mimicking phosphorylation of residues exposed to the solvent – which are the most common phosphorylation sites – is often a straight forward process. However, if the phosphoresidue is located in the protein core (i.e. Tyr68 of HuR), the addition of negative charges may induce conformational changes or impair correct protein folding. MD simulations at varying conditions (temperature, ionic force or pH) can inform and detect possible difficulties in protein expression and folding.

Finally, MD simulations, in a feedback-like loop, can validate and assess the gathered functional/structural information of the phosphorylated/phosphomimetic protein.

Cytochrome c phosphorylation

A clear example where it was possible to verify the poor choice of glutamate to mimic pTyr is in Cc. The specific Cc kinase is unknown and a broad study on the functional implications of phosphorylation is unachievable with purified phosphorylated Cc because of low yields obtained from cell extracts [106], [113], [144], [154], [155]. The proposed methodology (see Section 3.1) was used to analyze different phosphomimetic variants (mutation by glutamate or substitution with pCMF) for the study of Cc Tyr48 phosphorylation (pY48) [144].

Initial MD calculations of WT Cc and modelled mutants Y48E, Y48pCMF and pY48 Cc revealed that the Y48pCMF mutation induces a localized distortion of the structure similar to that computed for the phosphorylated species. Radius of gyration (RG) values of phosphomimetic Y48pCMF and phosphorylated Cc displayed a similar behavior, while Y48E profile resembles that of WT Cc (Fig. S2). Indeed, surface accessible to solvent (SAS) analysis revealed that the solvation of the heme group in the ncAA mutant (128.1 ± 18.2 Å2) mirrors more closely the solvation experienced by the phosphorylated variant (123.6 ± 22.7 Å2) as compared with the Y48E mutant (61.3 ± 13.4 Å2) (Fig. 3 and Fig. S2). The loop containing residues 20–30 showed the highest root mean square fluctuations (RMSF) for all mutants, and were similar to values previously reported for Thr28 and Thr58 mutants [113], [155], [156] (Fig. S2). RMSF values of residues surroundings position 48 from Y48pCMF Cc were slightly higher than those from pY48 Cc, and RMSF values from Y48E Cc were slightly lower in comparison with those from WT Cc. Notably, the RMSF values at the loop between residues 80–91, where one of the axial ligands (Met80) of the heme group is located, were similar for both phosphomimetic Y48pCMF and phosphorylated pY48 Cc species, whereas WT and Y48E mutant displayed analogous values. In summary, the Y48pCMF variant showed a behavior closer to the one displayed by phosphorylated cytochrome c with respect to the overall properties calculated by MD simulations.

Fig. 3

Cytochrome c conformational changes induced by Tyr48 phosphorylation. Surface representation of human Cc (A, WT in white; B, phosphorylated at Tyr48 in orange; C, Tyr-to-Glu mutant in gold; D, Tyr-to-pCMF mutant in pale rose) after MD simulation. Heme group is colored in green and position 48 in purple. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

This data was used as a starting point to characterize the Y48pCMF Cc mutant [114]. This mutant impairs Cc diffusion between respiratory complexes, boosts hemeprotein peroxidase and reactive oxygen species (ROS) scavenging activities, and hampers caspase-dependent apoptosis [114]. Further experimental NMR relaxation measurements and hydrogen–deuterium exchange experiments on the Y48pCMF variant revealed that the structure of the ncAA mutant exhibits enhanced dynamics around the pCMF-containing loop (also observed in MD simulations), which could explain the changes observed in its functionality (Table 4). Upon cardiolipin binding, Cc undergoes conformational changes which allows H2O2 to access the heme crevice [157]. In the absence of cardiolipin, Y48pCMF Cc displayed higher peroxidase activity compared to the Y48E phosphomimetic variant. This disparity in activity may be explained by the difference in the heme group SAS, which is two-times higher for the ncAA variant. Another explanation could be the increased dynamics of the Met80-containing loop for the Y48pCMF variant. The Y48H Cc variant (naturally occurring mutation related to thrombocytopenia) also showed a higher peroxidase activity in the absence of cardiolipin with heightened dynamics in the Met80-containing loop [158].

Table 4

Comparison of biophysical and functional properties of two phosphomimetic cytochrome c mutants with respect to the native phosphorylated protein.

	in vivo pY48 Cc[110], [96]	Y48E Cc[154], [126]	Y48pCMF Cc[144], [116]
Detection by anti-pTyr	Yes	No[154], [126]	No[144], [116]
E’0, pH 7.0 (mV)	–	192.0 ± 5.0[154], [126]187.3 ± 0.8[162], [130]	208.5 ± 0.3[144], [116]
Alkaline transition (pKa)		7.0[162], [130]	6.3[144], [116]
Oxygen consumption, Vmax (% of diminished capacity)	50%	30%[154], [126]	60%[114], [100]
CcO activity, Km values (μM)	3.0	3.7[154], [126]	–
Caspase activity	–	Unable to activate caspase-3 [154], [126]Decrease by 60% (ca.) compared to WT[162], [130]	Decrease by 47.7% compared to WT[114], [100]
Bound to cardiolipin	–	30% lower affinity[154], [126]	27% more affinity[114], [100]
Peroxidase activity	–	Cc peroxidase activity only inducible at high cardiolipin concentration[154], [126]Larger than WT in the absence of CL[162], [130]	In the absence of CL, the mutant shows 3 times more activity than WT [114], [100]In the presence of CL, the peroxidase activity increases twice as much as in the absence of CL[114], [100]

– : Not available.

HuR phosphorylation

We also applied the workflow detailed above to study HuR1-101 phosphorylation. Perturbation analysis of the MD simulation of phosphorylated Tyr68 (pY68) HuR1-101 revealed that phosphorylation at this position induces conformational changes at the α1 helix (residues 33–47), thereby increasing Tyr68 solvent accessibility (Fig. 4A, Fig. S3-4). SAS of position 68 in phosphorylated HuR1-101 (102,7 ± 28.2 Å2) was higher than WT HuR1-101 (−17.3 ± 3.3 Å2). Analogous MD simulations were performed to evaluate the behavior of phosphomimetic variants (Fig. S3-4). Substitution of Tyr68 by pCMF slightly destabilizes the α1 helix, as does Tyr68 phosphorylation (Fig. 4B, Fig. S3). RMSF values for this stretch are comparable between pY68 and Y68pCMF HuR1-101 (Fig. S4). However, the solvent accessibility does not increase (−4.7 ± 8.5 Å2), displaying similar values as the WT protein (Fig. S4). On the other hand, Tyr68-by-Glu68 substitution does increase the solvent accessibility of the Glu residue (35.0 ± 13.0 Å2) without affecting the secondary structure of HuR1-101 nor inducing significant conformational changes (Fig. S4). Interestingly, the Glu side-chain flips with respect to its initial conformation, popping out of the hydrophobic pocket (Fig. 4C). MD simulations of the Y68F mutant, in which Tyr68 is replaced with Phe, was carried out as a control to check that this substitution does not alter HuR1-101 structure (Fig. 4D, Fig. S3).

Fig. 4

HuR conformational changes induced by Tyr68 phosphorylation. A, D) Superimposition of ribbon representation of Tyr68 phosphorylated (red) or mutated by Phe (purple) and non-phosphorylated (light gray) HuR1-101 structures after 500 ns of MD simulation. B, C) Overlay of the structures mutated in Tyr68 by unnatural amino acid pCMF (blue) or by canonical Glu (green) amino acids with respect to the phosphorylated Tyr68 (light gray) HuR1-101 after 500 ns of MD simulation. E) Western blots using anti-HuR (upper panel) and anti-6xHis-tag (lower panel) antibodies from E. coli cell extracts previously transformed with plasmids encoding different HuR1-101 mutants at positions 63 and 68. Molecular mass markers are also depicted (lane M). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Based on the findings of this data, we decided to study Tyr68 phosphorylation using ncAA substitution, as it better mirrors the generated perturbation. Nonetheless, HuR1-101 Y68pCMF, like Y68A and Y68E mutants, could not be expressed in vitro. In fact, the only mutant at position 68 that was successfully expressed was Y68F HuR1-101 recombinant variant (Fig. 4E). Evolutionary conservation analysis performed with ConSurf web server [159] revealed that the most conserved residue at position 68 is Tyr (88%) followed by Phe (12%). Intriguingly, ncAA HuR1-101 mutant at position 63 (Y63pCMF) and the variant Y63F were all well expressed (Fig. 4E). A plausible explanation is that Tyr68 is located at the hydrophobic core of HuR1-101, and as an aromatic amino acid key in protein folding may further stabilize the protein. To test this hypothesis, we performed several MD simulations of WT, Y68F and Y68A HuR1-101 variants at two different temperatures (400 and 500 K) to analyze their respective stability (Fig. S5). Trajectory analyses indicated that the most stable mutant at position 68 is Y68F—in fact, the only one can be successfully expressed in vitro—and that Y68F correlation between native contacts and radius of gyration throughout the whole simulation is similar to that of WT HuR1-101. Altogether, these findings suggest a crucial role for the aromatic ring (provided by either Tyr or Phe) at position 68 for correct HuR1-101 folding.

A recent NMR analysis of HuR1-101 binding to an RNA oligonucleotide derived from adenine and uracil-rich element (ARE) I of the cyclooxygenase 2 mRNA revealed that the RNA binding site comprises Tyr63 among other residues in the β 1 and β 3 strands [114]. Moreover, Tyr63 has been found phosphorylated in Hodgkin’s lymphoma cell lines [160]. Further analysis performed by Yoon et al. [120] detected Y63 phosphorylation after in vitro kinase assays with recombinant purified JAK3 kinase and HuR. As the Y63pCMF mutant is recombinantly expressed, we wanted to explore how effectively it mirrors phosphorylation at this position. For Tyr63, similar MD simulations on phosphorylated (pY63) and phosphomimetic variants (Y63E and Y63pCMF) of HuR1-101 were carried out (Figs. S6–7). Tyr63 is located at the surface of the protein, where it contacts the unstructured N-terminal tail. Perturbation analysis of the Y63E HuR1-101 mutant revealed that substitution by Glu undermines the contacts between the N-terminal tail and the protein core. This mutant displayed higher values of RG, which is consistent with partial unfolding of the N-terminal (Fig. S6). This perturbation is not induced by phosphorylation nor substitution by ncAA pCMF (Fig. 5). Moreover, RMSF values of the area surrounding position 63 are higher than those for WT, pY63 and Y63pCMF HuR1-101. Analysis of the Y63E mutant could thus led to an overestimation of the effects of phosphorylation. For this reason, the functional implications of Tyr63 phosphorylation will be assessed using the Y63pCMF HuR1-101 mutant in future studies.

Fig. 5

HuR Tyr63 phosphorylation is mimicked better by pCMF than Glu substitution. A) Superimposition of ribbon representation of Tyr63 phosphorylated (red) and non-phosphorylated (light gray) HuR1-101 structures after 500 ns of MD simulation. B-C) Overlay of the structures mutated in Tyr63 by the unnatural amino acid pCMF (blue) or by canonical amino acid Glu (green) with respect to the phosphorylated Tyr68 (light gray) HuR1-101 after 500 ns of MD simulation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Concluding remarks

The structural and functional analysis of phosphorylated proteins can be achieved by several approaches. The purpose of the herein presented workflow is to help researchers choose the most suitable alternative when an extensive study of the functional effects of phosphorylation is unfeasible with phosphorylated species purified from cells and/or tissues. MD simulations have become more accessible due to the development of user-friendly software packages and the increase in computing power at a modest cost. However, computational approaches are often ignored by molecular biologists. For this reason, we believe that a combination of computational and biochemical approaches will help to better understand how phosphorylation events affect protein structure and function.

CRediT authorship contribution statement

Gonzalo Pérez-Mejías: Conceptualization, Methodology, Software, Visualization, Investigation, Writing - original draft, Writing - review & editing. Alejandro Velázquez-Cruz: Writing - original draft, Writing - review & editing. Alejandra Guerra-Castellano: Investigation, Writing - original draft, Writing - review & editing. Blanca Baños-Jaime: Resources. Antonio Díaz-Quintana: Software, Writing - review & editing. Katiuska González-Arzola: Writing - review & editing. Miguel Ángel De Rosa: Supervision, Writing - review & editing, Project administration, Funding acquisition. Irene Díaz-Moreno: Conceptualization, Supervision, Writing - review & editing, Project administration, Funding acquisition.