Computational Study of Symmetric Methylation on Histone Arginine Catalyzed by Protein Arginine Methyltransferase PRMT5 through QM/MM MD and Free Energy Simulations

Protein arginine methyltransferases (PRMTs) catalyze the transfer of the methyl group from S-adenosyl-l-methionine (AdoMet) to arginine residues. There are three types of PRMTs (I, II and III) that produce different methylation products, including asymmetric dimethylarginine (ADMA), symmetric dimethylarginine (SDMA) and monomethylarginine (MMA). Since these different methylations can lead to different biological consequences, understanding the origin of product specificity of PRMTs is of considerable interest. In this article, the quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) and free energy simulations are performed to study SDMA catalyzed by the Type II PRMT5 on the basis of experimental observation that the dimethylated product is generated through a distributive fashion. The simulations have identified some important interactions and proton transfers during the catalysis. Similar to the cases involving Type I PRMTs, a conserved Glu residue (Glu435) in PRMT5 is suggested to function as general base catalyst based on the result of the simulations. Moreover, our results show that PRMT5 has an energetic preference for the first methylation on Nη1 followed by the second methylation on a different ω-guanidino nitrogen of arginine (Nη2).The first and second methyl transfers are estimated to have free energy barriers of 19–20 and 18–19 kcal/mol respectively. The computer simulations suggest a distinctive catalytic mechanism of symmetric dimethylation that seems to be different from asymmetric dimethylation.


Introduction
Posttranslational methylation of histone proteins on their arginine residues is an epigenetic mark that plays a vital role in cell function and is related with cell disorders and diseases [1]. Protein arginine methyltransferases (PRMTs) catalyze the transfer of methyl group(s) from S-adenosyl-L-methionine (AdoMet) to the guanidine nitrogen of arginine residue, resulting in the reaction products containing methylarginine along with S-adenosylhomocysteine (SAH) [2]. In addition to the methylation of histone proteins, PRMTs could also modify a variety of other proteins [3]. Depending on types of PRMTs (I, II or III), the methylation products may contain asymmetric dimethylarginine (ADMA), symmetric dimethylarginine (SDMA) or monomethylarginine (MMA) (shown in Scheme 1) [1]. Scheme 1. Methylation of arginine by different types of PRMTs. Type I PRMT can produce both monomethylarginine (MMA) and asymmetric dimethylarginine (ADMA). Type II PRMT can produce both MMA and symmetric dimethylarginine (SDMA). Type III PRMT can only produce MMA.
The Type II PRMT5 to be investigated in this work catalyzes the transfer of methyl groups to the two different ω-guanidino nitrogen atoms on the arginine residue of the target protein, producing the ω-N G , N′ G symmetrically dimethylated arginine (SDMA) [4]. PRMT5 has a variety of substrates that include histones, transcription factors, splicesomal proteins and piRNA biogenesis related proteins, and this enzyme functions in both nucleus and cytoplasm [3]. SDMA may profoundly impact many biological processes including epigenetic control of gene expression [5], circadian rhythms [6,7], splicing regulation [8,9], germ cell development and pluripotency [10,11], and DNA damage response [12,13]. The Type II PRMT5, however, often share the common recognition sequence with the Type I PRMTs which add two methyl groups to the same ω-guanidino nitrogen atom (ADMA) [14]. Thus, the same target arginine may be either symmetrically or asymmetrically dimethylated. As such, SDMA and ADMA are isomeric protein posttranslational modifications with distinct and sometimes reversible biological effects [14]. One example is the methylation of arginine 3 on histone H4 (H4R3): symmetric dimethylation of H4R3 could repress gene expression [15,16], while asymmetric dimethylation of H4R3 is correlated with gene activation [17,18]. Recent studies have focused on understanding the enzymatic mechanisms that differentiate the two chemically isomeric but functionally antagonistic posttranslational modifications [14].
The crystal structures have been determined for several PRMTs, including PRMT1 [19], PRMT3 [20], PRMT6 [21], PRMT10 [22] of Type I, PRMT5 [23] of Type II and PRMT7 [24] of Type III. Computer simulations have been applied to investigate the catalytic mechanism of PRMT1 [25] and PRMT3 [26], and some important questions concerning the product specificity of Type I ADMA have been addressed. The computational approaches used in these earlier studies include molecular dynamics and free energy simulations (potential of mean force) with the hybrid quantum mechanical and molecular mechanical (QM/MM) potential that seem to be suitable to investigate the enzyme-catalyzed methyl transfer process and have also been widely used for some other methyltransferases [26][27][28][29][30][31][32].
In this study, the methylation reaction catalyzed by Type II PRMT5 is investigated by use of QM/MM MD and free energy simulations on the basis of experimental observation that the dimethylated product is generated through a distributive fashion. The simulations have identified some important interactions and proton transfers involving the active site residues. Similar to the cases of Type I PRMT1 [25] and PRMT3 [26], a homologous Glu residue (Glu435) in PRMT5 seems to function as a general base catalyst during the catalysis, and the corresponding proton transfer is found to be somehow concerted with the methyl transfer process. However, unlike Type I PRMTs which energetically favors a single ω-guanidino nitrogen (Nη2) of arginine as the target for the both 1st and 2nd methylations [25,26], PRMT5 is found to have an energetic preference of targeting Nη1 for the first methylation and then targeting a different ω-guanidino nitrogen (Nη2) for the second methylation. The first and second methyl transfers are estimated to have free energy barriers of 19-20 and 18-19 kcal/mol, respectively, from the simulations. These results are consistent with the existing experimental data [23]. Our computational study provides a better understanding of the symmetric di-methylation mechanism that is different from the asymmetric di-methylation.

Comparison of the Active Site Structures of Type II PRMT5 and Type I PRMT3
The two invariant glutamate residues from the "double-E" loop found in all PRMTs [33] are structurally conserved in both Type I PRMT3 and Type II PRMT5 active sites (E435 and E444 in Figure 1A; E144 and E153 in Figure 1B). As is shown in Figure 1, the carboxylate side chain of E444 in PRMT5 (E153 in PRMT3) forms a stable salt bridge with the Nε and Nη1 atoms of the substrate arginine. The side chain of E435 in PRMT5 (E144 PRMT3) interacts with Nη1 and Nη2 of the substrate. These two Glu residues are required for the enzymatic activities; the mutation of either of them could greatly decrease the enzymatic activity [14]. Previous computational studies have identified E144 as the general base to accept proton from the arginine during the methyl transfer catalyzed by Type I PRMTs [25,26]. One interesting questions is whether the corresponding E435 from PRMT5 would play a similar role during the catalysis. Four other residues that are conserved in the active site of PRMT5 are F327, K333, S578 and S439; the corresponding residues in PRMT3 are M48, R54, Y148 and H293, respectively. Although the mutation of S439 and S578 of PRMT5 diminishes the enzymatic activity significantly [14], the exact role of these residues is not clear; both S439 and S578 seem to be far away from the methyl donating AdoMet ( Figure 1A). Interestingly, the F327M mutant of PRMT5 could produce both ADMA and SDMA [14], indicating that F327 may occupy a key position in PRMT5 and its properties may be important in determining the product specificity. Moreover, K333 is in the vicinity of F327 and could form hydrogen bonds with both the carboxylate group of AdoMet and E435.
The previous computational studies have shown that for PRMT1 and PRMT3 the 1st and 2nd methyl transfers would be energetically more favorable with relatively lower barriers if Nη2 (see Figure 1B) is the methyl acceptor in each of the cases [25,26]. This is consistent with the fact that they are both Type I PRMTs and the methylation products contain asymmetric dimethylarginine (ADMA). In the active site of the PRMT5 structure ( Figure 1A), Nη1 appears to be in a much better position for acting as the methyl acceptor compared to Nη2; this is in contrast with the case of PRMT3 where Nη2 seems to be in a better position for accepting the methyl group [26]. Indeed, the distance between Nη1 and CM is 3.2 Å in the active site of PRMT5 as compared to 4.9 Å between Nη2 and CM.

The First Methylation Catalyzed by PRMT5
The acceptor site for the first methylation process catalyzed by PRMT5 is studied by comparing the energetics of the methyl addition to Nη1 and Nη2, respectively, on the arginine via QM/MM MD and free energy simulations. As shown in Figure 2B, the free energy barrier for the methyl transfer to Nη1 is 20.4 kcal/mol, which is 9 kcal/mol lower than that for the methyl transfer to Nη2. This suggests that the 1st methyl group is likely to be transferred to the Nη1 atom, as the corresponding process is energetically more favorable. This is in contrast with the 1st methylation catalyzed by PRMT3 which favors Nη2 as the methyl acceptor [26]. The average active-site structures of the reactant complex and near transition state for the 1st methyl transfer are shown in Figures 3 and 4 for the methyl transfer to Nη1 and Nη2, respectively. The interactions obtained from the simulations for the reactant complex generally resemble the interactions observed in the crystal structure. This seems to indicate that the conformation of the active site in the crystal structure is well maintained after initial 1.5 ns MD simulations. In the reactant complex for the methyl transfer to Nη1 ( Figure 3A), Nη1 is well aligned with the methyl group of AdoMet, and this good alignment presumably leads to a relatively low barrier for the methyl transfer. The arginine on the substrate appears to be well stabilized through the salt bridge with E444 as well as hydrogen bonding interaction with E435. K333 may also help to adjust the orientation of E435 through a hydrogen bonding interaction. F327 seems to be involved in the π-cation interaction with the arginine in the reactant state. Moreover, as mentioned earlier the bulk size of this residue may prevent the formation of a good alignment between the transferable methyl group and Nη2 and interfere with the methyl transfer to Nη2. Indeed, the structures in Figures 1A shows that the distance between the transferable methyl group (CM) and Nη2 can be as much as 4.9 Å in the reactant complex. By contrast, the corresponding distance in PRMT3, which has M48 at the location rather than a Phe residue, is much smaller (e.g., 3.9 Å in Figure 1B). It is of interest to note from Figure 4B that the distance between F327 and the arginine increases significantly near the transition state during the methyl transfer to Nη2.  The 1D free energy simulations only used the reaction coordinate for the methyl transfer. Nevertheless, the proton transfer occurred near the transition state of methylation ( Figures 3B and 4B). To better understand the relationship between the methyl transfer and the proton transfer process, the 2D free energy simulations were performed ( Figure 5). As is shown in Figure 5, the reaction path goes from point A (designated as the reactant complex) to point B (the transition state), and finally reaches point D as the product complex. Before reaching the transition state at point B, the proton has been basically transferred from Nη1 of the substrate arginine to the carboxyl oxygen of E435. And the free energy barrier is estimated to be 19-20 kcal/mol from the 2D free energy simulations, which is almost the same as that obtained from the 1D free energy simulations (20.4 kcal/mol). Moreover, point D in Figure 5 is about 6 kcal/mol lower than point C, suggesting that the deprotonated product is more stable in the active site than the protonated product and that the proton transfer is coupled with the methyl transfer during methylation.

The Second Methylation Catalyzed by PRMT5
With the first methyl group being transferred already, the 2nd methylation reaction can now be examined to determine the mechanism for the formation of the symmetrically dimethylated product. Previous experimental study on Caenorhabditis elegans PRMT5 (cPRMT5) has shown that the dimethylated product is generated through a distributive fashion [4] in which the peptide is released prior to rebinding to facilitate a second round of methylation. The distributive mechanism was further confirmed from the study of the human hPRMT5•MEP50 complex [34]. This is in contrast with the cases of some protein lysine methyltransferases for which the multiple rounds of methylation are believed to proceed processively without the release of the intermediates from the active sites. The kinetic data (Km and kcat) for cPRMT5, hPRMT5 and hPRMT5•MEP50 complex have also been obtained for a variety of un-methylated and mono-methylated substrates. For the hPRMT5•MEP50 complex for which the current investigation is based on, the differences in Km and kcat for the first and second methyl transfers are rather small and are beyond the accuracy of most of quantum mechanical approaches used in QM/MM studies. Although several experimental structures for PRMTs are available, the structure for the reactant complex of the second methyl transfer has not been determined. In our earlier investigations of protein lysine methyltransferases [27][28][29][30], we developed an approach based on the free energy simulations of the methyl transfer processes that allows us to determine whether or not the enzymes would be able to add methyl groups. This approach has been successfully applied to study the product specificity for a number of protein lysine methyltransferases, and its usefulness has also been confirmed on some other methyltransferases [26,31,32], including PRMT1 [26]. This approach is especially suitable for the current investigation because the Km values for the first and second methyl transfers are rather similar (see above), and the hypothetical non-existent methyl transfers can be presumably eliminated using this approach in the similar manner as the product specificity is determined [26][27][28][29][30]. In Table 1, the three different configurations for the mono-methyl arginine without destroying the salt bridge with E444 are given. The nitrogen atom already connecting to the methyl group is designated as Nη1 and that without the methyl group (i.e., NH2) designated as Nη2. As is demonstrated in Table 1, for each configuration there are two possible di-methylation products (SDMA and ADMA). The QM/MM free energy simulations have been performed for each of the six methylation processes. Table 1(I) shows that the free energy barriers were calculated to be 32.3 and 36 kcal/mol for the 2nd methyl transfer to the Nη1 and Nη2 atoms, respectively. Since either of these energy barriers seems to be too high, the 2nd methylation may not start from the configuration in Table 1(I). The similar argument can be made for the configuration in Table 1(III). Table 1. Free energy barriers for the 2nd methylation on either the mono-methylated nitrogen atom (designated as Nη1) or un-methylated nitrogen atom (designated as Nη2) of the mono-methylated arginine (MMA) with three possible configurations for the reactant complex without destroying the salt bridge involving E444. The only difference between Configurations II and III is whether the methyl group is pointing to the AdoMet side (III) or away from the AdoMet side (II). The conformation in (II) has the lowest energy barrier of 20.1 kcal/mol when the 2nd methylation occurs to Nη2 which will eventually produce the symmetrically dimethylated arginine (SDMA). Error bars are also determined, similar to those given in Figure 2 (within ± 0.2 kcal/mol). shows that the second methyl transfer has the lowest free energy barrier (20.1 kcal/mol) if it is transferred to Nη2 (i.e., the SDMA product), while the free energy barrier for the second methyl transfer to Nη1 (i.e., the ADMA product) based on the configuration in Table 1(II) is 11.2 kcal/mol higher.
The average structures of the reactant complex and near transition state for the 2nd methyl transfer to Nη2 and Nη1 are exhibited in Figures 6 and 7, respectively. The active-site structures of the reactant complex ( Figure 6A) show that the Nη2 atom is well aligned with the methyl group of AdoMet. The relatively high barrier for the 2nd methyl transfer to Nη1 appears due in part to the steric hindrance of F327 (as in the case of the 1st methyl transfer), although other factors may be involved as well. One interesting observation for the second methyl transfer is that the proton has not been transferred to the general base E435 at the transition state (see Figures 6B and 7B). This is in contrast with the case of the first methyl transfer (shown in Figures 3B and 4B) where the proton transfer occurs before the methylation reaches the transition state. This observation is also reflected in the 2D free energy map involving both the 2nd methyl transfer and the proton transfer ( Figure 8). As is shown in the reaction path illustrated in Figure 8, the reaction reaches the transition state at point B without the deprotonation of the monomethylarginine.  Another interesting observation from Figures 6B and 7B is that E435 forms a salt bridge to both Nη1 and Nη2 near the transition state. This salt bridge may be relatively strong compared to that in the reactant complex and help to stabilize the transition state; it may also play a role for the delayed proton transfer mentioned earlier. The free energy barrier is estimated to be 18-19 kcal/mol by the 2D free energy simulation (Figure 8), consistent with 20.1 kcal/mol from the 1D free energy simulations (Table 1).

Experimental Section
The simulation coordinates for the reactant complexes of methyl transfers were based on the crystallographic complex (PDB ID: 4GQB, 2.06 Å) of PRMT5 [23] that contains AdoMet analog and the H4 peptide with the methylation targeting substrate, Arg3. The coordinates of PRMT3 for comparison is from the crystal structure (PDB ID:1F3L, 2.06 Å) and modified based on the previous work [26]. The QM/MM MD and free energy (potential of mean force, PMF)) simulations were applied for monitoring the methylation processes and determining the free energy profiles with the CHARMM program [35]. A water sphere based on a modified TIP3P water model [36] with radius(r) of 30 Å, centered at CZ of Arg3, was pre-equilibrated to the system. A stochastic boundary with a Poisson-Boltzmann charge-scaling scheme [37] was applied for the model. The reservoir region had r > 22 Å, and the buffer region had r equal to 20 Å ≤ r ≤ 22 Å. The reaction region had r ≤ 20 Å. The -CH2-CH2-S + (Me)-CH2-part of AdoMet, the side chain of substrate Arg3/monomethylated Arg3, the side chains of E435 and E444 were treated by QM and the rest of the system by MM. The resulting systems contained around 5500 atoms with about 800 water molecules. The all-hydrogen potential function (PARAM27) [38] was used for the MM region, the self-consistence charge density functional tight binding (SCC-DFTB) [39,40] method was used for the QM region. The link-atom approach [41] was applied to separate the QM and MM regions.
The initial structures for the entire resulting system were optimized by use of the steepest descent (SD) and adopted-basis Newton-Raphson (ABNR) methods. The systems were gradually heated from 50 to 310.15 K in 50 ps. A 1 fs time step was applied to integrate equation of motion. 1.5 ns QM/MM MD simulations were initially executed for each of the reactant complex. The reaction coordinate was defined as a linear combination of r(CM ... Sδ) and r(CM ··· Nη1/2), which is R = r(CM ··· Sδ) − r(CM ··· Nη1/2). The umbrella sampling method [42] in the CHARMM program along with the weighted histogram analysis method (WHAM) [43] was used to determine the free energy(PMF) change as a function of the reaction coordinate(s). The Monte Carlo Simulation bootstrapping integrated with WHAM [43] was applied 100 times to estimate the errors (shown in Figure 2 and indicated in the legend of Table 1) of the PMF profile. 20-22 simulation windows were saved for each methyl transfer process. And for each window 50 ps production runs were performed after 50 ps equilibration. The force constants of the harmonic biasing potentials used in the PMF simulations were 50-500 kcal·mol −1 ·Å −2 . During the 50-ps production run of each window, structural coordinates were saved every 50 steps (1 fs/step), and consequently there are 1000 frames saved for each window. The distances shown in Figures 3, 4, 6 and 7 are the means of corresponding distances with stand deviation (S.D.) collected from those 1000 frames. The 2D free energy (PMF) map were also determined with umbrella sampling method and two-dimensional-WHAM. The time step for simulations is 1 fs. The horizontal reaction coordinate for the 1st methylation on Nη1, Rx = r(CM ... Sδ) − r(CM … Nη1), was used to describe the methyl transfer process. And the vertical reaction coordinate, Ry = r(HH12 ... Nη1) − r(HH12 … OE1), was used to describe the proton transfer process. Similarly, for the 2nd methylations, the methyl transfer was explained by Rx = r(CM ... Sδ) − r(CM … Nη2) and the proton transfer by Ry = r(HH22 ... Nη2) − r(HH22 … OE1). Approximated 300-400 windows were used in the construction of 2D free energy map, with 50 ps production run following 50 ps equilibration for each window. The force constants for each window were in the range of 100-800 kcal·mol −1 ·Å −2 .

Conclusions
In this study, the Type II PRMT5 was investigated by the QM/MM MD and free energy simulations. E435 was found to function as a general base catalyst for accepting a proton from the substrate arginine during the 1st and the 2nd methylation. This residue seems to play a similar role as E144 in Type I PRMT3. The simulations provide the detailed mechanism for the symmetric di-methylation by the enzyme and suggest the possible role for some other key residues as well during the catalysis. The kcat for the arginine methylation catalyzed by PRMT5 is measured to be between 10 to 50 h −1 , indicating that the activation energy barrier is about 18 kcal/mol based on transition state theory [23]. Thus, our estimates of the free energy barriers for the methyl transfers (19-20 and 18-19 kcal/mol for the first and the second methylation, respectively) are reasonable compared with the experimental observations. Furthermore, the simulations suggest that the symmetric di-methylation by PRMT5 seems to be energetically favorable, in agreement with the fact that PRMT5 is Type II PRMT. The proposed mechanism here based on the simulation results is different from the asymmetric dimethylation [25,26].