KRAS K104 modification affects the KRASG12D-GEF interaction and mediates cell growth and motility

Mutant RAS genes play an important role in regulating tumors through lysine residue 104 to impair GEF-induced nucleotide exchange, but the regulatory role of KRAS K104 modification on the KRASG12D mutant remains unclear. Therefore, we simulated the acetylation site on the KRASG12D three-dimensional protein structure, including KRASG12D, KRASG12D/K104A and KRASG12D/K104Q, and determined their trajectories and binding free energy with GEF. KRASG12D/K104Q induced structural changes in the α2- and α3-helices, promoted KRAS instability and hampered GEF binding (ΔΔG = 6.14 kJ/mol). We found decreased binding to the Raf1 RBD by KRASG12D/K104Q and reduced cell growth, invasion and migration. Based on whole-genome cDNA microarray analysis, KRASG12D/K104Q decreased expression of NPIPA2, DUSP1 and IL6 in lung and ovarian cancer cells. This study reports computational and experimental analyses of Lys104 of KRASG12D and GEF, and the findings provide a target for exploration for future treatment.

Calculation of the free energy of binding. The hybrid topologies of the systems for the free energy calculations were constructed using the pmx tool 22,23 . This tool automatically generates hybrid structures and topologies for amino acid mutations that represent the two physical states of the system. The double system (combining both branches of the thermodynamic cycle) in a single simulation box was set to keep the system neutral at all times during a transition 22 . After obtaining this hybrid structure, we used GROMACS to perform free energy simulations. One hundred snapshots were extracted from the 10-ns equilibrated trajectories, and a rapid 100-ps simulation was performed starting from each frame. Lambda (λ) ranged from 0 to 1 and from 1 to 0 for the forward and backward integrations, respectively, thus describing the interconversion of the WT to the mutant system and of the mutant to the WT system, respectively. Finally, the pmx tool was used to subsequently estimate free energy differences and integrate the multiple curves by the fast-growth thermodynamic integration approach, which relies on the Crooks-Gaussian Intersection (CGI) protocol 24 . In such a double-system/singlebox setup, the estimated free energy corresponds to the ΔΔG of binding.
Western blotting and immunoprecipitation. Protein was extracted from cells in lysis buffer that contained a protease inhibitor cocktail (Roche Applied Science Reverse transcription (RT)-PCR. RNA was isolated from cell lines using TRIZOL reagent (Invitrogen).
cDNA was synthesized using a cDNA Synthesis kit (Promega, Madison, WI, USA). Quantitative PCR was performed with 10 ng of cDNA using Power SYBR Green PCR Master Mix (Applied Biosystems), and the cRNA was analyzed with an RT-PCR System kit (StepOne, Applied Biosystems, Darmstadt, Germany); 18S was used as the reference gene for normalization. All values were calculated using the 2 −ΔΔCT method.

Results
Modeling of the KRAS and KRAS/GEF complex. The KRAS model was constructed using the published human KRAS crystal structure (PDB ID: 3GFT) as the template (Fig. 1A). Of the two replaced residues, Gly12 is located at the GDP-binding site, and Lys104 is located at the C-terminal end of the α3 helix, which is near the GEF-binding domain. Figure 1B shows the KRAS/GEF complex that was constructed by using the crystal structure of a ternary Ras:SOS (PDB ID: 1XD2) from Homo sapiens as the template. GEF activates monomeric KRAS by stimulating the release of GDP to allow binding of GTP.

Structural changes in the WT and mutant KRAS/GEF complex in conventional MD simulations.
The trajectories from MD simulations of the WT and mutant KRAS/GEF complex in explicit solvent models were calculated. The RMSD values of the backbone atoms for the WT and mutant KRAS/GEF complex with regard to the initial structures were plotted ( Fig. 2A), and the quality and convergence of the simulated dynamic trajectories can be estimated by assessing the RMSD values. With regard to the simulations of the WT and mutant KRAS/GEF complex, the data obtained demonstrated that after a rapid increase during the first 0.2 ns, the trajectories were stable, with average values of 4.36, 4.14, 3.84 and 4.44 Å for the WT, G12D, G12D + K104A and G12D + K104Q KRAS/GEF complexes, respectively. In contrast, the RMSD values showed a large increase as a function of time for the G12D + K104Q KRAS/GEF complex ( Fig. 2A). The results reveal that the conformational changes of the G12D + K104Q KRAS/GEF complex were significant when compared with other complexes.
To understand the structural changes of KRAS in the GEF-binding region caused by these mutations, we detected the conformational difference between the WT and mutant KRAS/GEF complex structures by measuring the dihedral angle of the α2-helix and α3-helix (Cα in Met67, Thr74, Lys104 and Thr87) (Fig. 2B). Figure 2C shows the dihedral angle of the α2-helix and α3-helix as a function of time. In Fig. 2D, histograms of the dihedral angle during the last 30 ns of the simulations are depicted; the average dihedral angles are − 55.3°, − 55.8°, − 50.3° and − 62.0°. As illustrated in Fig. 2D, the most populated dihedral angle of the α2-helix and α3-helix in the WT www.nature.com/scientificreports/ (blue) and G12D (red) KRAS/GEF complexes shows a similar distribution of the degree. The most populated dihedral angle of the G12D + K104A KRAS/GEF complex (green) is ~ 5° larger and more widely distributed than that in WT KRAS (blue). That is, the mutation causes the region of α2 and α3 helices of the G12D + K104Q KRAS/ GEF complex to be more flexible. Interestingly, a different phenomenon was observed for the G12D + K104Q KRAS/GEF complex (purple), where a left-shifted distribution of the dihedral angle of the α2-and α3-helices was detected by comparison with the WT KRAS/GEF complex (blue), with the middle value being ~ 7° smaller in the G12D + K104Q KRAS/GEF complex (purple) (Fig. 2D). The results suggest that more structural changes of the α2-and α3-helices were prevalent in the G12D + K104Q KRAS/GEF complex, which is consistent with the RMSD results ( Fig. 2A). The K104Q mutation induced an additional structural change in the α2-and α3-helices, which may be involved in the support and stabilization of GEF binding. Based on these analyses, the KRAS K104Q mutation is predicted to affect GEF binding.

Residue fluctuation of the WT and mutant KRAS/GEF complex structures.
The results of B-factor calculations for each residue revealed that the atomic fluctuations of the G12D/K104Q mutant were significant at the KRAS/GEF interaction regions (especially in the α25-helix of the GEF protein) when compared with WT and with the G12D and G12D/K104A mutants (Fig. 3). This is consistent with the left-shift distributed dihedral angle of the G12D/K104Q mutant relative to WT and the other two mutants. The data from MD simulations of the KRAS/GEF complex suggest that the K104Q mutation leads to a considerable decrease in the affinity of KRAS for GEF.

Analysis of KRAS/GEF interactions.
KRAS/GEF binding free energy differences between the G12D, G12D + K104A and G12D + K104Q mutants were calculated in a series of alchemical free energy simulations using the CGI protocol 24 . The GEF binding affinity differences ( G Mutation Binding ) were calculated according to ΔG 1 (GEF-bound) -ΔG 2 (GEF-free), as shown in the thermodynamic cycle that was distributed between two simulation boxes for λ = 0 (green) and λ = 1 (blue; Fig. 4A). Thus, the change in the net charge during the alchemical simulation remains zero. Notably, stabilizing mutations have negative ΔΔG values. We adopted the CGI protocol www.nature.com/scientificreports/ to calculate the binding affinity differences between KRAS and GEF and the single point mutations. The resulting binding free energy differences ( G Mutation Binding ) between G12D KRAS and G12D + K104A KRAS and between G12D KRAS and G12D + K104Q KRAS were − 10.80 and 6.14 kJ/mol, respectively (Fig. 4A). The results reveal that the KRAS-GEF interaction may be disrupted in the G12D + K104Q mutant (i.e., ��G Mutation Binding > 0 ), and thus, GEF may not be able to activate this KRAS mutant by stimulating the release of GDP to allow binding of GTP.
In addition to the alchemical free energy calculation, we examined the level of protein-protein interaction between GEFs, effectors and KRAS by experimental methods. We used immunoprecipitation-western blotting and pull-down assays and detection of active Ras when KRAS WT , KRAS G12D , KRAS G12D/K104A and KRAS G12D/K104Q were individually expressed in cell lines. KRAS G12D/K104A and KRAS G12D/K104Q plasmids were generated with pBabe-Kras G12D (Addgene) using site-directed mutagenesis. Active Ras binds to the Ras-binding domain (RBD) of Raf1, leading to Ras activation. The results of the analysis showed increased interaction of KRAS G12D/K104Q with Ras GTPase-activating protein 1 (RASA1, RasGAP) but decreased interaction with GEFs (Fig. 4B) as well as effectors, reducing the binding to RBD of Raf1 and leading to Ras inactivation (Fig. 4C,D) compared with KRAS G12D and KRAS G12D/K104A . We also analyzed the level of RAS in the total MCAS and H1299 cell lysates, RAS downstream pathways (phospho-AKT) and acetylation. According to the results, the levels of RAS of KRAS WT , KRAS G12D , KRAS G12D/K104A and KRAS G12D/K104Q were the same. However, ac-lysine of RAS was not detected. KRAS G12D and KRAS G12D/K104A exhibited increased levels of p-AKT compared to WT and KRAS G12D/K104Q (Fig. 4E,F).

K104 modification affects KRAS G12D -mediated cell growth and motility. To investigate whether
K104 modification affects KRAS G12D , we used KRAS WT , KRAS G12D , KRAS G12D/K104A and KRAS G12D/K104Q and analyzed the biological function of these mutations in lung and ovarian cell lines. Individual expression of KRAS G12D , KRAS G12D/K104A and KRAS G12D/K104Q was induced by transfection of plasmid DNA in H1299 and MCAS cancer cell lines. We then analyzed their effects on cell growth (Fig. 5A,B), wound healing (Fig. 5C,D) and invasion (Fig. 5E,F). KRAS G12D and KRAS G12D/K104A promoted cell growth, wound healing and invasion, whereas KRAS G12D/K104Q attenuated these processes. Moreover, compared to KRAS G12D , KRAS G12D/K104A significantly increased the rate of cell growth (P = 0.019) and migration ability in H1299 (P < 0.0001) and MCAS (P = 0.02) cell lines. These results suggest that K104 modification regulates KRAS G12D -mediated cell growth and motility in lung and ovarian cancer cell lines.

KRAS G12D/K104Q mediates NPIPA2, DUSP1 and IL6 expression.
Our previous study found that KRAS G12D/K104Q affects cellular biological functions by associating with GTP binding. However, the related downstream genes were not clear. Therefore, we sought to identify downstream genes of KRAS G12D/K104Q that are reduced compared with KRAS G12D/K104A using a whole-genome cDNA microarray to screen KRAS G12D (> twofold), KRAS G12D/K104A upregulation (> twofold) and KRAS G12D/K104Q downregulation (> twofold) of gene expression in the lung cancer cell line H1299 and the ovarian cancer cell line MCAS (Fig. 6A). Further selection of the top three and four assessments of the microarray results showed that NPIPA2, DUSP1 and IL6 were consistently identified as upregulated genes in KRAS G12D -expressing cells and KRAS G12D/K104A -expressing cells and downregulated in KRAS G12D/K104Q -expressing cells relative to control cells (Fig. 6B). To confirm the microarray results, NPIPA2, DUSP1 and IL6 expression was examined by qRT-PCR, and expression was induced by KRAS G12D and KRAS G12D/K104A and decreased by KRAS G12D/K104Q in cancer cell lines (Fig. 6C,D). These results were obtained in both cell lines, with the three genes showing the same pattern in both. These results suggest that KRAS G12D/K104Q decreases expression of NPIPA2, DUSP1 and IL6 in lung and ovarian cancer cells.

Discussion
Based on the crystal structure of KRAS bound to GDP (PDB code 4LPK) 33 , K104 is located on the C-terminal end of the α3-helix, and the positively charged side chain directly interacts with a carbonyl group at the negative pole of the α2-helix. These interactions have been predicted to play a key role in the structural stability of the KRAS-GEF complex 15 . Previous studies have indicated that K104 modification may affect both guanine nucleotide exchange and GTP hydrolysis rates and influence the downstream signal transduction pathways involved in the control of cell survival. Marcus et al. showed that HRAS K104Q decreases intrinsic and catalyzed hydrolytic behavior but that K104A does not affect intrinsic hydrolysis 34 . Based on the NIH3T3 proliferation assay, Yang et al. observed that the G12V and G12V/K104A forms of KRAS enhanced cell survival but that G12V/ K104Q did not 15   www.nature.com/scientificreports/ exhibits defects in both GEF-mediated exchange and GAP-mediated GTP hydrolysis and also indicated that KRAS K104Q did not alter steady-state GTP-bound levels or the ability of the oncogenic KRAS G12V mutant in NIH3T3 mouse fibroblasts 35 . From these findings, it is clear that these investigations focused on the K104 modification of RAS-WT or RAS-G12V; however, very few attempts have been made to target RAS-G12D. In this study, we focused on the KRAS-G12D mutant and investigated how KRAS K104 modification affects the KRAS G12D -GEF complex interaction and mediates cell growth and motility. Based on our results, K104Q is able to reduce GTP-bound levels of KRAS-G12D and the oncogenic effects of the KRAS-G12D mutant in H1299 and MCAS cells, but K104A cannot. www.nature.com/scientificreports/ MD simulations of the KRAS/GEF protein complex structure revealed that the conformational changes of the G12D + K104Q mutant were significant at the α2-and α3-helices when compared with the WT, G12D and G12D + K104A KRAS/GEF complexes (Fig. 2). Moreover, the K104Q mutation may induce additional fluctuations at the KRAS/GEF interaction regions (especially in the α25-helix of the GEF protein) (Fig. 3). As mentioned above, the α2-and α3-helices play important roles in the binding of GEFs; therefore, we postulate that such fluctuations may promote instability in this region, which consequently reduces the binding ability between KRAS and GEFs. The data obtained from the binding free energy differences indicate that the binding of GEF with the K104Q mutant is less favorable than that of GEF with the K104A mutant (Fig. 4A). These results are consistent with previous studies 15, 28,35 showing that KRAS K140Q decreases the SOS-catalyzed nucleotide exchange rate.
In this study, we found that the effect of the mutant G12D/K104A is higher than that of G12D. This phenomenon was confirmed in both our computational and experimental analyses. Our computational studies showed that the G12D and G12D/K104A forms of the KRAS-GEF complex have different effects on the interactive regions (α2-and α3-helices). According to RMSD results, the most populated dihedral angle of the G12D + K104A KRAS/GEF complex (green) is ~ 5° larger and more widely distributed than that in G12D KRAS (red) (Fig. 2D). Additionally, the data obtained from the binding free energy differences indicate that the binding of GEF with the G12D/K104A mutant is more favorable than that of GEF with the G12D mutant (ΔΔG = − 10.80 kJ/mol) (Fig. 4A). This result is also consistent with our immunoprecipitation assay (Fig. 4B). Moreover, our experimental studies showed that G12D/K104A slightly increased the active form of KRAS compared with the G12D mutant (Fig. 4C,D). Therefore, we hypothesize that the G12D/K104A mutant allows GEF to have a greater ability to activate KRAS by removing GDP than the G12D mutant. As a result, the G12D/K104A mutant may have a more aggressive oncogenic phenotype than the G12D mutant. Regardless, the detailed mechanisms need to be further investigated.
We also show that KRAS G12D/K104Q induces structural changes in the α2-/α3-helix to block binding with GEF and mediate tumor formation. KRAS G12D/K104Q also induced expression of the downstream target genes NPIPA2, DUSP1 and IL6. A previous study suggested that Lys104 of KRAS is easily acetylated by posttranslational modification 15 to regulate cellular biological function. KRAS G12V K104 acetylation (as represented by K104Q) represses downstream signaling and inhibits cell growth through HDAC6 and SIRT2 in NIH 3T3 cells 36 . Previous studies as well as our own observed that KRAS G12V and KRAS G12D have the same function with respect to regulating the KRAS pathway and cellular function when K104 is acetylated. Although KRAS G12D can activate cancer cells, progression to malignancy requires additional genetic lesions 37 . Among cancers with the highest mortality rates, KRAS-activating mutations occur in ~ 90% of pancreatic cancers, 40% of colon cancers and 30% of lung cancers. Therefore, it is important to consider cancer therapies that target KRAS mutations. Research has predominantly focused on EGFR mutation for lung cancer-related drug development, and there are many EGFR antagonists that have been used for clinical treatment. Unfortunately, the mortality rate for lung cancer remains high throughout the world. Indeed, epidemiological statistics indicate that EGFR and KRAS mutations are present in ~ 50% of lung adenocarcinomas, of which ~ 20% are EGFR mutations and 26% KRAS mutations 38 . KRAS mutations occur in 3.7-36.4% of endometrioid ovarian cancers 39,40 , though therapies targeting KRAS have rarely been used to treat ovarian cancer. Hence, it is very important to focus on KRAS mutations for the development of future therapies to treat lung and ovarian cancer.
Although expression of NPIPA2, DUSP1 and IL6 was decreased after acetylation of KRAS, it is still unclear whether this is related to the ubiquitin system. This possibility will be pursued in our future studies. Finally, our findings suggest that KRAS G12D/K104A regulates GEF to activate KRAS signaling and promote cell growth, invasion and migration in lung and ovarian cancer cell lines. Mutationally activated KRAS (through G12D and G12V) is still a difficult pharmacological target. Next, we will analyze K104 acetylation of KRAS G12D using histone deacetylase inhibitors (HDACis), which regulate cellular functions, including cell growth, invasion, migration and apoptosis, by repressing expression of HDAC, leading to increased acetylation of lysine residues in target proteins 41,42 . The HDACi trichostatin A (TSA) is a naturally derived hydroxamic acid, which has been studied as a potential nontoxic anticancer drug 43,44 . It is hoped that an HDACi such as TSA-either alone or in combination with other chemotherapy approaches-may provide a viable strategy for clinical targeting options in lung and ovarian cancer.