1. Introduction
The facilitated diffusion of ammonia across cellular membranes by the Amt/Mep/Rh family is key for acid-base homeostasis. The ability to transport ammonia by the four structures from the Rh/Mep/Amt superfamily was determined using x-ray crystallography of AmtB from Escherichia coli [1] [2], Amt-1 from Archaeoglobus fulgidus [3], Rh50 from Nitrosomonas europaea (Rh50NE) [4], and RhCG from humans [5]. The five human Rhesus proteins fall into two different groups: ammonia transporting glycoproteins (RhAG, RhBG, and RhCG) and the non-transporting Rh proteins (RhD and RhCE) [5]. RhBG and RhCG are found in a variety of tissues including the brain, liver, gut, and kidney. Epithelial expressions of RhBG and RhCG in the kidney are involved with ammonium secretion in the connecting tubules and collecting ducts [6]. While RhCG has been reported to be found in the basolateral and apical membranes; RhBG has not yet been detected in those areas [7].
The previous computational studies for the AmtB mechanism [8] showed that the ammonium initially comes in contact with the charged residues at the periplasmic loops of Escherichia coli. Then,
follows the negatively charged oxygen atoms of the side chain of the carboxyl groups and arrives at the aromatic matic cage, which is composecage at the bottom of the periplasmic vestibule. This arod of F103, F107 and W148, stabilizes the charged substrate at the entrance of the hydrophobic pore. AmtB from Escherichia coli has two tightly stacked phenylalanine rings at the periplasmic entrance of the pore. The synchronized opening and closing of these two phenylalanine rings allows charged ammonium to enter the pore and release its proton to the bulk of water, becoming a neutral NH3 in the transfer [9].
Many studies have been done to decipher if ammonia is transported into the cell by RhCG as
or as NH3. Measurement of apical plasma membrane
and NH3 permeability in the collecting duct of mice with global Rhcg deletion found decreased NH3 permeability without any change in
permeability [10]. Other studies reconstituted purified human RhCG into liposomes and demonstrated increased NH3 permeability but no differences in
permeability [5] [11]. Thus, the bulk of evidence suggests Rhcg/RhCG functions as a facilitated NH3 transporter [12] [13].
The structural differences between AmtB and Rh proteins suggest a variation on this Am conduction mechanism [5]. First of all, the RhCG loops are bigger and contain more negatively charged residues, which could possibly yield better conduction rate. In addition, the W148 is replaced by L170, and F103 is completely missing as it does not have comparable residue in RhCG. The absence of an aromatic cage at the bottom of the vestibule raises the question of recruitment site. At the bottom of the periplasmic vestibule, A144 of AmtB is replaced by a conserved E166 in RhCG, which looks like a prime candidate for a charged substrate recruitment site. The crystal structures of Rh50 from Nitrosomonas europaea and human RhCG also show that the two PHE rings at the periplasmic vestibule exit are not stacked in Rh proteins (see Figure 1). Additionally, the cytoplasmic exit of the pore is not guarded by a PHE ring as in E. Coli AmtB; Human RhCG and Rh50NE have LEU and CYS residues instead. These structural differences at the cytoplasmic exits may result in differences in water’s ability to access the central pores of AmtB and Rh proteins.
While this study focuses mainly on RhCG, we also included results from the simulations for Rh50NE to try and investigate general ammonia conduction mechanisms in Rh proteins.
2. Results and Discussion
2.1. E166 Recruits Charged Ammonia
Human RhCG and Rh50NE periplasmic vestibules do not have an aromatic cage [4] [5] [14] as E. coli AmtB does [2]. The main piece of the AmtB aromatic cage,
Figure 1. Overlapped view of significant residues in AmtB (red), human RhCG (yellow), and Rh50NE (blue).
W148, is replaced by a LEU in Rh proteins (L126 for Rh50NE, L170 for human RhCG). The proposed mechanism used by AmtB relies on recruitment of
by the aromatic cage [9], as the absence of W148 signals that Rh proteins use a different mechanism for ammonia recruitment than similar ammonia transporters. At the bottom of the periplasmic vestibule, A144 of AmtB is replaced by a conserved Glu in RhCG (E166) and Rh50NE (E122), which looks like a prime candidate for a charged substrate recruitment site.
Our equilibrium MD simulations showed that
is readily recruited by E166 on RhCG. The loops of RhCG were particularly effective for
recruitment. Within 5 ns of the equilibration MD simulations two of the
molecules were recruited by E166 from 10 Å away from the lipid phosphate group plane. This recruitment site attracted the molecules so strongly that
molecules did not leave this site for the next 40 ns of the MD simulations. In one of the observed recruitment cases, another
molecule was also observed to arrive at the E166 site, and stayed there for almost 20 ns. Our simulations clearly suggest that conserved E166, which is replaced by A144 in AmtB, is a
recruitment site for RhCG proteins.
We also verified the local energy minima status of E166 site for
by using umbrella sampling simulations. Umbrella sampling simulations also revealed an energy minimum for
at the phenyl gate. However, we did not observe
readily moving into the phenyl gate area at any point in 450 ns of a monomeric structure simulation. The energy barrier of around 4 kcal/mol between the E166 residue and the phenyl gate is significant and is likely the main reason that
does not spontaneously move to phenyl gate energy minima (Figure 2).
One reasonable claim might be that
is being recruited so strongly at the E166 site that it would prevent any other
molecules from entering the vestibule. However, during our MD simulations, we observed that is not the case. Although
is connected to E166 with hydrogen bonds, another
molecules come as close as 3.9 Å to the recruited
(Figure 3). This means that the effective charge of the site is not enough to repel other
molecules away from the recruitment site as one might expect; instead the charge of the GLE and the surrounding residues as well as the water molecules “screen” the
charge and allow other
molecules to approach.
Figure 2. Free energy of systems with
and NH3 as the molecules travel into RhCG protein with reference to E166 and F235, the bottom phenyl of the gate.
Figure 3. Two
molecules attracted to GLU-166 with a distance of 3.87 Å between them.
Although the equilibrium MD simulations have not observed an
arriving at the phenyl gate, NH3 has been observed to position itself comfortably there, shown in Figure 4. F130 and F235 create an excellent NH3 recruitment site for RhCG. At this position there is no water within 5 Å of NH3. The free energy profile derived from the umbrella sampling simulations with NH3 also supports this result, as the energy minima of NH3 lies just above the bottom phenyl of the gate (Figure 2).
Figure 5 shows the normalized probability distributions for NH3,
and CO2 molecules through the channel. This distribution, obtained from high density equilibrium simulations, clearly identifies the energy minima locations of E166 for
and PHE rings for NH3. During the long equilibration simulations, the CO2 molecules did not show any tendency to diffuse the periplasmic vestibule, same as the AmtB and Rh50NE simulations.
Figure 4. NH3 settled in phenyl rings with nearest water molecules.
Figure 5. RhCG probability distributions for three substrates, corrected with available space (radius) measured via hole. The z-position of primary residue(s) each substrate is attracted to have been labeled.
2.2. No Synchronized Phenyl Motion for Rh Proteins
Throughout the known structure of the Amt/Mep/Rh protein family the periplasmic vestibule is separated from the pore with two phenylalanine residues. For the Amt protein, the phenylalanine residues are stacked parallel on top of each other and completely cover the channel. Rh proteins, on the other hand, have phenylalanines that create a right angle with each other with one in a vertical position while the other lies horizontal. In their experimental study, Javelle et.al found out that mutating both phenyls kills the function of AmtB [8]. On the other hand, for human RhCG, single mutation of the phenyl pair substantially stops the NH3 conduction, but the double mutation keeps the channel almost as active as the WT [15]. The difference between these experimental studies suggests that the role of the phenyl gate is different between Rh proteins and AmtB.
The conserved twin phenyl residues that keep the gate of the hydrophobic pore closed are observed to flip open in a synchronized manner for E. Coli AmtB [9]. Additionally, the dihedral angle Ca-Cb-Cg-Cd1 distribution shows two energetically favorable configurations of AmtB phenyls. The π-π interaction between the phenyl side chains allows these two stable cases [9].
During over 1 μs MD simulations performed on both Rh50NE and human RhCG, we did not observe synchronized phenyl motion. This may suggest that the deprotonation process of Rh proteins is achieved through a different mechanism than that of AmtB.
The
recruitment location of the GLU is around 10 Å away from the phenyl gate; in our previous study we found that the phenyl gate actively takes part in the deprotonation of
. This is another indication that the Rh proteins have different conduction mechanisms for both
and NH3 than AmtB.
Our analysis on dihedral angles Ca-Cb-Cg-Cd1 shows that when there is no substrate in the vicinity of the phenyl gate, the F130 takes dihedral angles of 16 ± 23 or −170 ± 14 while F235 takes 87 ± 14 or −93 ± 12. These numbers represent the mean and the sigma of the Gaussian fit. Having exactly 180 degrees between these angles clearly shows that there is an only one energy minimum for phenyl pairs on human RhCG in the absence of substrate. However, the distribution of the same dihedral angles with the presence of
at the E166 recruitment site shows that the presence of
changes the dihedral angles F130 favors. When
is in the vicinity, F130 takes two dominant angles: −42 ± 14 and 76 ± 14 (Figure 6). The existence of charged substrates seems to change F130 dynamics while F235 is not affected.
Although the crystal structures of Rh50NE and human RhCG manifest similar positions for phenyl pairs, the dynamics of the ring create differences in position frequency. When there is no substrate in the vicinity of the phenyl pairs, the F86 of Rh50NE takes dihedral angles of −61 ± 16 and 60 ± 28 while F194 distribution yields −133 ± 26 and 66 ± 26. Having 180 degrees between the dihedral angles of F194 means the bottom of the phenyl pairs on both human RhCG and Rh50NE have only one energetically favorable position. On the other hand, the top phenyl of Rh50NE (F86) clearly shows two distinct favorable angles with almost 90 degrees between them.
Figure 6. RhCG dihedral angles for F130 (top) and F235 (bottom) (left column), RhCG F130 (top) and F235 (Bottom) dihedral angles with
is recruited by the protein.
When we extend the analysis to the cases when the substrate (
or NH3) is close to the phenyls, the top phenyl shows different angle preference. The F86 residue prefers to keep the Ca-Cb-Cg-Cd1 dihedral angle at −61 ± 16 minima, while F194 seems unaffected by the existence of the substrate.
2.3. Water Leaks into the Central Pore of Rh Proteins
The crystal structures of Rh50NE [4] [14] (1.3 Å resolution) and human RhCG (2.1 Å resolution) have some distinct pore differences when compared to E. coli AmtB. The conserved phenyl rings at the pore entrance from the periplasmic vestibule are not stacked as tightly as in AmtB. The cytoplasmic exits of the pores on these structures do not have a phenyl gate either. This allows water molecules to have hydrogen bonds with central histidine residues, thus shortening the length of the hydrophobic pore to 11 Å. By contrast, the distance between the water molecules on both sides of the pore for AmtB is more than 20 Å.
During long equilibrium MD simulations (equivalent of 540 ns simulations for a monomer) on a trimeric Rh50NE model, water molecules were observed to diffuse through pore in both directions. Although the F194 appears to keep the periplasmic entrance of Rh50NE pore closed, there is still enough space for water molecules to diffuse through. Throughout the simulations, a total of 16 water molecules leaked into the pore. Four of these water molecules diffused from the cytoplasmic vestibule to the periplasmic vestibule, and four others diffused in the opposite direction. The remaining 8 water molecules were observed to leak into pore from cytoplasm; after spending some time in pore they moved back to the cytoplasmic vestibule. The average time it took for a water molecule to diffuse from periplasm to cytoplasm was measured to be 2.4 ns. The cytoplasmic entrance of the pore is even less protected for Rh50NE proteins. The four water molecules conducted to the periplasmic vestibule of Rh50NE from cytoplasmic vestibule in average of 1.1 ns. However, the rest were seen going back to the cytoplasmic vestibule after spending an average of 1.5 ns within the pore. On two different occasions two water molecules entered together to Rh50NE pore.
The crystal structure of Rh50NE shows that the bottom of the periplasmic vestibule is divided into two parts by the I198 and F86. The narrow pathway created by A195, T196, A197, and I198 allows water molecules to create a water bridge between the pore entrance and the periplasmic vestibule. The water conduction from periplasm to pore is mediated through this narrow pore. Otherwise the water molecules located at the entrance of the pore have not been observed to be accessible from the F86 side of the periplasmic vestibule. Chérif-Zahar et al. [16] found that Rh50NE regulates NH3 uptake, but that it appears to lack the primary proposed ammonium binding site found in AmtB and thus seems to function better in a higher pH environment. A computational study by Hub et al. [17] also reports that the Rh50NE should be partially active and suggests that its water conduction rate is 40 times slower than its conduction rate of NH3. The 540 ns monomeric MD simulations reported here clearly observed that water molecules being conducted by Rh50NE.
In regards to the human RhCG model, Baday et al. [18] believed the pore to be hydrophobic because of a conserved phenylalanine. Our simulations show that, although the water molecules can enter the channel, the human RhCG does not conduct water molecules. After approximately 450 ns of monomeric MD simulations performed with a trimeric human RhCG model, we observed 23 water molecules entering the pore from the cytoplasmic vestibule. All of these water molecules leaked into the pore from the cytoplasmic vestibule and returned back to cytoplasm. During the water leakages to the pore, the water molecule(s) tended to stay in the pore for 5 ns or less (Figure 7). As is shown by the crystal structure, there is a water molecule at the top of the pore that has hydrogen bonds with H185, G179, G290 and N236 and does not leave. Unlike Rh50NE, this water molecule positioned at the pore and periplasm interface does not have any direct connection to the bulk of the water in the periplasmic vestibule. However, the cytoplasmic entrance of the human RhCG pore allows water molecules to leak through more freely. On many occasions during 450 ns monomeric MD simulations, two water molecules were observed to enter the pore together. In these cases, the water molecules were more stable inside the pore. Interestingly, in one occasion we observed a transiently stable chain of water molecules flooding the pore from cytoplasm all the way to F235 (Figure 8).
3. Conclusions and Discussion
In an attempt to better understand the mechanism of transport in RhCG transporters, we performed various simulations that affected different molecules in the system to better identify how they interact.
Although the transport mechanism is still a mystery, our simulations suggest some viable theories as to what the process is. Because
is attracted to E166 and NH3 is attracted to the phenylalanines, it is logical to assume that there is a place of deprotonation somewhere between the two points. Another
Figure 7. Positions of the water molecules interacting with the RhCG pores over time. The two lines represent F235 (left) and H334 (right).
Figure 8. Water chain observed in the periplasmic pore of RHCG. The water molecules enter from cytoplasm but they cannot pass to periplasm. W490 was detected at this location on the x-ray structure of the protein [5].
molecule has been observed to come very close to the
molecule already attracted to the E166. Since there are no proton acceptors around it, this suggests the deprotonation happens before the rings, possibly between E166 and the phenyl gate. The new molecule may serve as a catalyst to kick out the already attracted
molecule, causing it to be forced over the energy barrier.
Because of the large barrier of around 4 kcal/mol beyond E166, it is reasonable to believe that the
molecule loses a proton before making the jump, considering it would be more energetically favorable to break through the barrier in an uncharged state. It would be very difficult for it to get over the barrier without doing so. After it loses the proton and becomes NH3, it is attracted to the phenylalanines at the top of the pore.
The phenylalanines are constantly moving and turning in the simulations. This might give the NH3 molecule the opportunity to slip into the channel and navigate to the pocket with the structural water molecule. Once in the channel, we know that the NH3 will act as the water does because it is a similar size, and will jump down past the histidines into the periplasm (Figure 9). During our simulations we substituted a NH3 molecule entering the pore in place of the second water molecule, and observed that the NH3 acted similarly to water and had jumped down the pore to the cytoplasmic side.
Our simulations have shown that the pore of RhCG is much less hydrophobic than that of AmtB. During our long equilibrium MD simulations (the equivalent of 740 ns for a single monomer) in an earlier study on trimeric AmtB, the author of this study reported only one water molecule entering the hydrophobic pore from cytoplasmic vestibule [9], unlike another study reporting hydrophilic AmtB pore [19]. In this study, we found that both the Rh50NE and human RhCG crystal structures had water molecules at the cytoplasmic exits of the pore. Furthermore, Zidi-Yahiaoui et al. report water conduction by human RhCG in HEK293E cells [14]. The reported conduction rate is ten times less than AQP1 and does not change with F130A/F235V double mutation. However, Gruswitz et al. reports that stopped flow assays conducted on proteoliposomes do not show any water conduction by human RhCG. During our 450 ns of monomeric MD simulations we did not see any water conduction from one vestibule to the other. These results are in agreement with Gruswitz et al. [5].
Further simulations and investigation will be necessary to illustrate the complete transport mechanism in RhCG. Nonetheless, our simulations suggest that
is brought in and is converted into NH3, which then slips into the pore and diffuses in the same way the water molecules do. Figure 10 summarizes our proposed mechanism of human RhCG ammonia conduction.
Figure 9. Our simulations showed NH3 acting like a water molecule when we placed it in the channel. On the left, it is stuck in the pore with the first water molecule and on the right it is shown to move down the histidine complex as expected.
Figure 10. Proposed mechanism of
deprotonation and NH3 conduction by RhCG. (a)
is recruited by E166; (b) A second
molecule is attracted to E166, pushing the first
towards the phenyl rings; (c) The first
reaches the energy barrier, sheds a proton to become NH3 and settles into the phenyl rings at its energy minimum; (d) As the rings fluctuate, the NH3 works its way into the pore, and settles near the water shown in the crystal structure in the mildly unstable pocket; (e) NH3 acts in a similar manner to water, and rushes out of the pore down the channel; (f) NH3 is flushed out of the channel.
4. Materials and Methods
Model building. The trimeric human RhCG model was based upon the X-ray structure of RhCG at 2.1 Å resolution (PDB ID = 3HD6) [5]. The periplasmic loops with 18 residues (residue IDs 35 - 52, and 22 residues (residue IDs 362 - 283) were constructed by using COOT [20] [21]. The trimeric Rh50 from Nitrosomonas europaea model was based upon the X-ray structure with 1.3 Å resolution (PDB ID = 3B9W) [14]. The trimeric AmtB model was based upon the X-ray structure of AmtB at 1.35 Å resolution (PDB ID = 1U7G) [2]. Three mutations (F68S, S126P, and K255L) and SeMet residues were changed back to their native residues.
Hydrogen atoms were added to the structures with XPLOR-NIH [22]. The protonation states for each His residue were modeled based on their hydrogen bond environment. Phosphatidylethanolamine (POPE) was used to create a model of a lipid bilayer with dimensions of 120 Å × 120 Å × 39 Å. Trimeric protein was inserted into the center of the lipid bilayer and overlapping POPE molecules were removed. The z-axis is aligned in the direction of the hydrophobic pore of the protein and normal to the lipid bilayer plane. The lipid-protein complex was solvated by adding 25 Å TIP3P [23] water slabs on each side of the lipid bilayer. The model used in all simulations contains 22,718 water molecules and more than 142,000 atoms, with dimensions of 120 Å × 120 Å × 86 Å.
All the MD simulations were performed with the NAMD program [24] using the CHARMM36 force field [25] [26]. Periodic boundary conditions were applied in all three dimensions. A Langevin thermostat and the Langevin piston model [27] were used to keep the temperature and pressure constant at 310 K and 1 atm, respectively. The Van der Waals interaction cutoff was set at 12 Å in conjunction with switching at 8 Å. Full electrostatics was employed using the Particle Mesh Ewald (PME) method [28], with a grid spacing less than 1 Å. All simulations used an integration time step of 1.0 fs. The total charge of the system was neutralized by adding Cl− and Na+ ions.
Equilibrium molecular dynamic simulations. To investigate the water diffusion into protein vestibules, separate initial models were built for each protein. Before the long MD simulations, each system was first minimized for 2000 steps and then equilibrated for 2 ns with no constraints on the system. Total of 246 ns, 166 ns, and 164 ns of MD simulations were performed on AmtB, Rh50NE, and RhCG, respectively.
To determine how water diffuses into these proteins, the probability distributions were calculated, by tracing the 22,718 water molecules inside a hypothetical cylinder, perpendicular to the plane of the lipid bilayer. The height of the cylinder extends the thickness of the periodic boundary box, and with 70 Å diameter it accommodates the trimeric protein. The number of solute molecules that appeared inside the cylindrical volume was recorded at every 5 ps frames by monitoring the location of the oxygen atom of each water during the MD simulations. The probability distribution curves were derived from the observed number of solute molecules per ns × Å3 along the axis of the cylinder. The probability distributions were corrected with respect to the radius of the channel at that z value. The radius of the channels was determined by using HOLE program. For each protein, the total probability distributions within the cylinder were normalized to 100.
Umbrella sampling simulations. The umbrella sampling methods were used to calculate the free energy profile of H2O conduction through proteins. During umbrella sampling simulations, each monomer was sampled simultaneously with a 0.5 Å window size. In each sampling window atom coordinates were recorded at every 0.5 ps, for 450 ps. The first 150 ps were used for equilibration and the last 300 ps were defined as the production run. The starting point for each window was determined by the positions of the atoms from the last frame of the previous window. This way sufficient relaxation of the system was assured. Å biasing harmonic potential with force constants of 10 kcal/(mol × Å2) was applied to the O of the H2O molecule. The solute molecules were restrained only on the z axis and free to move on the xy plane. The data were unbiased and a free energy profile was calculated by using Weighted Histogram Analysis (WHAM) method [29] [30].
Acknowledgements
The authors would like to thank XSEDE for computational time, Coe College, and the McElroy Foundation for continuous support. This work has been supported by NSF grants PHY1950337 and DMR1746230.