Mechanistic insights into multiple-step transport of mitochondrial ADP/ATP carrier.

The ADP/ATP carrier (AAC) is crucial for mitochondrial functions by importing ADP and exporting ATP across the inner mitochondrial membrane. However, the mechanism of highly specific ADP recognition and transport by AAC remains largely elusive. In this work, spontaneous ADP binding process to the ground c-state AAC was investigated through rigorous molecular dynamics simulations of over 31 microseconds in total. With improved simulation strategy, we have successfully identified a highly specific ADP binding site in the upper region of the cavity, and this site exhibits selectivity for ADP over ATP based on free-energy calculations. Sequence analyses on adenine nucleotide transporters also suggest that this subgroup uses the upper region of the cavity, rather than the previously proposed central binding site located at the bottom of the cavity to discriminate their substrates. Identification of the new site unveils the unusually high substrate specificity of AAC and explains the dependence of transport on the flexibility between anti and syn glycosidic conformers of ADP. Moreover, this new site together with the central site supports early biochemical findings. In light of these early findings, our simulations described a multi-step model in which the carrier uses different sites for substrate attraction, recognition and conformational transition. These results provide new insights into the transport mechanism of AAC and other adenine nucleotide transporters.

Introduction

Mitochondria are the powerhouse of eukaryotic cells that synthesize ATP through oxidative phosphorylation (OXPHOS) [1], [2]. The generated ATP in the mitochondrial matrix is exported to the cytosol where it is hydrolyzed to ADP to fuel cellular processes, and ADP in the cytosol is imported back to the matrix for ATP regeneration. Both ATP export and ADP import are carried out by mitochondrial ADP/ATP carrier (AAC, or ANT for adenine nucleotide translocase), the most abundant protein in inner mitochondrial membrane [3]. AAC mediates ADP/ATP exchange through alternating between cytosolic-open (c-) state and matrix-open (m-) state (Fig. 1A). The ground c-state AAC specifically recognizes and imports ADP, and m-state AAC specifically recognizes and exports ATP [4]. However, mechanism of transport process of AAC remains largely elusive.

Fig. 1

A scheme of AAC function and the workflows of this study. (A) A schematic diagram to show the function of AAC. The structures of AAC in c-state (PDB entry: 1OKC) and m-state (PDB entry: 6GCI) are shown. Residues of the previously proposed central binding site (K22, R79, G182, I183 and R279) are highlighted in blue spheres. (B) Workflow of the current work. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

AAC is a paradigm of mitochondrial carrier family (MCF), the biggest solute carrier subfamily that composes of 53 members in human. MCF transporters are special for their small size and for having three homologous domains instead of two. Each domain is about 100 aa long and contains highly conserved MCF motif: Pxx[D/E]xx[K/R]xR-(20–30 residues)-[D/E]Gxxxx[Y/W/F][K/R]G [5], [6], and these motif residues determine the common structural scaffold of this family. Despite similar structures of these transporters, their substrates are extremely diverse, ranging from very small proton to very large Coenzyme A (CoA) and NAD+. Therefore, members of this family are intriguing subjects to identify structural determinants for specific and selective transport, which are of special significance for rational drug design.

AAC is the only MCF member whose atomic-level crystal structures have been solved up to date [7], [8], [9]. The structure of c-state AAC exhibits three-fold pseudo-symmetry, and each homologous domain consists of two transmembrane (TM) helices connected by one short matrix helix. The six TM helices twist into a right-handed coiled coil and surround a cone-shaped transport cavity (Fig. 1A) [7]. From the crystal structure it is hard to infer the ADP binding sites because the co-crystallized inhibitor carboxyatractyloside (CATR) is much bulkier and structurally different from ADP. Analysis of the structure suggested that the featured tyrosine ladder (Y186 Y190 Y194) and four basic patches from entrance to the bottom of the cavity (K198R104, K91K95R187, K22R79R279 and K32R137R234R235 respectively) could be important in ADP binding [7], [10]. Based on the crystal structure, theoretical studies including molecular modeling [11], symmetry analysis [12] and molecular dynamics (MD) simulations [13], [14] identified a central substrate binding site located at the bottom of the cavity. It was further proposed that this central site is the single substrate-binding site of AAC and all the other MCF members [9], [15], [16]. However, early biochemical experiments from different groups consistently demonstrated that AAC has two specific substrate binding sites: one high-affinity site and one low-affinity site [17], [18], [19], [20], and a pair of high-affinity and low-affinity sites interacts with one CATR [21]. Moreover, the transport process was characterized as a two-step event: nucleotide binding to a specific site, followed by the vectorial process of transport [21], and the transport is a slow process with one ADP/ATP exchange cycle estimated to be 12–15 ms (turnover number around 2000–2500/min at 37 ℃) [22]. It’s apparent that discrepancies exist in previous studies on the number of ADP binding sites and the kinetics of the ADP binding process. Therefore, more structural details, especially the structural dynamics information are needed to clarify the ADP transport process of AAC.

Starting from the static crystal structures, all-atom MD simulation is currently the only technique to generate dynamic trajectories for proteins at atomic level. It’s especially powerful in elucidating structure–function relationship for those functionally important model molecules like AAC, for which a huge body of previous experimental results can serve as validations and provide guidance to better interpret the simulation results. Previously we have successfully used MD simulations to improve our understandings on the conserved MCF motif and cardiolipin binding sites of AAC [6], [23]. In the current work, spontaneous ADP binding process to c-state bovine AAC1 was investigated through rigorous MD simulations of over 31 μs in total. Improved simulation strategies enabled us to describe a multi-step process of ADP binding to AAC. Most significantly, a new highly specific ADP binding site was identified in the upper region of the cavity. This new site was also supported by sequence analysis among the five human adenine nucleotide transporters. The workflow of the current work is illustrated in Fig. 1B.

Materials and methods

MD simulations

MD simulations on the c-state apo-AAC

Three independent 3-μs all-atom MD simulations (traj-1 ∼ 3) were first carried out to sample relaxed conformations of apo-AAC. The system is composed of 70,769 atoms, including one AAC molecule, 219 POPC lipids, 12,226 water molecules, 23Na+ and 42Cl-. The GROMACS 4.5.5 package [24] is used for MD simulations. Please refer to our previous work for more details on the system setup and simulation protocols [6]. The same trajectories have been used to analyze interactions among MCF motif residues [6].

Large scale short MD simulations on ADP binding

Snapshots were taken from 2 μs till the end of the above three apo-AAC trajectories with an interval of 100 ns. This leads to total 33 snapshots of different conformations, and for each snapshot an ADP was put near the entrance of the carrier in a random orientation and the waters or ions overlapped with the added ADP were removed. The GROMACS 4.5.5 package [24] is used for MD simulations, with the force field CHARMM36 [25] applied to the protein and refined CHARMM lipid parameters for the POPC molecules [26]. As each system was already fully equilibrated except the solvent around the added ADP, we ran simple equilibration steps. Two rounds of minimizations were first carried out, with and without positional restraints on the protein heavy atoms. Then each system was heated from 50 K to 310 K with positional restraints applied on Cα atoms of the protein in a time step of 2 fs, and simulation time is 50 ps. The equilibration was followed by a 10-ns production MD simulation in an NTP ensemble (310 K, 1 bar). The temperature of the system was maintained through the v-rescale method [27] with a coupling time of 0.1 ps, and the pressure was maintained using the Berendsen barostat [28] with τ of 1.0 ps and compressibility of 4.5 × 10-5 bar−1. The SETTLE [29] and LINCS [30] constraints were applied to the hydrogen-involved covalent bonds of water and other molecules respectively. The Particle-Mesh Ewald (PME) algorithm [31] was applied for electrostatic interactions calculation. Minimization and equilibration were carried out in local workstation carrying 32 CPUs (Intel Xeon E5-2650 2.00 GHz), and production simulations were carried out on TianHe-1(A) at National Supercomputer Center in Tianjin, China. The coordinates for each system were saved every 10 ps.

Microsecond-timescale simulations on ADP binding and ATP binding

Six independent MD simulations were carried out on ADP binding (ADP-1 ∼ ADP-6). For comparison purpose, six simulations were also carried out on ATP binding (ATP-1 ∼ ATP-6). The twelve simulations started from the same conformation of the carrier taken from the last snapshot of one trajectory in our previous work [6], but with different orientations in ADP or ATP. Remaining simulation protocols and parameters are equal to the ones employed in the shorter simulations, as described above. When we ran ADP-1 and ATP-1, the system corrupted and therefore only ADP-2 ∼ ADP-6 and ATP-2 ∼ ATP-6 are included in further analysis in this work.

MD simulation to promote ADP movement from the specific site S2 toward the central site S1

The simulation (ADP-mod) started from the last snapshot of ADP-3, with each atomic charge on the diphosphate group of ADP reduced to 70% of original values in the ADP topology file (the.itp file in Gromacs) and the charges of other atoms kept as the original, and therefore the total charge of ADP was reduced to −2.3. Minimization and equilibrations were first carried out with restraints placed on both ADP and AAC, followed by a 2.1-μs production simulation. Then we continued the simulation from the checkpoint file (.cpt file) till 3.2 μs, but used the original topology file (.tpr file) in which the charges of the atoms constituting ADP were set to the original values. The simulation was performed using GROMACS 4.5.5 with a time step of 2 fs. Other protocols are the same as above.

Trajectory analyses and free energy calculation

Most analyses were carried out with programs provided in the GROMACS package. Trajectories were viewed with VMD [32] and structural graphics were produced by the PyMOL Molecular Graphic system (version 2.0 Schrödinger, LLC). Time evolutions of the pocket volume of the c-state apo-AAC were calculated with POVME [33]. The line plots, edge bundling plot and free energy landscape were created in R studio with ggplot2 and ggraph packages. Binding free energies were calculated with the MM/PBSA method of AMBER 14 package [34], with the GROMACS trajectories (.xtc files) converted to the AMBER format (.nc files) through the mdconvert command in the MDTraj toolkit and the GROMACS topology files converted to the AMBER format through ParmEd (https://github.com/ParmEd/ParmEd).

Sequence analyses among adenine nucleotide transporters

Multiple sequence alignment was carried out on amino acid sequences of all 53 human mitochondrial carriers obtained from UniProt, and then the maximum-likelihood phylogeny tree was constructed in Mega 6, with support for the nodes calculated using Bootstrap method for 500 replications. For each of the carriers AAC, SLC25A42, GDC (Graves disease carrier protein, or SLC25A16), SCaMCs (short Ca2+-binding mitochondrial carrier, or Mg-ATP/Pi carrier) and SLC25A41 that fall into the same clade with AAC in phylogeny analysis, multiple sequence alignments were carried out on 41 reviewed sequences of AACs from 23 organisms, 231 SLC25A42 sequences from 146 organisms, 185 GDC sequences from 130 organisms, 450 SCaMC sequences from 131 organisms and 84 SLC25A41 sequences from 58 organisms respectively. For the carriers excluding AACs, unreviewed sequences were also included for sequence alignment, as the number of reviewed sequences is limited. Then based on the multiple sequence alignment results, the sequence logos were created through WebLogo [35].

Results

Microsecond-timescale MD simulations provide more relaxed conformations of the c-state apo-AAC

As the inhibitor CATR is much bulkier than ADP, crystal structure of CATR-bound AAC was proposed to open much wider than the apo conformations [36]. Therefore, before working on ADP binding process, we first ran long time MD simulations (traj-1 ∼ 3, each 3 μs in length) to sample more relaxed conformations of the c-state apo -AAC. Results show that the pocket entrances become narrowed (Fig. 2A) and the pocket volumes decreas in all the three independent trajectories (Fig. 2B). Principal component analysis (PCA) shows that the most significant movements are observed in cytoplasmic loop 2 (CL2) and cytoplasmic halves of H1 and H6 (Fig. 2C).

Fig. 2

Relaxed conformations of c-state apo-AAC are obtained from three 3-μs MD simulations. (A) Different shapes of the pocket entrance in the crystal structure and in the structure at the end of the simulation traj-1. (B) Time evolution of the pocket volume of AAC in the three trajectories. (C) The first component of the PCA analysis on the three trajectories. (D) The difference of the salt bridges near the entrance of AAC in the crystal structure and apo-AAC simulations. Occupancies of the salt bridges are averaged over three trajectories and are shown in red numbers (in percentage). The occupancies of these salt bridges in each trajectory are listed in Table S1.

The salt bridges near the entrance of apo-AAC also show big differences from those in the CATR-bound crystal structure (Fig. 2D). In the crystal structure of bAAC1, only one salt bridge (D195:K198) is observed near the entrance of the cavity. While in apo AAC, D195 forms a highly stable salt bridge with R104 (occupancy: 97%), and K198 is free majority of the time. Although the salt bridge D195:R104 is not observed in the crystal structure of bAAC1, it appears in the crystal structures of the yeast AACs (D212:K119 in yAAC2 and D201:K108 in yAAC3). This justifies our simulation results and implies sensitivity of salt bridges at the lipid-water interface to the dehydrated crystallization conditions.

Although the relaxed conformations of apo AAC deviate from the CATR-bound crystal structure, the three trajectories show considerable variations in pocket size (Fig. 2B), major movements (Fig. 2C) and salt bridge distributions near the entrance (Table S1). This highlights the conformational heterogeneity of apo AAC and suggests a necessity to start from a variety of initial AAC conformations for ADP binding simulations.

ADP is first attracted to the first (K198K205) and second (K91K95R187) basic patches

Large-scale short MD simulations (S1 ∼ S33, each 10 ns long) were then carried out to reveal initial ADP binding steps in an unbiased way. The initial conformations of AAC were picked evenly from the last 1 μs of the three apo-AAC trajectories, with ADP placed in a random orientation above the pocket entrance in each simulation. Results show that in 30 out of 33 simulations, ADP was quickly (within 1 ns on average) attracted to the positively charged residues in the upper region of the cavity, with the highest first-contact frequencies observed on residues K205, K198, K95 and K91. Once in contact, none of ADP left the carrier or moved further toward the bottom of the cavity, but the bound ADP moved dynamically between these positively charged residues and two major moving trends were observed: K91/K95 → K91K95 → K91K95R187 and K198/K205 → K198K205 (Fig. 3A and Table S2). Calculations over all the 33 trajectories also show that K198K205 and the second basic patch K91K95R187 exhibit the highest ADP binding occupancies (Fig. 3B). K198R104 was reported as the first basic patch based on the crystal structure of c-state AAC, while in our simulations R104 formed a very strong salt bridge with D195 and did not get involved in ADP binding (Fig. 3B). Formation of strong R104:D195 salt bridge and its potential significance was also reported in previous MD simulation on much shorter time scale [47]. Therefore, these results suggest that the first basic patch is composed of K198K205 instead of K198R104.

Fig. 3

Initial ADP binding revealed by large-scale short MD simulations. (A) The movement of ADP phosphate moiety between the positively charged residues near the entrance of carrier in 31 10-ns MD simulations (in S19 and S26 ADP left the carrier in the beginning). The moving direction of each step is shown in gradient line from blue to red. The number of the simulations in which the first contact is observed is provided within the parentheses. (B) The ADP binding occupancies calculated on the 33 10-ns simulation trajectories are shown in colored scale on the 3D structure of c-state AAC. All positively charged residues are highlighted in spheres. (C) The trajectory of the ADP α-phosphate group in one 200-ns extended simulation (e-S33). The simulation time is shown in colored spectrum along the moving trajectory. Initial and ending positions of the ADP α-phosphate group are highlighted in blue and red spheres respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Based on the simulations, ADP bound to the first basic patch was more dynamic than ADP bound to the second basic patch, and they exhibited a moving trend from the first patch toward the second patch (Fig. 3A). To confirm this, we picked five simulations in which the ADP phosphate moiety bound to the first basic patch at the end of the 10-ns simulations (S1, S15, S24, S29, S33), and extended each of them till 200 ns. The results support our speculations: ADP in all the five extended simulations left the first basic patch (at 10 ns, 105 ns, 76 ns, 176 ns and 15 ns respectively) and moved on to bind with the second basic patch (Fig. S1). The above 10-ns simulations together with 200-ns extended simulations suggest initial steps of ADP binding: when approaching from intermembrane space, ADP is first attracted by the first basic patch (K198K205) and then quickly relayed to the second basic patch (K91K95R187), and this process is primarily driven by the non-specific electrostatic interactions between basic patches and phosphate moiety of ADP.

ADP attached to the second basic patch converges to a highly specific binding mode (site S2)

To understand the subsequent transport process of ADP, five long-time MD simulations (ADP-2 ∼ ADP-6, ranging from 800 ns to 6 μs) were carried out, starting from the same relaxed conformation of apo-AAC but different orientations of ADP. In ADP-2, the phosphate moiety of ADP first bound with K294 at 1.1 ns, and then K91 joined in and brought the phosphate moiety to the K91K95 patch at 2.3 ns. R187 joined in at 17.3 ns, from when the phosphate moiety was anchored stably to the second basic patch (K91K95R187) till the end of the 6-μs simulation. In ADP-3, the phosphate moiety of ADP became anchored to K205K294 within 1 ns, and then bound with K91K294 from 395 ns. From 610 ns, the phosphate moiety was relayed to the K91K95R187 patch via K91, and it bound stably with the second patch till the end of 3-μs simulation.

In ADP-4, ADP quickly bound to the K91K95 patch (within 3 ns), and from 340 ns, K294 joined in and bound with ADP together with K91 and K95. The adenine base was very flexible in the aromatic ladder region (Y186Y190Y194 and F191) until 686 ns when it stacked with Y194 and its –NH2 group formed electrostatic interaction with the backbone carboxyl group of S212 (Fig. S3A). This conformation promoted ADP to leave K91 and K95 (at 690 ns), and K294 brought the phosphate moiety of ADP to bind together with K205. From 730 ns, K198 started to join in and meanwhile the hydroxyl group of S212 flipped to bind with –NH2 group of the ADP, and this conformation (Fig. S3B) was stably maintained till the end of 2-μs simulations. In ADP-5, ADP first bound with K95 at 3 ns, with K91 joining in at 15 ns. From 60 ns till the end of the 800-ns simulation, the phosphate moiety was anchored to K91K95K294, and the adenine base was quite flexible (Fig. S3C). In the simulation ADP-6, the phosphate moiety of ADP quickly bound with K91K95 patch within 2.2 ns, while the adenine base stacked with Y194. Then the base moved to stack stably with both R187 and Y190 from 7.7 ns, and R187 also joined in to bind with the phosphate moiety at 502 ns (Fig. S3F). This conformation was well maintained till 1730 ns, when the adenine base started to flip around and the diphosphate group became anchored to K91K95K294 (Fig. S3G).

Therefore, consistent with the above large-scale short MD simulations, ADPs in these long-time simulations are also quickly attracted to the positively charged residues in the upper region of cavity with none of them leaving the carrier or moving further toward the bottom of the cavity. Long residence time on the second basic patch also highlights its strong electrostatic attraction to ADP (Fig. S2).Moreover, these long-time simulations also indicate the importance of stacking interactions between adenine base of ADP with the aromatic ladder (Y186Y190Y194 and F191).

Of special significance, a highly specific ADP binding mode was consistently observed in two independent simulations (ADP-2 and ADP-3): the diphosphate group of ADP anchors to the second basic patch, and the adenine base forms an exquisite H-bonds network with both N115 and the backbone carbonyl oxygen of R187 (Fig. 4A). Meanwhile, the adenine base forms π-π stacking with Y190 side chain and π-cation stacking with R187 guanidinium group, leading to a sandwich stacking structure (Fig. 4A). In addition, 2′-OH in the ribose ring occasionally (occupancy: 9%) forms H-bond with Y190 side chain. Therefore, the five residues (K91, K95, N115, Y190 and R187) form extensive and highly specific interactions with ADP, and every part of ADP is recognized. For convenience, this newly identified specific site is denoted as site S2, and the previously proposed central binding site at the bottom of the cavity is denoted as site S1 hereafter in this work, although residues involved in site S1 reported in the previous studies show some differences [16], [13].

Fig. 4

Long-time MD simulations identify a highly specific ADP binding site S2. (A) The highly specific ADP binding mode with site S2. Interaction occupancies are shown in red numbers. The favorable stacking structure is highlighted in the inset, with the bound sodium ion shown with a green sphere. (B) The consensus pathway of ADP locating to the specific binding mode at site S2 in ADP-2 and ADP-3. (C) Free energy landscape built from the χ angle and the z-component of the N6 Cartesian coordinates of ADP in ADP-2. (D) The non-specific ATP binding to the site S2. The favorable stacking to Y194 is highlighted in the inset, with the bound sodium ion shown in a green sphere. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Consensus pathway of ADP converging to the specific binding mode at site S2

Then we analyzed time evolutions of important parameters during ADP binding process. Very intriguingly, ADP in ADP-2 and ADP-3 converged to the same specific binding mode with site S2 via the same pathway (Fig. 4B). The phosphate moiety of ADP was first attracted and anchored strongly to the second basic patch (Fig. S4A-C, S5A-C), while the adenine base formed dynamic stacking interactions with the tyrosine ladder (Y186 Y190 Y194), F191 and the guanidinium group of R187 (Fig. S4D, S5D). These stacking partners help the adenine base of ADP sample enormous conformations and facilitate it converging to the specific binding mode. Of special significance, the two simulations converge to the specific binding mode via the same intermediate conformation. In this intermediate conformation, ADP adopts the relaxed anti glycosidic conformer and the adenine base forms T-shape stacking interactions with both tyrosine ladder and F191 (Fig. 4B middle). From this intermediate conformation, ADP in the two trajectories underwent the same changes to form the specific binding mode with site S2 (supplementary Movie S1 and Movie S2): the adenine base first experienced a rotamer change around the glycosidic bond, which led to the unfavorable syn glycosidic conformation and brought the partially negative O4′ and N3 atoms to the phosphate moiety. Then a sodium ion quickly moved in to stably bind with these negative atoms. With the glycosidic conformation fixed in the syn conformation, ADP quickly formed the specific H-bond interactions with R187 and N115 after subtle adjustment. Time evolutions of important parameters in ADP-4 ∼ ADP-6 are provided in supplementary Fig. S6.

In the free energy landscape built from χ angle and z-component of coordinates of N6 atom in ADP, the specifically bound conformation forms the deepest basin centered at ( = -47°, z(N6) = 5.13 nm) and the intermediate conformation form the second free-energy minimum centered at ( = -128°, z(N6) = 6.03 nm) (Fig. 4C). A significant population of the intermediate conformation and the consistency within ADP-2 and ADP-3 suggest that the conformational changes described above might represent a predominant pathway of ADP locating to the specific mode at site S2. Although Y186, F191 and Y194 are not directly involved in specific binding with ADP (Fig. 4A), these resides are crucial in the pathway of ADP locating to the specific mode and therefore are also important part of site S2. As the intermediate conformation and specific binding mode do not reach equilibrium in the limited simulation time, free energy landscape cannot be used to compare energetic difference of the two conformations. Here we use the MM-PBSA method. Although MM-PBSA usually gets lower values than the real binding free energy since the method does not include the (unfavorable) entropy contribution to binding, this method is effective to compare two different binding modes or different ligands binding to the same site [23], [37], [38]. Our MM-PBSA calculations show that the binding affinity in the specific binding mode (ΔGbind = -65.1 kcal/mol) is more favorable than that in the intermediate conformation (ΔGbind = -63.4 kcal/mol, P = 0.01) (Table 1).

Table 1

Binding free energy (kcal/mol) of ADP and ATP with the new site estimated by the MM-PBSA method.

Trajectory	Site	ADP conformation	Region used for calculation	VDWAALS	EEL	EPB	ENPOLAR	ΔGgas	ΔGsol	ΔGbinda
ADP-3	S2	Intermediate,anti	1687–1770 ns	−20.4 ± 0.4 b	−1273.5 ± 2.3	1234.8 ± 2.2	−4.4 ± 0.0	−1293.9 ± 2.2	1230.4 ± 2.2	−63.4 ± 0.6
ADP-3	S2	Specific,syn	2350–3000 ns	−16.9 ± 0.2	−1288.9 ± 2.2	1244.7 ± 2.0	−3.9 ± 0.0	−1305.8 ± 2.2	1240.7 ± 2.0	−65.1 ± 0.3
ADP-mode	S1	anti	2350–3265 ns	−27.6 ± 0.2	−1516.6 ± 2.1	1452.4 ±1.6	−5.2 ±0.0	−1544.3 ± 2.1	1447.2 ± 1.6	−97.1 ± 0.6
ATP-6	S2	syn	1000–2000 ns	−15.3 ± 0.2	−1506.6 ± 2.3	1467.5 ± 2.1	−4.3 ±0.0	−1521.9 ± 2.3	1463.2 ± 2.1	−58.7 ± 0.3

VDWAALS van der Waals energy; EEL Electrostatic energy; EPB Polar solvation energy; ENPOLAR Non-polar solvation energy; ΔGgas Absolute binding free energy in gas phase; ΔGsol Absolute binding free energy of solvation; ΔGbind Absolute binding free energy.

Here the binding free energy refers to the difference between the energy of the complex and the total energy of free carrier and free ligand, i.e. E(Complex) - E(Ligand + Carrier).

Errors represent standard errors of mean.

To demonstrate the role of the stably bound sodium ion (Fig. 4A), we carried out two additional simulations starting from the last snapshot of ADP-3. In one simulation we removed all sodium ions from the system, and ADP quickly (within 1 ns) lost its specifically bound conformation (Fig. S7A-F). In another simulation we heated the system to 368 K, and the specifically bound sodium ion left ADP at 330 ns and then the adenine base quickly lost the specific binding mode(Fig. S7G-L). These results highlight the importance of the stably bound sodium ion in maintaining the syn glycosidic conformation and the specific ADP binding mode at site S2.

ATP binds with site S2 in a non-specific way with lower binding affinity

ATP can also be imported by AAC, but at much lower efficiency. To understand mechanism of strong selectivity of c-state AAC for ADP over ATP, we also ran five MD simulations on spontaneous ATP binding (ATP-2 ∼ ATP-6, lasting from 200 ns to 2 μs). The results show that in the simulation ATP-2, the phosphate moiety of ATP first bound with K95 at 1.3 ns, ad then K91 joined in at 3.6 ns and R187 joined in at 7.4 ns. The phosphate moiety kept anchored to K91K95R187 till 77 ns, and during this period of time the adenine base pointed upward and fliped around but did not stack with tyrosine ladder or any other aromatic residue (Fig. S8A). At 77 ns, the phosphate moiety left R187 and bound with K91K95K205K294 (Fig. S8B). At 136 ns, K205 and K294 left ATP, and the phosphate moiety kept anchored to K91K95 till the end of the 200 ns simulation. In ATP-3, the phosphate moiety first bound with K91K95 at 0.6 ns, then it moved dynamically between K91, K95 and K294 till 154 ns when ATP left the carrier and moved into the bulk solvent. In ATP-4, the phosphate moiety first bound with K95 at 0.5 ns, and then K91 joined in at 6 ns and K294 joined in at 28 ns. From 80 ns, the adenine base stacked with Y194 (Fig. S8C). At 236 ns, K205 also joined in to bind with the phosphate moiety, and this conformation (Fig. S8D) was maintained till the end of the simulation. In the ATP-5, the phosphate moiety first bound with K95 at 0.4 ns and then was relayed between the basic residues in direction of K95→K95K198→K198K205→K198K205K294→K205K294→K205K294K91→K205K294K91K95 in the beginning 12 ns. At 104 ns, K205 left the phosphate moiety, and at 138 ns the adenine base fliped down to stack with Y194 (Fig. S8E). From 992 ns till the end of 1052-ns simulation, the base group moved down to stack with Y190, with the triphosphate group kept anchored to K91K95K294 (Fig.S8F). Therefore, these results show that ATP moves more dynamically between the positively charged residues near the cavity entrance than ADP does (Fig. S2), and ATP in each of these simulations predominantly adopts anti conformation with less favorable binding free energy (Fig. S9). Potential reason could be that ATP has longer triphosphate tail to bind with more Lys residues located on different helices, and therefore is more sensitive to the thermodynamic fluctuations of these helices.

In ATP-6, the phosphate moiety of ATP first bound with K95 at 2.7 ns, and then K91 joined in at 11 ns and R187 joined in at 28 ns. In contrast to the other four simulations, the phosphate moiety of ATP in ATP-6 attached stably with the second basic patch from 28 ns till the end of the 2-μs simulation. Intriguingly, from 146 ns, the adenine base fliped to stack with Y194 or Y190, and ATP bound to the site S2 in a conformation quite similar to the specifically bound ADP (Fig. 4D). Here ATP also adopted the same non-canonical syn conformation which was stabilized by a sodium ion (Fig. 4D, Fig. S9O). However, because of the longer triphosphate tail, when the adenine base of ATP folds back it cannot form specific H-bonds with R187 or N115. Therefore, our simulations demonstrate the sensitivity of site S2 to the length of triphosphate tail of the adenine nucleotides, and ATP can only bind with site S2 in a non-specific way. MM-PBSA free energy calculations show that the ATP binding at site S2 is less favorable than the specific ADP binding by a difference of 6 kcal/mol (Table 1, Fig. S9E), which is consistent with previous experiments on AAC showing lower affinity for ATP [39], [40]. These results imply that only ADP could bind with site S2 in a highly specific mode, and c-state AAC uses this site to screen different nucleotides and specifically recognize its natural substrate ADP.

Increased internal energy may facilitate the specifically bound ADP to leave the second basic patch

During ADP conformational search, the adenine base sometimes reached very deep positions in the cavity (Fig. S10), but it did not guide ADP moving forward to the bottom of the cavity. Instead, the deeply positioned adenine base always flipped up because of the strong attachment of the diphosphate group to the second basic patch. This implies that ADP may inevitably form the specific binding mode before moving forward, and to leave the site S2, the phosphate moiety of ADP needs to break away from the second basic patch first.

To explore potential mechanism for the specifically bound ADP to move forward, we first studied thermodynamic stability of the site S2. As shown in Fig. S11A, the whole region around site S2 is highly stable and is not affected by ADP binding to it. Potential stabilizing factors may include the short C1 loop, the strong D195:R104 salt bridge that attaches C1 loop to H4 (Fig. S11B-C) and the aromatic ladder (F90Y94F98) located at the lipid-facing interface between H2 and H3. This result suggests that ADP breaking away from the second basic patch cannot be induced by the protein local conformational changes around this site. In fact, stability of this region could be very crucial for substrate screening and specific recognition without causing obvious disturbance to the structure of the carrier. Then we calculated time evolutions of the binding free energy and the internal energy of ADP with the MM-PBSA method. The results show that when ADP binds specifically with site S2 in ADP-2 and ADP-3, the internal energy of ADP is increased and keeps increasing at the end of the simulations (Fig. S4L, S5L). The increased internal energy is mainly caused by distorted non-canonical syn conformation (Fig. S4F, S5F) in which the adenine ring is almost at the same plane of the ribose ring (Fig. 4A). The results also show that when the specific binding mode first forms, the binding free energy still fluctuates heavily at this moment (Fig. S4K, S5K). Based on our simulations, dynamic interactions between N115 and R187 and between D92 and K91/K95 may contribute to the fluctuations. In ADP-5, the adenine base is quite flexible, and when ADP adopts syn conformation from 262 ns to 354 ns and from 436 ns to 543 ns (Fig.S6J), its internal energy is also increased (Fig. S6L) and in these conformations a sodium ion also binds simultaneously with O4′ atom, N3 atom and phosphate moiety of ADP (Fig. S3D, E). In the simulations ADP-4 and ADP-6 in which ADP does not adopt syn conformation, internal energy of ADP remains relatively stable (Fig. S6F, R), and the binding free energy of the predominant conformations are less favorable than that of the specifically bound ADP at site S2 (Fig. S6E, Q). Our simulation results suggest that syn conformation is a prerequisite for ADP to form specific binding mode, and vice versa, specific binding with site S2 also helps ADP fix in the syn conformation and further leads to even higher internal energy in ADP. The above results also indicate that the increased internal energy of the specifically bound ADP could facilitate the diphosphate group breaking away from the second basic patch.

An intermediate ADP conformation bound with both site S1 and R187 of site S2

The specifically bound ADP is stably maintained in our simulations although the internal energy of ADP keeps increasing before the end of the simulations. Provided that one ADP/ATP exchange cycle was estimated to be 12–15 ms [22], we expect this specific ADP binding mode would last much longer than current MD simulation time scale. As described in 3.6, to leave the site S2 the phosphate moiety of ADP needs to be released from the second basic patch first. Therefore, to speed up ADP movement in the limited simulation time, we ran a 3.2-μs simulation (ADP-mod) starting from the last snapshot of ADP-3, with partial charges of the phosphate moiety of ADP reduced to 70% of normal values [41] before 2.1 μs and turned back to normal values after that. The simulation revealed a stable ADP conformation that binds with both site S2 and site S1 (Fig. 5A): the adenine base of ADP stacks with R187 of site S2, while the phosphate moiety forms extensive salt bridge network with the six positively charged residues (K22R79R137R234R235R279) of site S1 (Fig. 5B). This conformation formed at 2342 ns and was stably maintained till the end of the simulation. Before forming this stable conformation, the phosphate moiety of ADP first binds with K22K79 at 1311 ns and then with K22R79R137R279 at 1746 ns. During this process ADP adopts the relaxed anti glycosidic conformation, and the adenine base keeps stacking with R187. This stacking interaction pulled the side chain of R187 down to a position quite similar to that in m-state crystal structure (Fig. 5C). This implies that this ADP conformation could be an important intermediate conformation during c-state to m-state transition. MM-PBSA calculation shows that the ADP binding energy in this intermediate conformation (from 2350 ns till the end) is about −97.1 0.6 kcal/mol, much stronger than specific ADP binding mode at site S2. This is caused by more positively charged residues involved in binding with phosphate moiety of ADP. Of special interest, drastic conformational change in the C1 loop near site S2 was observed (Fig. 5D), and we expect this was induced by the strong electrostatic interactions between phosphate moiety of ADP and basic residues of site S1. This is similar to what happens in receptors: ligand binding induces drastic conformational changes at the other side of the transmembrane helices. For AAC, the whole region around site S2 is very stable and not affected by direct ADP binding, but is very sensitive to ADP binding to the distant site S1. We speculate that such thermodynamic feature of site S2 could be crucial for substrate selection and recognition in the highly dynamic transport process. The conformational change of C1 loop during state transition was also demonstrated by cysteine labeling experiment [42].

Fig. 5

An intermediate conformation of ADP binding with the central binding site S1 and R187 of site S2. (A) A stable intermediate ADP conformation bridging site S2 and site S1. Residues of site S2 and site S1 are shown in yellow and violet spheres respectively, and ADP is shown in green sticks. (B) The salt bridge network formed between the phosphate moiety of ADP and the positively charged residues of site S1. (C) Conformations of R187 in different structures. (D) Superposition of the snapshots at the end of traj-2 (silver) and ADP-mod (violet). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Presence of specific site S2 in the upper region of the cavity is supported by sequence analysis among adenine nucleotide transporters

To validate the new specific ADP binding site S2, we compared sequences of transporters that are closely related to AAC. For those closely related transporters, residue variability in the substrate binding area confers specificity for their structurally similar yet different substrates. Therefore, specific substrate binding area could be predicted through identifying variable residues among the closely related paralogs. To this end, we first did phylogeny analysis on protein sequences of all 53 human mitochondrial carriers to identify the closely related paralogs of AAC. As shown in Fig. 6A, AACs fall into the same clade with SLC25A43, SLC25A42, GDC (Graves disease carrier protein, or SLC25A16), SCaMCs (short Ca2+-binding mitochondrial carrier, or Mg-ATP/Pi carrier) and SLC25A41. The substrate of SLC25A43 remains unknown yet, and this carrier is not included in analysis. Both c-state SLC25A42 and GDC specifically recognize and import CoA into the mitochondrial matrix. C-state SCaMCs and SLC24A41 specifically recognize and import Mg-ATP in exchange for intramitochondrial phosphate. ADP, CoA and Mg-ATP are all adenine nucleotides, yet quite different in size of tails and charge distributions (Fig. 6B). Therefore, their binding sites most likely show variability in the charged residues accordingly. To identify the variable charged positions in the five transporters, we first did multiple sequence alignment on the protein sequences from different organisms for each transporter (Fig. S12, S13), and then picked those variable positions that face the central cavity and have conserved charged residue in at least one of the five transporters. As shown in Fig. 6B, totally 12 variable charged positions have been identified. When mapping these variable positions on the structure of c-state AAC, ten of them are located at the upper region of the cavity, while only two (corresponding to D231 and R235 in AAC) are located at the bottom (Fig. 6C). Moreover, D231 and R235 are commonly shared by AACs, SLC25A42 and GDC (Fig. 6B), and hence they cannot be used to discriminate ADP and CoA. These sequence analyses indicate that AAC and other adenine nucleotide transporters use the upper region, rather than the previously proposed bottom region of the cavity, to discriminate and specifically recognize their substrates. This strongly supports our MD simulation results on identifying a highly specific binding site S2 in the upper region of the cavity.

Fig. 6

Adenine nucleotide transporters use the variable positions at the upper region of the cavity to discriminate substrates. (A) AACs fall into the same clade with SLC25A43, SLC25A42, GDC, SCaMCs and SLC25A41 in phylogeny analysis. (B) Sequence logo presentation of the variable positions among adenine nucleotide transporters and their corresponding substrates. Only variable positions that face the central cavity and involve charged residues in at least one of the five transporters are shown. Residues belonging to the specific ADP binding site S2 of AAC are highlighted with red arrows. The residues are numbered based on hAAC1, hSLC25A42, hGDC, hSCaMC1, hSLC25A41 respectively. Complete sequence logos on the odd-numbered and even-numbered helices are provided in Figs. S12, S13 respectively. (C) Mapping the variable positions (shown in spheres) on the c-state structure of AAC.

Among the five residues of site S2 that specifically binds with ADP, four of them (K91, N115, Y190, R187) appear at the variable positions among adenine nucleotide transporters, and these four residues are actually unique for AAC (Fig. 6B). K95 and the other three aromatic residues (Y186, F191 and Y194) of site S2 are commonly shared by adenine nucleotide transporters. Our simulations highlight the central role of R187 in specifically recognizing ADP: its side chain binds with ADP phosphate moiety, its backbone binds with N6 atom in the adenine base of ADP, and meanwhile its side chain stacks with ADP adenine base. Actually, it forms a perfect interlocked structure with ADP, and we expect no other nucleotide could form similar structure with R187. R187 of AAC is replaced with alanine in the other four adenine nucleotide transporters (Fig. 6B). In these transporters, an arginine consistently appears at the same position equivalent to N115 of AAC, and we expect that this arginine also plays a central role in recognizing the adenine group of the substrate Mg-ATP or CoA.

Discussions

In the current work, starting from relaxed apo conformations of AAC and using improved MD simulation strategy for ligand binding, we have identified a new specific ADP binding site S2 in the upper region of the cavity and successfully explained the unusually high substrate specificity of AAC at atomistic level. This new site together with the previously identified central binding site located at the bottom of the cavity supports the early biochemical findings about existence of two specific substrate binding sites [17], [18], [19], [20]. Meanwhile, previous biochemical experiments also characterized the ADP transport process as a two-step event: nucleotide binding to a specific site, followed by the vectorial process of transport [21]. Moreover, conformational transition is linked to the vectorial reaction of transport, rather than substrate recognition [21]. This implies that substrate recognition and conformational transition occur at the two different sites. These previous findings provide a framework to better interpret our simulations results, and the proposed model of ADP binding and transport is illustrated in Fig. 7.

Fig. 7

The proposed mechanism of ADP recognition and transport in AAC. ADP is shown in green, with the adenine base, ribose ring and phosphate moiety represented with rectangle, pentagon and circle respectively. The specific site S2 is shown in yellow and the central site S1 is shown in violet. The first basic patch is shown in cyan. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Our MD simulations suggest that ADP coming from intermembrane space is first attracted to the first basic patch (K198K205) of AAC (Fig. 7A), a conserved non-specific site that was also reported in the closely related SCaMC in previous NMR experiments [43]. Then phosphate moiety of ADP is quickly relayed to the second basic patch (K91K95R187) of the carrier (Fig. 7B). The second basic patch together with the aromatic tyrosine ladder (Y186Y190Y194) including F191 match with the negatively charged phosphate moiety and aromatic base group of nucleotides, therefore this region may show non-specific attractions to a variety of nucleotides. However, only ADP could converge to the highly specific binding mode after rigorous conformational searches (Fig. 7B-C), and therefore this region could function as efficient substrate screening site (site S2). It’s apparent that both electrostatic interactions and stacking interactions are important in ADP recognition process. Significantly, the specifically bound ADP adopts the syn glycosidic conformation, and we speculate that the increased internal energy of ADP could facilitate its diphosphate group breaking away from the second basic patch. Then the diphosphate group of ADP binds to the positively charged residues in the third and fourth basic patches (central binding site S1), with the adenine base still stacking with R187 of site S2. In this intermediate conformation, the strong salt bridge interactions between the diphosphate group and site S1 induce drastic structural changes to the region around site S2 (Fig. 7D), which we expect will further lead to the closure of the cytoplasmic network (Fig. 7E) and then the opening of the matrix network (Fig. 7F).

The most significant result of the current work is the identification of the new specific ADP binding site S2 that is composed of the second basic patch, the tyrosine ladder, F191 and N115 (Fig. 4A), and successful explanation of the unusually high substrate specificity of AAC at atomistic level. The second basic patch and the tyrosine ladder represent the most prominent asymmetric structural elements in the upper region of the cavity in c-state AAC, matching with the highly asymmetric structure of ADP. Our MD simulations demonstrate how AAC utilizes these asymmetric structural elements to specifically recognize ADP and show selectivity over ATP. These results are reliable in that the same highly specific ADP binding mode was consistently observed in two independent long-time MD simulations, and ADP in the two simulations converge to the specific binding mode through the same pathway and via the same highly populated intermediate ADP conformation. Presence of the specific site S2 in the upper region of the cavity is also supported by sequence analyses among adenine nucleotide transporters (Fig. 6). Functional significance of some residues of site S2 has already been validated in previous experiments: it was reported that mutations on R187 or Y190 destroy ATP synthesis (10); double mutations of Y203 and Y207 in yeast AAC (equivalent to Y186 and Y190 in bovine AAC) were reported to affect cellular growth (7); it was also reported that Y203 and F208 in yeast AAC are important for transport activity [45]. Moreover, this specific site S2 of AAC is also supported by previous NMR experiments on SCaMC, in which D360 of SCaMC, corresponding to Y190 of site S2 in AAC, was demonstrated responsible for substrate selectivity [43]. The consistency among the current work and with previous experiments strongly suggest that the highly specific ADP binding site S2 in the upper region of the cavity in c-state AAC has been identified unambiguously.

Of special interest, five (A89D, L97P, D103G, A113P, A122D) of all the nine documented pathological mutations of AAC are located near the site S2 (Fig. S14), and previous experiments showed that these mutations lead to impaired transport of nucleotide [45]. Presence of the disease-causing mutations nearby may also serve as a justification to the newly identified binding site S2, and vice versa, identification of this new site provides potential mechanistic interpretation of these disease-causing mutations. Judging from the structure, we expect that the A122D and A89D mutations might disrupt the site S2 through forming salt bridges with R187 and K91 respectively. The L97P and A113P mutations might introduce helical backbone deformation and respectively affect K95 and N115 of site S2. The D103G mutation might destabilize the R104:D195 salt bridge that is important to stabilize the whole region around site S2.

In the m-state structure, F191 of site S2 together with F88 and L287 form a hydrophobic plug [9] which effectively disrupts the site S2 (Fig. S15). R187 and Y186 are located at the bottom of the cavity and are the only residues of site S2 that are accessible from the matrix side. Positively charged residues of site S1 are distributed in the middle region of the cavity in the m-state. Therefore, it’s possible that when AAC changes to the m-state, the adenine base of ADP may still stack with R187 and phosphate moiety could still partially bind with the positively charged residues constituting S1. That is to say, the orientation of ADP shown in Fig. 5A&B could be maintained during the c-to-m transition, with the adenine base pointing toward the intermembrane space and the phosphate moiety pointing toward the matrix. We predict that before the m-to-c transition, ATP may also adopt similar conformation: the adenine base stacks with R187, and the triphosphate group binds with part of the positively charged residues of site S1.

As mentioned above, early biochemical experiments demonstrate AAC has two specific ADP binding sites of different binding affinities [17], [18], [19], [20]. Our calculations indicate that the newly identified site S2 corresponds to the low-affinity site, while the central binding site S1 is the high-affinity site. Therefore, it’s energetically favorable for ADP to move from site S2 to S1. However, due to the strong electrostatic attraction from the second basic patch of site S2, such movement is not easy kinetically. AACs solve this problem through linking specific binding at site S2 with the high-energy syn conformation of ADP, which could dramatically lower the energy barrier of ADP moving from site S2 to S1. Our results also showed that the internal energy of the specifically bound ADP still keeps increasing at the end of the simulations (Fig. S4L, S5L). Therefore, we expected that provided longer time, the interplay between the specifically bound ADP with site S2 will help ADP cross the energy barrier and move its phosphate moiety from site S2 to S1. The current work together with previous experiments suggest that the highly specific site (S2) in the upper region of the cavity is utilized for substrate selection and recognition, while the less specific central binding site (S1) at the bottom of the cavity is mainly used for triggering conformational transition. Different from receptors or enzymes, substrate binding with transporters is coupled with translocation of substrate, “movement” of binding site and more drastic conformational changes of the protein. Using the secondary binding site for substrate recognition and the central binding site only for conformational transition could be a smart strategy for transporters to cope with substrate recognition problem in the highly dynamic transport process.

Conclusions

Rigorous MD simulations of over 31 microseconds have demonstrated that ADP binding to c-state AAC is a slow multi-step process, and AAC uses the featured tyrosine ladder and the second basic patch to specifically recognize ADP.

Data availability

Data and material are available upon reasonable request to the corresponding authors.

Funding

Research was supported by the National Natural Science Foundation of China (Grant No. 32171241 to X.C.) and the Natural Science Foundation of Zhejiang Province (Grant No. LY18C050002 to X.C.).

CRediT authorship contribution statement

Shihao Yao: Investigation, Formal analysis, Writing – original draft. Qiuzi Yi: Investigation, Formal analysis, Writing – original draft. Boyuan Ma: Data curation, Visualization. Xiaoting Mao: Data curation, Visualization. Ye Chen: Data curation, Visualization. Min-Xin Guan: Writing – review & editing, Supervision. Xiaohui Cang: Conceptualization, Formal analysis, Writing – review & editing, Supervision, Project administration, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.