A conserved salt bridge network stabilizes the hepatic organic anion transporters OATP1B1 and OATP1B3

TL;DR

This study identifies a salt bridge network that stabilizes liver transporters OATP1B1 and OATP1B3, influencing their function and conformational changes.

Contribution

The discovery of a conserved salt bridge network centered around E185 that stabilizes the structure and function of OATP1B1 and OATP1B3.

Findings

A salt bridge network centered around E185 is crucial for the uptake activities of OATP1B1 and OATP1B3.

The salt bridge network stabilizes the inward cavity and connects the N- and C-bundles of the proteins.

The network changes with conformation and does not involve ligand coordination.

Abstract

The organic anion transporting polypeptide (OATP)-1B1 and -1B3 are liver-specific transporters that govern the uptake of numerous endogenous molecules and drugs before their metabolism and excretion by the hepatocytes. Structurally, these two transporters are members of the major facilitator superfamily, operating by the alternating access mechanism that facilitates the movement of solutes between extracellular and intracellular compartments. Given their dynamic nature, salt bridges often modulate the conformations of transporters and participate in the orchestration of conformational changes. In this study, we identified and characterized a network of salt bridges within the internal cavities of OATP1B1 and OATP1B3 by cell-based uptake assays, uptake kinetics, and molecular dynamics simulations. These experiments revealed that a salt bridge network centered around E185 is crucial for uptake activities in these two proteins, as it stabilizes the inward cavity of the proteins and bridges the N- and C- bundles of the protein. Interestingly, this salt bridge network changes as a function of conformation. Furthermore, the residues studied do not participate in ligand coordination in the published structures nor in our simulations. These findings advance our understanding of the elaborate network of ionic interactions that govern the structure and dynamics of OATP1B1, OATP1B3, and other MFS transporters.

Figures

Full text

The organic anion transporting polypeptide (OATP) family of plasma membrane transport proteins is responsible for the multispecific uptake of numerous endogenous molecules, food constituents, drugs, and toxins. To date, more than 300 OATP members have been identified in both vertebrates and invertebrates, classified as the SLCO family in the solute carrier classification (1). The 11 human OATPs are divided into families, subfamilies, and individual genes according to phylogenetic relationships and amino-acid sequence identities (2). Some members of the OATP family display wide tissue distributions and/or ligand multispecificities, while others are organ-specific and/or display high levels of substrate specificity. Many of the identified substrates of the OATPs are structurally divergent, ranging in charge from anionic, cationic, neutral, or zwitterionic with diverse chemical attributes and sizes (3). Due to their multispecificities and wide distribution in tissues that govern drug absorption, distribution, metabolism, and excretion such as the intestine, liver, kidney, and brain; the OATPs influence the pharmacokinetics of several highly prescribed drug classes such as statins and antivirals (4) and are involved in several drug-drug interactions that elevate systemic substrate concentrations, at times causing life-threatening conditions (5, 6). Furthermore, OATPs are abnormally expressed in a number of cancers, such as colon (7), prostate (8), pancreatic adenocarcinomas (9), breast cancer (10, 11), and in multiple cancer cell lines (reviewed in (12)). Moreover, sulfated hormones that promote the proliferation of hormone-dependent breast (13, 14) and prostate cancer (15) are thought to enter these malignancies through overexpressed OATPs both in vitro and in vivo. Finally, multiple chemotherapeutics such as methotrexate, doxorubicin, paclitaxel, docetaxel, irinotecan, and some tyrosine kinase inhibitors, are transported by various OATPs into cancerous lesions as well as into hepatocytes, affecting both their therapeutic effects and clearance rates (16).

Despite the functional importance of the OATPs, there were no published structures of any member of this family until late 2023 when the Locher and Zhang groups published two articles a few days apart (17, 18). This is yet another testament to the difficulties in obtaining high-resolution structures of mammalian membrane proteins. Prior to these publications, homology models (1, 19, 20, 21, 22) were used to generate structure-based hypotheses. Structurally, the OATPs belong to the major facilitator superfamily (MFS) fold with 12 transmembrane (TM) helices that are organized into two bundles of 6TM helices, the N- and C-terminal bundles (Fig. 1).Figure 1OATP1B structures and internal symmetry. A, an OF structure of OATP1B1 (PDB ID 8k6l, left) and an IF structure of OATP1B1 (PDB ID 8hnd, right) with the N- and C-bundles colored yellow and cyan, respectively. B, a superposition of the N- and C-bundles of OF OATP1B1 with an RMSD of 3.7 Å. C, superposition of the transmembrane (TM) domains of seven OATP1B structures. All-against-all matrix of RMSD values for the superposition of the TM domains is shown. The PDB ID codes are provided for each structure, and “x” values along the diagonal indicate comparison of the structure against itself. Heatmap colors go from green-yellow (low RMSD values below 1.2 Å) to red (high values above 4 Å) for panels C–E. D and E, superposition of only the N-bundle domains (panel**D) or C-bundle domains (panel**E). Table 2 lists the PDB IDs of the solved structures, and their conformations. IF, inward-facing; OATP, organic anion transporting polypeptide; OF, outward-facing; PDB, Protein Data Bank.

During the transport cycle, MFS proteins undergo rocker-switch motions in which the two bundles rearrange around a centrally located ligand binding site that is alternatively exposed to the extracellular space (outward-facing (OF)) and intracellular space (inward-facing (IF)) (left and right structures in Fig. 1A, respectively). The two bundles are related by a pseudo two-fold symmetry axis that runs through the center of the transporter and perpendicular to the plane of the membrane. Superposition of the TM domain of one bundle on the other yields a similar helical arrangement with RMSD values around 5 Å for both OATP1B1 and OATP1B3 (Fig. 1B). The large-scale isomerization rocker-switch events are often coupled with additional asymmetric local, gating rearrangements within each bundle (intrabundle) and between the two bundles (interbundle) that are often controlled by salt bridges (reviewed in (23, 24)). There are now 10 deposited structures of OATP1B1 and OATP1B3 in various states with and without ligands and small molecules bound (17, 18, 25), but they generally cluster broadly into these two different IF and OF isomerization states, as exhibited by the all-against-all RMSD value matrix for the TM regions shown in Figure 1C. The first seven OF structures in the matrix cluster together with RMSD values less than 1, while the second group of IF conformations share low RMSD values of 1 to 1.2 Å. However, comparing one group to the other provides large RMSD values above 3.8 Å (red rectangle in Fig. 1C). Interestingly, OATP1B3 has only been solved in the IF state (Protein Data Bank (PDB) ID 8pg0), yet its RMSDs to IF OATP1B1 structures are very low, 1.1 to 1.2 Å. Finally, the individual bundles adopt nearly identical conformations within the IF and OF states, as revealed by bundle only superpositions in panels D (N-bundles) and E (C-bundles). For instance, the N-bundle from the IF OATP1B3 structure superposes onto the N-bundle of all OF OATP1B1s with RMSD values 0.7 to 1.0 Å which is just as good as the IF OATP1B1 bundles (bottom row of matrix in Fig. 1D).

Salt bridges form important interactions, that are often conformation dependent, to control the integrity and dynamics of transporters. OATP1B1 and OATP1B3 have a number of salt bridges in the cavity. In this study, we characterized a conserved complex salt bridge network in the N-bundle involving E185 and R181 from TM4 and K41 from TM1 by molecular dynamics (MD) simulations in the presence and absence of transported molecules, ligand uptake measurements, and site-directed mutagenesis. The intrabundle K41-E185 salt bridge appears to govern the stability of the N-bundle, dictates discrete helical movements during the transport cycle, and controls the intracellular gate, while according to our MD simulations, E185 also forms interactions with R580 from the C-bundle only in the OF conformation, coupling the two bundles together. Our findings highlight the complexity and dynamics of salt bridge interactions and possible involvement of interbundle and intrabundle salt bridges as molecular determinants of transport through members of the OATP family and likely other MFS transporters.

Results

The functional importance of charged residues in the OATP1B family has been the focus of various studies (26, 27, 28, 29) due to the propensity of charged residues to form interactions with ligands, or repel molecular entities of similar charges. Furthermore, charged residues could potentially take part in salt bridges that control the stability and interactions of helices in the transporters, gating, and isomerization between different conformations. As an example, in the lactose transporter LacY, E126, and R144 form an intrabundle salt bridge in the OF conformation that breaks upon lactose binding to these two residues, supporting the transition of LacY to an IF conformation (30). In OATP1B1 and 1B3, the conserved E185 is located in the middle of TM4. On the basis of sequence conservation and structural models that were later verified by cryo-EM structures, we hypothesized that this charged residue could be involved in a salt bridge network with close-by positively charged residues to stabilize the protein. As such, we investigated the possible involvement of E185 with K41 and R181 in the formation of a complex salt bridge network, the influence of mutations of these residues on uptake levels, kinetics, and protein dynamics.

Conservation of residues in TM1 and TM4 in the OATP1 family

According to the OATP superfamily nomenclature, proteins with amino acid sequence identities of ≥40% belong to the same family (that is denoted by a numeral, i.e. OATP1) while proteins with a sequence identity of ≥60% belong to the same subfamily (that is denoted by a capital letter, i.e. OATP1B) (1). A ClustalOmega sequence alignment (31) of the eleven human OATPs reveals a conserved glutamate (E185) and arginine residue (R181) in TM4 and a conserved lysine (K41) in TM1 in all four OATP1 members (Fig. 2A). Table S1 lists the exact primary sequence position of these three residues in each OATP1 member. A more comprehensive alignment of the ∼180 OATP1B members curated in Uniprot shows that E185 (1B family numbering) and R181 are conserved in ∼98% of the cases, whereas a positively charged residue is preferred for position 41 with K at 72% of the sequences and R at 28% according to Jalview analysis (32) (Fig. 2B and Supplementary File 1).Figure 2Sequence alignments of human OATPs and representative OATP1B subfamily members. A, sequence alignment of the eleven human OATPs around the residues studied in this article (marked in hexagons above the alignment). K41, R181, and E185 (OATP1B1 and OATP1B3 numbering) are conserved in the OATP1B subfamily. B, sequence alignment of a representative sample consisting of human OATP1B1 and the top 10 similar proteins in the Uniprot database show a preference for positively charged residues, either K (72%) or R (28%), at position 41, while positions 181 and 185 display a strict selection for R181 and E185. The logo representation at the bottom of this panel corresponds to the complete sequence alignment of ∼180 members of the OATP1B family. Supplementary file 1 contains the full alignments of 176 proteins that share a >60% identity with OATP1B1 and OATP1B3. Sequences with >50% identity to OATP1B1 and OATP1B3 were obtained from Uniprot and trimmed to the 60% threshold. Alignment was generated by ClustalOmega (31) and displayed using the Clustal color scheme. The logo representation generated by JalView (32) is displayed at the bottom, showing the consensus residues at each position. OATP, organic anion transporting polypeptide.

Salt bridges

Our work started prior to the release of the cryo-EM structures of the OATPs. A published OATP1B1 model (21), as well as models later released by AlphaFold2 (AF2) (22) and other studies (19), suggested that E185 in OATP1B1/1B3 faces a putative cavity at the center of the protein. Two positively charged and conserved residues are in proximity to E185 and potentially form salt bridges with it (20). With the publication of the OATP1B1/3 cryo-EM structures, Shan et al. explicitly described the K41-E185 salt bridge in the ligand binding cavity, as well as K49-D70 and R580-E74 (Fig. 3) (17). We then compared the E185 and K41 side chain distances in the accompanying structures. Some support the formation of a salt bridge, such as the IF conformation of OATP1B1 with estrone 3-sulfate (E1S) bound in which the K41-E185 side chain distance is 2.7 Å (17), while other structures in the same publication, and the two structures published by Ciuta et al., have distances up to 5.8 Å (18) (Table 1). These different structures adopt both inward and outward-facing conformations, and interestingly, two of the IF conformations, both bound to the same E1S substrate, provide different K41-E185 distances in OATP1B1 with the Ciuta et al., structure being quite far at 5.8 Å, and the structure from Shan et al., being much shorter at 2.7 Å (17). These differences arise, in part, because of different side chain placements in each structure, but E185 is found at a distorted part of TM4, suggesting that TM4 might be ideally positioned to dictate conformational changes. Furthermore, the map density of many of the side chains is quite poor, which lowers the confidence in the exact placements of these side chains in the cryo-EM structures. Table 1 lists the K41-E185 side chain and Cα distances in the published OATP1B structures.Figure 3Salt bridges in the binding cavity. A, the IF DCF-bound structure of OATP1B1 is shown in the membrane plane. TM5 was removed to better show the area involved in the formation of the salt bridges studied. B, a zoomed-in view of the area enclosed by a box in a. TM numbering is shown as white numbers on top of helices. C, a tilted view from the cytoplasm of the exposed cavity. In all panels the N-and C-bundles are plotted as yellow and cyan cartoons and salt bridges identified in the structure are show as black lines. DCF, 2′,7′-dichlorofluorescein; IF, inward-facing; OATP, organic anion transporting polypeptide; TM, transmembrane.Table 1K41-E185 side chain and Cα distances in published OATP1B structuresProteinPDB IDConformationBound moleculeK41-E185 side chain distance [Å]K41-E185 Cα distance [Å]OATP1B18K6LOFDCF3.6, 3.77.3OATP1B18HNBOFApo3.66.3OATP1B18HNCOFBilirubin2.96.5OATP1B18HNHOFSimeprevir2.86.9OATP1B18HNDIFE1S2.77.1OATP1B18PHWIFE1S5.87.2OATP1B38PG0IFBicarbonate5.56.6OATP1B19CY1OFApo4.56.3OATP1B19CY3OFAtorvastatin3.56.5OATP1B19CY4OFCyclosporin3.56.6

In order to differentiate if E185 exclusively forms a salt bridge with K41 or also with R181 and probe the influence of this salt bridge on the activity of OATP1B1 and 1B3, we employed uptake and kinetic studies in HEK293T cells overexpressing wild type, E185, K41, and R181 mutants of either proteins. First, we mutated the three residues to alanines and probed the uptake of the fluorescent ligands, 2′,7′-dichlorofluorescein (DCF) for OATP1B1 and 4′,5′-dibromofluorescein (DBF) for OATP1B3, into adherent HEK293 cells in order to probe the contribution of each side chain to the transporter activity. The membrane abundance of all OATP1B1 mutants studied ranged between 73 to 110% relative to the expression levels of WT protein and 87 to 105% for OATP1B3 (Table S2 for OATP1B1 and Table S3 for OATP1B3). E185A and K41A mutations of OATP1B1 and OATP1B3 resulted in a complete loss of activity whereas R181A retained 65% activity for OATP1B1 and 73% for OATP1B3 (Fig. 4). All uptake data displayed are normalized to plasma membrane protein levels and representative blots used for quantification are shown in Figure 4B.Figure 4Ligand uptake levels of WT and alanine mutants of E185, K41, and R181. A, DCF uptake levels for OATP1B1 (white-filled bars) and DBF for OATP1B3 (gray-filled bars) were monitored. The alanine mutations resulted in complete loss of activities for E185A and K41A while the R181A mutation retained 65% and 74% of activities for OATP1B1 and OATP1B3. Results are calculated as mean ± standard deviation of three separate experiments done in triplicate with data displayed as fractional activity relative to that of the WT protein. All data are normalized to cell surface expression levels. Significance was determined through ANOVA analysis with ∗∗∗∗ indicating a p < 0.0001 and ∗∗∗ indicating p = 0.0003 relative to WT proteins. B, representative western blots used for cell surface quantification. Top: anti-FLAG blot of biotinylated cell surface proteins from HEK293T cells over expressing various OATP1B1 constructs. Each lane is for an independent experiment. Bottom: anti-CD71 blot of same lysates as a loading control. The lower prominent bands are for monomeric CD71 while the top most band is for dimeric CD71. The combined intensities of both bands were used for the normalization of OATP expression levels. DBF, 4′,5′-dibromofluorescein; DCF, 2′,7′-dichlorofluorescein; OATP, organic anion transporting polypeptide; WT, wild-type.

The loss of activity for E185A and K41A prompted us to test uptake with more conservative mutations at these positions that preserve the charge character of the original residue (Fig. 5). Maintaining a positive charge for K41 with the K41R mutation that is also common in the OATP1B family (Fig. 2B) retains transport, but uptake levels are significantly suppressed at 19% of WT for OATP1B1 and 20% for OATP1B3. The E185D mutation led to complete loss of activity, as well as the E185Q mutation, suggesting that the negative charge and space occupied by the native glutamic acid side chain are important for transport. The single point mutants K41A, E185A, and E185Q all disrupt the salt bridges while also leaving an uncompensated charge in the area. This local charge could potentially disrupt the activity of the protein and explain the loss of activity. To distinguish if the loss of activity is due to the disruption of the salt bridge or unbalanced charge, we generated two double mutants that do not alter the local net charge: K41A-E185Q and K41Q-E185Q. Both of these mutants led to a complete loss of activity, supporting the notion that the loss of activity is not due to unpaired charges. Finally, we wanted to test if swapping the charge pair with K41E-E185K that keeps the paired charges in this site albeit on the opposite residues would impact activity. This swap also results in a complete loss of activity, similar to results in a swap mutant of LacY (33). We hypothesize that the loss of activity for the swap mutant may be related to the need for E185 to salt bridge to other basic residues, as expanded upon later.Figure 5Ligand uptake levels of WT, conservative, charge neutralization, and swap mutants for K41 and E185. DCF uptake levels for OATP1B1 (white-filled bars) and DBF for OATP1B3 (gray-filled bars) were measured. K41R mutations exhibited 19% of WT uptake levels for OATP1B1 and 20% for OATP1B3. Conserving a negative charge with E185D did not rescue transport levels. A single charge neutralization mutant, E185Q, and double charge neutralizing mutants, K41A E185Q and K41Q E185Q also resulted in loss of activity. A swap mutant K41E E185K also did not exhibit any uptake. Results are calculated as mean ± standard deviation of three separate experiments done in triplicate with data displayed as fractional activity relative to the activity of the WT protein. All data are normalized to cell surface expression levels. Significance was determined through ANOVA analysis with p < 0.0001 levels for all samples relative to WT proteins. DBF, 4′,5′-dibromofluorescein; DCF, 2′,7′-dichlorofluorescein; OATP, organic anion transporting polypeptide.

Further, we performed a kinetic analysis of ligand uptake for the conditions that resulted in detectable uptake levels. Namely, WT, R181A, and K41R of OATP1B1 (Fig. 6A), and WT and R181A of OATP1B3 (Fig. 6B). DCF uptake by WT and R181A OATP1B1 are described by Michaelis–Menten curves with almost identical Michaelis constants (Km) of 4.83 μM (WT) and 4.41 μM (R181A). But, R181A exhibited a 40% reduction in V_max_ as compared to WT protein (Fig. 6A. and Table 2), in agreement with the reduction seen in the measured uptake levels (Fig. 4). The 1.90 μM km of K41R was slightly lower than WT OATP1B1, but with a V_max_ reduced by 36.8% as compared to WT. The low levels of uptake for the E185 mutants in both proteins and K41 in OATP1B3 precluded us from performing a thorough kinetic analysis. The kinetic characterization of OATP1B3 for the WT protein and R181A mutant mirror the results of OATP1B1 - the Km values of WT and R181A toward DBF are 2.06 and 1.00 μM, respectively, with the mutant exhibiting a 26.4% decrease in V_max_.Figure 6Uptake kinetics of WT, R181A, and K41R. A, uptake kinetics of WT (○), R181A (∇), and K41R (◻︎) OATP1B1. B, uptake kinetics of WT (○), and R181A (∇) OATP1B3. Data are calculated as mean ± standard deviation of two separate experiments each done in triplicate and is displayed as fractional activity relative to WT protein. Data are normalized to membrane expression levels with relative amounts of all OATPs determined by cell surface expression assay and relative abundance averaged over three independent experiments. OATP, organic anion transporting polypeptide.Table 2. Kinetic parameters for OATP1B1 WT, R181A, K41R and OATP1B3 WT, R181AConstructKm [μM]Km CI [μM]Vmax [%]Vmax CI [%]Relative abundance1B1-wild type4.833.78–6.17100107.2–124.811B1-R181A4.413.57–5.4560.657.0–64.70.99 ± 0.031B1-K41R1.901.43–2.5036.834.4–39.40.94 ± 0.021B3-wild type2.061.72–2.46100105.8–115.111B3-R181A1.000.74–1.3673.671.5–80.20.96 ± 0.05

Taken together, our uptake and kinetic data strongly suggest that the exact chemical makeup and spatial arrangement of the side chains of E185 and K41 are crucial for transport via OATP1B1 and OATP1B3, while R181 is not critical for ligand uptake. Wang et al. (28) previously studied TM1 of OATP1B1 and obtained similar results for K41. Namely, the uptake levels of taurocholate and E1S were significantly reduced as compared to the WT protein, and K41R displayed higher uptake then K41A.

Molecular simulations

To probe the roles of the charged residues K41, R181, and E185 on the stability of the proteins, their dynamics, and ability to bind substrate we performed a series of MD simulations on three different systems (IF OATP1B1, OF OATP1B1, and IF OATP1B3). Simulations were carried out under four different conditions: (1) apo, (2) E185A apo, (3) DCF bound and, (4) DBF bound, although not all conditions were simulated for each protein system (see Table 3).Table 3. List of molecular dynamics simulationsStructureConditionSimulation timeOF OATP1B1 (8k6l)apo3 × 500 nsOF OATP1B1 (8k6l)E185 A3 × 500 nsOF OATP1B1 (8k6l)+DCF3 × 500 nsOF OATP1B1 (8k6l)+DBF3 × 500 nsIF OATP1B1 (8hnd)+DCF3 × 500 nsIF OATP1B3 (8pg0)apo5 × 500 nsIF OATP1B3 (8pg0)E185 A5 × 500 nsIF OATP1B3 (8pg0)+DCF6 x ∼1–400a nsIF OATP1B3 (8pg0)+DBF3 x ∼2–300a nsaSome simulations were halted before the “full” 500 ns run after the ligand was observed to have unbound.

First, we analyzed the Cα distances between residues 41 to 185 and 181 to 185 in the WT and the E185A mutants of the OF conformation of OATP1B1 and the IF conformation of OATP1B3 (Fig. S1). In all cases, there were minimal distance fluctuations in good agreement with the distances derived from the cryo-EM structures, and no differences between the WT and E185. This suggests that at the secondary structure level, the E185A mutation does not change the local secondary structure of TM4 and the local tertiary structure relationship between TM1 and TM4 on the 500 ns timescale. Then, using the angle between the N- and C-bundles as a metric for the spatial relationship between the two bundles (Fig. 7A), a striking feature emerged. Several of the IF OATP1B3 WT simulations open their inner gate even farther than in the cryo-EM structure (Fig. 7B), while the inner gate begins to close in several simulations when the E185A salt bridge breaking mutation is introduced (Fig. 7C). Five repeat simulations were run in each case, and all of the WT simulations either keep the opening angle constant (ice blue curve in panel B) or increased it by over 10° compared to the initial angle prior to equilibration (dashed horizontal line). In the case of the E185A mutant runs, the original glutamate was mutated in silico to alanine, but otherwise the system was identical. Three of the replicates showed less drastic opening on the 500 ns timescale with a magnitude of 5° instead of 10°, while two of the replicates even closed the gate reducing the angle in one case by almost 10° (yellow curve in panel C). Next, we visualized the inner gate from snapshots of each OATP1B3 replicate at the end of the simulations, viewed from the cytoplasm, to confirm the tendency of the E185A structure to close the intracellular cavity compared to the WT protein (Fig. 7, D and E). As suggested by the bundle angles, all WT simulations show that a significant separation distance between the N-bundle (dark gray) and C-bundle (light gray) lobes of the protein was maintained. Specifically, this can be seen by noting that the two lobes are not in physical contact, except for the upper edge of 1 snapshot (second from right in panel D), and there is a significant empty white space between the lobes in the center of the transporter (water has been removed for clarity and the visualization plane clipped so the contact on the extracellular side of the membrane is absent). This empty white space corresponds to an open inner cavity ready to exchange ligands with the cytoplasm. However, the E185A snapshots reveal a different picture. The first and third images show completely closed inner gates, and the two lobes are in physical contact. The fourth and fifth snapshots remain open, but the central cavity is more narrow than most WT images in panel D. Meanwhile, the second snapshot in panel E is the most open of the E185A simulations. So there is a close correspondence between the bundle angle and the inner gate closure, and the mutant transporter shows a tendency for closing while the WT does not. We also realized that inner leaflet lipids (yellow in panels D and E) gain access to the vestibule when the gate remains open. When the polar headgroups of the lipids penetrate, the gate cannot close due to sterics, but the cavity remains well solvated and open for exchange with the cytoplasm.Figure 7E185A induces closure of the inner gate, but not of the outer gate. A, cartoon of the N- and C-bundles with vectors defined within each bundle that are used to calculate the angle between the domains. This angle is used to examine the degree of gate closure. B and C, five independent replicates of the IF OATP1B3 WT structure (PDB ID 8pg0, panel**B) and E185A mutant (panel**C) plotting the bundle angle in time. The vertical dashed lines are the initial crossing angle prior to equilibration. D and E, final snapshot of IF OATP1B3 WT (panel**D) and E185A mutant (panel**E) trajectories from panels b and c, respectively, viewed from the cytoplasm with the N-bundles in dark gray, the C-bundles light gray, and the lipids yellow. Waters are not shown, and the scene is clipped in the middle of the membrane plane so that the inward gates and central cavity (white space) can easily be viewed. F, cartoon from panel A, with identical definitions, but overlayed on an OF state. G and H, three independent replicates of the OF OATP1B1 WT structure (PDB ID 8k6l, panel**G) and E185A mutant (panel**H) plotting the bundle angle in time. The vertical dashed lines are the initial crossing angle prior to equilibration. Solid traces in B, C, G, and H are smoothed with respect to the transparent data, gray instantaneous values. IF, inward-facing; OATP, organic anion transporting polypeptide; OF, outward-facing; PDB, Protein Data Bank.

We also performed WT and charge neutralizing simulations of E185 on the OF OATP1B1 structure (8k6l) after removing the bound DCF molecule. Using the same metric to determine the relative angle between the N- and C-bundles as the previous simulations, the crossing angle now reports negative values indicating an outward-open structure with the binding site facing the extracellular space (Fig. 7F). In this case, the outward-open structures start at −26° and −27°, respectively (dashed horizontal lines in panels G and H), and the gates open farther as the angles become more negative. Each of the three replicates was stable with a slight tendency to open further by 2° to 6° (Fig. 7, G and H). Hence, our simulations suggest that the E185A charge neutralization has its greatest impact by destabilizing the IF state, while having little-to-no impact on the stability of the OF state.

To determine if the motions described in Figure 7 represent an unbiased interpretation of the protein conformational changes, we performed principal component (PC) analysis on the aggregated MD data (Fig. S2, Movies S1,2). The dominant mode across all systems is indeed a rocker-switch motion (describing 77.3% of the variance), and further inspection confirms our earlier interpretations: (1) the WT IF state opens wider than the OF state (Fig. S2A), and (2) the IF E185A mutant’s movement along this direction is damped, exhibiting a smaller angular change (32° versus 39°). Most importantly, this PCA confirms that E185A mutation causes a larger shift in the conformational landscape in the IF than in the OF state.

To identify the specific helices that shift during the simulations, we performed an RMSD analysis for individual TM helices (Fig. 8). For WT OATP1B3 in the IF conformation, nine of the helices are quite static, while TM2, 5, and 11 display fluctuations in some or all replicates. Comparing the starting model (in gray) to the final snapshot in one of the trajectories that displayed pronounced helical movements (in yellow) shows that TM4, and TM5 are mainly responsible for the small opening of the vestibule (Fig. 8A). For the E185A mutant, seven of the TMs are static, while TM2, 4, 5, 8, and 11 display higher levels of variability in their positions throughout the simulation time. A structural comparison of the starting model (in gray) and final snapshot (in green) for one of the most dynamic simulations shows the closing of the intracellular cavity mainly by a constriction of the ends of TM4 and 11.Figure 8Opening and closing of the intracellular cavity ofIFOATP1B3 in WT****(A), and E185A (B). Per-helix RMSD plots for five replicate MD simulation trajectories relative to the starting structure (PDB ID 8pg0 with modeled loops). Snapshots in the right hand sides show the final frames of trajectories with the most pronounced helix movements (yellow trace for 1B3-IF-WT, light green trace for 1B3-IF-E185A), aligned to the starting structure (in gray). Note how TM4 and TM5 move out of the cavity in WT, and inward in E185A. Transparent lines are raw data, solid lines show a 10 ns running average. IF, inward-facing; MD, molecular dynamics; OATP, organic anion transporting polypeptide; PDB, Protein Data Bank; TM, transmembrane.

We also analyzed the salt bridge interactions between K41-E185 and R181-E185 as shown in Figure 9 for the apo WT OF OATP1B1 and IF OATP1B3 simulations. The dotted red lines are the initial values in the cryo-EM structures, while the dotted black lines at 2.6 Å in each panel are consistent with a very strong salt bridge (generally the nitrogen to oxygen distance between residues participating in a salt bridge should be under 4.0 Å (34)). The plots show that in simulations of both states a tight salt bridge interaction between K41 and E185 is maintained, despite that these two residues do not form a salt bridge in the IF structure (panel C). The fluctuations across the three independent trajectories, however, are lower in the OF OATP1B1 state, suggesting that this interaction is stronger in the OF state (Fig. 9A). Meanwhile, the R181-E185 salt bridge is only maintained in the IF OATP1B3 simulation (Fig. 9D), and is broken in the OF OATP1B1 state (Fig. 9B). This is consistent with the starting structures showing R181 to E185 N-to-O distances of 3.46 Å for IF OATP1B3 and 6 Å for OF OATP1B1. Thus, this analysis suggests that the IF state the transporter forms salt-bridging interactions between K41-E185 and R181-E185 simultaneously, but in the OF state only the K41-E185 interaction exists, and it becomes stronger since there is no longer competition between K41 and R181 for E185.Figure 9Sa****lt bridge formation of WT simulations of the IF and OF states. A and B, OATP1B1 WT simulations of the OF state (PDB ID 8k6l) plotting the side chain distance between the Nζ atom for K41 and the closer of the Oε1 or Oε2 distal oxygen for E185 (panel**A) and a η nitrogen of R181 and an Oε atom of E185 (panel**B). C and D, OATP1B3 WT simulations of the IF state (PDB ID 8pg0) plotting the side chain distance between the Nζ atom of K41 and the closer of the Oε1 or Oε2 distal oxygen for E185 (panel**C) and either Nη1, or Nη2 of R181 and Oε1 or Oε2 of E185 (panel**D). The dashed red lines are the initial distances from the cryo-EM structures, and the dashed black lines are 2.6 Å reference values. IF, inward-facing; OATP, organic anion transporting polypeptide; OF, outward-facing; PDB, Protein Data Bank.

To provide a better understanding for the loss of activity caused by the K41E-E185K swap mutant, we analyzed the local environment of the two residues in our simulations and noticed that the C-bundle R580 is in close proximity to E185 only in the OF conformation, but it forms a salt bridge with E74 on the N-bundle in both conformations. The E74-R580 salt bridge was also described in the Shan et al., structural paper (17). This bond is closer and better positioned in the OF state due to a large scale movement of TM2 (in which E74 resides) between the OF and IF states. In the OF state a kink develops right above R580 to open the outer gate, and it appears as if E74-R580 form a “pin” to aid in holding the N and C lobes together. The simulations show the E74-R580 interaction is primarily mediated via the ε nitrogen and one of the distal η′s leaving the other η nitrogen free to also salt bridge to E185 (Fig. S3, A and B). This interaction is possible because E185 no longer contacts R181 in the OF state and the N and C lobes have moved closer than 13 Å Cα-Cα. Meanwhile, in the IF state E185 is sequestered in the N-lobe by K41 and R181, thus R580 breaks contact with E185, and the N-to-C lobe separation increases to 16.6 Å compatible with inner gate opening (Fig. S3A). To quantify the conformation-specific salt bridges mentioned above we generated fingerprint-type plots that describe the range of distances between E185-K41, E185-R181, and E185-R580. We also plotted the Cα distances for E185-R580 to highlight the separation between the two residues in the IF conformation and the subsequent closeness in the OF conformation. The plots show that, indeed, in the IF state, E185 always interacts with K41 and R181, and never with R580. This lack of interaction with R580 is consistent with a large (>12.5 Å) Cα-Cα distance between E185 and R580. In contrast, E185 switches back and forth between R181 and R580 in the OF state. The latter interactions corresponds with shorter Cα-Cα distances between E185 and R580 in the OF conformation (Fig. S3 bottom). The notion that the IF state was more sensitive to the E185A mutation, while lacking the E185-R580 salt bridge suggests that the loss in transport activity from this mutant is not likely to be caused by interbundle coupling between the N- and C- lobes, but rather through the intrabundle salt bridge network involving K41 and R181. This agrees with previously published mutagenesis data showing that R580A mutation does maintain some transport activity (35).

Finally, we performed simulations on the WT structures in the presence of the transported substrates DCF and DBF to probe the determinants of substrate binding. The only available structure that resolved the binding pose of DCF is the OF OATP1B1 structure (8k6l), so we superposed the helices that comprise the binding pocket on to the other two structures and moved DCF into those systems. This did not result in any clashes between the protein and DCF, but the rotamers of several sidechains close to the bound DCF were adjusted to better match that seen in the 8k6l structure (see Experimental procedures). For the solved OF OATP1B1 (8k6l) structure, three replicates were run for 500 ns each, and the molecule remained stable with an RMSD to the starting pose of ∼1 Å, except for a partial unbinding event in replicate 2 at ∼450 ns resulting in a rise in the RMSD to 4 Å that subsequently rebound just before 500 ns (Fig. 10A). The starting placement of DCF determined in the OF OAT1B1 structure is pictured in pink in the middle and right panels, and the yellow box in the center panel delineates the zoomed in region to its right showing the final snapshot of replicate 2 at 500 ns. This snapshot has a 1.2 Å RMSD typical of all three traces throughout the trajectories, and the deviation from the starting pose (black molecule) is minor. Meanwhile, DCF remained stably bound in all three IF OATP1B1 simulations even though the molecule was hand placed in this conformation as it was not initiated from a solved structure (Fig. 10B). This stability arises because the binding site remains intact during the OF-to-IF transition and DCF’s interaction fingerprint with the protein is largely unchanged, explaining why the molecule remains bound. As in the OF trajectories, the small deviations are characteristic of thermal motions in a stable binding site.Figure 10Stability of DCF and DBF in OATP1B1 and OATP1B3. A, time series traces of the RMSD of DCF to the initial DCF position in the experimentally determined OATP1B1-DCF co-structure (PDB ID 8k6l) in the OF conformation. Three independent replicates were run. Middle image, initial starting model with DCF in pink. Yellow box indicates the zoomed in region pictured to the right. Zoomed in figure shows the final DCF snapshot (pink) for replicate two compared to the starting DCF position (black). Key charged residues K41, R181, and E185 are drawn as licorice in all panels. B, time series data from three independent simulations, as in panel**A, for DCF placed in the IF OATP1B1 structure (PDB ID 8hnd). Middle image, initial starting model with DCF in pink. Yellow box indicates the zoomed in region pictured to the left. Zoomed in figure shows the final DCF snapshot (pink) for replicate 2 compared to the starting DCF position (black). (C) DCF RMSD over time from six independent simulations of DCF placed in the IF OATP1B3 structure (PDB ID 8pg0). Cartoon models show the final DCF snapshot (pink) for each replicate compared to the starting DCF position (black). D, DBF RMSD over time from three independent simulations of DBF placed in the IF OATP1B3 structure (PDB ID 8pg0). Cartoon models show the final DBF snapshot (pink) for each replicate compared to the starting DCF position (black). For panels C and D, the zoomed in area is similar to that shown for the IF OATP1B1 structure in panel**B. In addition, the three charged residues K41, R181, and E185 adopt the same positions as those labeled in panel**B. DBF, 4′,5′-dibromofluorescein; DCF, 2′,7′-dichlorofluorescein; IF, inward-facing; OF, outward-facing; PDB, Protein Data Bank.

Next, as a control, we initiated multiple simulations of DCF in the IF OATP1B3 structure with the small molecule placed in the pocket as in panel B. DCF has a much lower affinity to OATP1B3 then to OATP1B1 and unlike the IF OAT1B1 simulations, DCF was not stable in the site. The most stable trajectory, 2, resided in the site for 125 ns, but then quickly jumped to 8 Å RMSD, while trajectory 1 exited in under 50 ns, and trajectory 3 at 100 ns (Fig. 10C). These simulations were terminated prior to 500 ns at some point after substrate reorientation. Since DCF is not transported by OATP1B3, it was unsurprising that it unbound in all three cases, but we performed three additional simulations to determine if the site would adjust resulting in a long-lived pose. In all three replicates, 4 to 6, DCF left the pocket. The final snapshot of DCF (pink) from all six replicates is shown on the left of panel B with the initial placement in black. All but trajectory 5 have undergone substantial reorientation breaking any initial contacts with the protein. The binding site residues are mostly conserved between OATP1B1 and 1B3 with key contacts N544, R633, and F356 strictly conserved along with a conservative change Y352F (Y in OATP1B1 and F in OATP1B3). We suspect that the two changes L545S and F386V are the determining factors, as the former places a polar residue against a hydrophobic portion of the molecule, and the volume reduction of the latter fails to provide van der Waals contacts for DCF in the OAT1B3 site.

There is no solved structure of OATP1B3 with DBF bound, but since it is chemically similar to DCF it was straightforward to position DBF in the putative binding pocket of the IF OATP1B3 structure based on the placement of DCF. From here we minimized, equilibrated, and ran three simulations lasting 250 to 300 ns. After equilibration (time zero in Fig. 10D), the molecules had already adjusted resulting in arise in RMSD to just above 2 Å. DBF unbound in all three replicates showing substantial movement in the pocket in the initial 50 ns, while trajectory 3 jumped to 10 Å RMSD on its way toward the cytoplasmic space at 200 ns. Each of these final DBF configurations (pink) compared to the starting template DCF placement (black) are substantial, and the DBF poses themselves are distinct from each other suggesting that an alternate DBF binding pose in OATP1B3 was not captured in our simulations. It could be that the starting DBF configuration was the correct binding site in this transporter, but that unlike the IF OATP1B1 state, the IF OATP1B3 state might be a weak binding state.

Finally, we wanted to determine if the charged residues K41, R181, or E185 played any role in the stability of the substrates in our WT simulations. To aid in this analysis, all of the zoomed-in final snapshots in Figure 10 highlight the three charged residues. For the most part, the molecules do not contact the charged residues except for replicate two in panel C. In this one example, the distance between E185 and DCF drops to ∼3 Å forming a stable interaction for about 100 ns at the end of the trajectory.

Discussion

Proteins fold due to the combination of multiple factors including hydrophobic interactions, disulfide bonds, hydrogen bonds, and ionic interactions. Once a protein reaches its energetic minimum fold, ionic interactions could form and break to control distinct conformational changes that drive isomerization events crucial for the protein function. Our uptake, kinetic, and computational studies, combined with previous modeling and recent structures of OATP1B1/1B3 show that E185, K41, R181, and R580 form a complex salt bridge network in OATP1B1/1B3 that is conformation-dependent (and may as well be involved in the conformational switches of the proteins), and that this interaction is important for transport and dynamics. It was previously shown that OATP1A2 nonsynonymous polymorphs showed significantly lower uptake levels and an order of magnitude slower uptake (Vmax) in the E172D polymorph (36), which is in the equivalent position of E185 in OATP1B1/1B3. Meanwhile the involvement of K41 in the uptake of bromsulphthalein, pravastatin, and E1S were previously studied and mutations of K41 to noncharged amino acids strongly reduced transport activity, while the K41R mutation had uptake levels similar to WT levels (28, 37), although the uptake data were not normalized to cell-surface expression levels. Our data show similar trends in which K41A and E185A lost all uptake activity, while K41R retains some uptake activity with a reduction of Vmax to levels that are almost a third of the WT protein.

Most MFS family members are composed of two 6-TM helical bundles, the N-terminal (TMs 1–6) and the C-terminal (TMs 7–12) bundles that are related by pseudo-2-fold symmetry that runs through the center of the protein perpendicular to the plane of the membrane. The K41-E185 salt bridge links TM1 and 4 of the N-terminal bundle of OATP1B1 and OATP1B3. However, each of the two bundles is made up of two 3-TMs (TMs 1–3 and 4–6 for the N-terminal bundle, and TMs 7–9 and 10–12 for the C-terminal bundle) related by yet another 2-fold symmetry axis that runs through the center of the 6-TM bundle parallel with the membrane plane (38). These domains are believed to be the evolutionary origins of the MFS family (39). The K41-E185 salt bridge holds together the two 3-TM domains of the N-terminal bundle and could be a staple helping to rigidify this 6-TM bundle, thus making it possible for the unit to move in concert during the transport process.

Our simulations do not immediately present an obvious mechanism for why the inner gate remains open in the WT, but the destabilization of the TM1-4 salt bridge in the E185A mutant, at the interface of the broken helical part of TM4, could lead to a larger kink of TM4 and to the closure of the intracellular gate. It is also compelling to note that in many MFS transporters TM1 and TM4 are gating helices that block the bound substrate from the extracellular or cytoplasmic spaces (23). In addition, in the K41E-E185K swap mutant, the reversal of the charges places a positive charge in position 185 that disrupts the native interaction of E185 with R580, leading to loss of activity. R580 is fully conserved in all human OATPs. Weaver and Hagenbuch previously studied the effects of R580 mutations (R580A, R580K, and R580H) on the uptake and kinetics of OATP1B1 (35). The authors hypothesized that R580 might be important for stabilizing a conformation that is required for normal substrate translocation, which is in agreement with our findings for mutations that affect the salt bridge network in which R580 is part of—even though it only interacts with E185 in the OF state. R580 was also studied in the context of OATP1B3 by Glaeser et al. that found small increases in Km values but decreased Vmax for the transport of bromsulphthalein in R580K and R580G (37). Furthermore, the C-bundle R580 in turn also forms a salt bridge with E74 from the N-bundle (Fig. 3). This site has been suggested to be a pivot point in the isomerization of the protein between the OF and IF conformations (17). Hence, the E185A mutation disrupts and influences the 185-41, 185-181, 185 to 580, and indirectly 580-74 charge pairs that are important for the normal isomerization of the protein.

A caveat that is common to all studies that disrupt charged residues is the possibility of altering the water solvation in the binding cavity, by the introduction of neutralizing mutations. Two of our experiments suggest that the loss of activity, rather than the effects of altered hydration is the cause of electrostatics: E185 is conserved in the OATP1B family. E185D would most probably retain its interactions with water molecules in the cavity, but it caused a complete loss of uptake for both OATP-1B1 and OATP-1B3 (Fig. 5). Moreover, the K41E-E185K swap mutant would most probably have a similar cavity hydration as the WT protein, but this mutant also lost all activity, albeit this could have been due to the disruption of the 185 to 580 interaction.

Apart from participating in a salt bridge, the charged K41 and E185 could potentially interact with ligands in the translocation cavity of the OATPs. Earlier MD simulations using homology models, prior to the determination of the OATP structures, showed that the equivalent positions in OATP1A2, K33 and E172, interact with neurosteroids (40). In addition, simulation work on 17β-estradiol 17-glucosiuronic acid binding to OATP1B1 also showed that the ligand interacts with both K41 and E185 (20). In contrast, in the existing OATP structures, the interactions with the four structurally resolved ligands, which are chemically distinct from each other, are primarily hydrophobic with a small number of hydrogen bonds. K41 and E185 do not interact with the ligands in any of these structures (17, 18), and our simulations that are based on these experimental structures are consistent with this even as the ligands unbind and release to the cytoplasm.

Several studies have shown that mutations of certain residues in OATP1B1 differentially affect uptake levels in a substrate-dependent manner. For example, W259A displayed normal uptake and kinetic values for E1S while exhibiting reduced uptake of taurocholate with an 8-fold increase in the Km value compared with that of the WT protein (41), and H54Q resulted in reduced transport and affinity for bromsulphthalein but in higher transport levels for E217βG (26). These differential effects have been noted in other multispecific transporters such as the organic cation transporter 1 where the uptake levels for the F244A mutant varied across a screen of 144 substrates (42). Here, we measured uptake for OATP1B1 with DCF and OATP1B3 with DBF. Although they showed similar uptake levels and kinetic trends across the mutants studies, these two molecules show a high degree of chemical similarity. Nonetheless, our simulations suggest that K41 and E185 impact transport at sites away from the binding sites. Furthermore, we analyzed the pattern of salt bridge pairing for K41, R181, or E185 in the IF and OF small molecule bound states, as we did of the apo structures in Figure 9. In each trajectory, the same pairings occurred as observed for the apo states with E185 interacting with both K41 and R181 simultaneously in the IF states (not shown). Thus, we hypothesize that the impact of K41 and E185 mutants will be insensitive to substrates of diverse chemical structures and would decrease the uptake levels of most OATP1B substrates.

One of the common methods to capture transporters in certain conformational states or to shift the equilibrium of states, either for functional or structural studies, involves the introduction of mutations that alter ionic interactions in the protein. Our approach of simulating salt bridge mutants was successful in capturing large-scale conformational shifts and we intend to employ this in follow up studies to in silico pre-screen mutants for equilibrium conformational shifts for OATPs and other transporters and follow with the promising leads in in vitro studies.

Experimental procedures

Creation of SLCO1B1 and SLCO1B3 mutants

FLAG tagged human SLCO1B1 and SLCO1B3 were restriction cloned into the pIRES2-EGFP vector. Sequences were verified using Eton Biosciences. All mutants were created using site-directed mutagenesis and sequences were verified by Eton Biosciences.

Cell culture and transfection

HEK293T cells (American Type Culture Collection) were cultured in high glucose Dulbecco’s modified Eagle medium (DMEM) (Thermo Fisher Scientific) containing no phenol red and supplemented with 10% fetal bovine serum (FBS) (GeminiBio) at 37 °C and 5% CO_2_. Forty-eight hours prior to transfection, 0.1X10^6^ cells were seeded into each well of a 12-well plate (Costar plates, Corning) that had been precoated with poly-D-lysine (Corning). Transfections were carried out by complexing plasmid DNA and polyethylenimine (linear 25k) (Polysciences) for 15 mins in FBS free DMEM, then adding the mixture to the proper well. Following transfection, cells were allowed to grow for 48 hrs.

OATP1B1 cellular uptake assay

Uptake assays were preformed similarly to methods previously reported (43) with some alterations. Briefly, 48 h post transfection media were removed and cells were washed with DMEM lacking FBS. Cells were then incubated for 5 mins in DMEM lacking FBS at 37 °C and 5% CO_2_. Media were then replaced with fresh media (no FBS) containing DCF (Sigma-Aldrich). OATP1B1 activity assays were carried out using 50 μM DCF for 30 mins while kinetic assays were carried out for 15 mins with varying amounts of DCF. Ligand containing media were then aspirated, and cells were washed 3 times with 1 ml ice-cold phosphate buffered saline (PBS) followed by the addition 250 μl of lysis buffer (1% n-Dodecyl-β-D-Maltoside in PBS). Plates were placed on an orbital shaker at room temperature overnight. The following day, 100 μl of each sample was placed in a 96-well plate and fluorescent measurements (501 nm excitation, 524 nm emission) were carried out using a BioTek Synergy 4 plate reader (BioTek).

OATP1B3 cellular uptake assays

Following transfection, cells were washed with DMEM containing no FBS. After a 5 min incubation at 37 °C and 5% CO_2_ in DMEM, media containing DBF (Acros Organics) was added to the cells (43). OATP1B3 activity assays were conducted with 50 μM DBF for 5 mins while kinetic assays were done for 1 min using different concentrations of DBF. Following the incubation period, media were removed, and cells were washed three times with 1 ml ice-cold PBS. Subsequently, 250 μl of 0.1 mM NaOH was then added, and plates were placed on an orbital shaker at room temperature overnight. Fluorescent measurements were than conducted in a similar fashion to the OATP1B1 samples.

Cell surface expression assay

Cell surface expression of the various OATP constructs was determined as previously described (44) with slight modifications. Briefly, HEK293T cells were seeded at a density of 0.2 X 10^6^ cells into 6-well plates precoated with poly-D-lysine. Cells were transfected as described above. Forty-eight hours post transfection, cells were washed with 2 ml ice cold PBS twice, then incubate for 1 h at 4 °C with 1 ml of 1 mg/ml EZ-Link Sulfo N-hydroxysuccinimide-SS-biotin (Thermo Fisher Scientific). Cells were then washed thrice with 2 ml ice-cold PBS containing 100 mM glycine and then incubated in 2 ml PBS/100 mM glycine for 10 mins at 4 °C. Buffer was then removed and replaced with 1 ml of lysis buffer (50 mM Hepes, 200 mM NaCl, and 1% n-dodecyl-β-D-maltoside) and placed on an orbital shaker for 5 mins until cells detached. The buffer, with the detached cells, was then transferred to 1.5 ml microfuge tubes and placed on an end-over-end rotator for 1 h at 4 °C. Samples were then centrifuged for 2 mins at 10000g, and the supernatant was transferred to fresh 1.5 ml microfuge tubes containing 100 μl of preequilibrated streptavidin-agarose beads (MedChemExpress). Samples were placed on an end-over-end rotator for 1 h at room temperature. Samples were then centrifuged at 850g for 1 min, supernatant was aspirated, and beads were washed thrice with ice-cold lysis buffer. Beads were then incubated with 100 μl of lysis buffer containing 100 mM DTT for 30 mins at room temperature, followed by centrifugation at 850g for 1 min. Supernatant was collected and transferred to fresh 1.5 ml microfuge tubes containing 6x Laemmli buffer. Western Blot analysis was conducted to assess relative OATP expression levels. All OATP constructs were detected using mouse anti-DYKDDDK monoclonal (mAB) IgG (Cell Signaling Technology) (1:1000 dilution) followed by horseradish peroxidase conjugated goat anti-mouse mAB IgG (Novus Biologicals) (1:10,000). CD71 was used as a loading control for normalization of OATP expression and was detected by mouse anti-CD71 mAB IgG (Santa Cruz Biotechnology) (1:1000) followed by horseradish peroxidase conjugated goat anti-mouse mAB IgG (1:10,000). Clarity Western ECL substrate (Bio-Rad) was used for detection with blots imaged by a ChemiDoc MP imaging system (Bio-Rad). Band intensities were quantified using Bio-Rad’s Image Lab software without any downstream modification or processing.

Data processing

Uptake data are presented as mean ± standard deviation. All of the data plotted and analyzed were subtracted with background fluorescence levels of mock transfected cells. Uptake levels were normalized with plasma-membrane OATP levels. Significance between WT and mutant proteins was determined using ANOVA analysis with GraphPad Prism version 10.02 with p values less than 0.05 being interpreted as significant. Kinetic analysis and data fitting were also preformed using GraphPad Prism.

MD simulations

General setup

Simulation systems were prepared using CHARMM-GUI (45, 46). Experimental structures of OF OATP1B1 bound to DCF (PDB ID: 8k6l), IF OATP1B1 (PDB ID: 8hnd), and IF OATP1B3 bound to bicarbonate (PDB ID: 8pg0) (residues 26–281 and 324–653) were embedded in pure 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer (see Table 3). Missing residues (125–136 and 145–169 in 8k6l; 123–137 and 150–159 in 8hnd; 89–90, 124–168, 201–206, 279–281, 449–450, 495–497, 509–517 and 605–609 in 8gp0) were added based on 8hnd (the “most complete” structure) directly by alignment (for small, <5 residue gaps) or with the aid of AlphaFold by using ColabFold (47) with 8hnd as the template structure. Missing termini were not added, nor was the loop joining the N- and C-bundles (residues 282–323). Neutral terminal capping was used, and mutation of E185 to A performed using CHARMM-GUI as appropriate. The systems were solvated with neutralizing 150 mM NaCl in a box size of about 120 x 120 x 136 Å^3^.

DCF was placed in the binding pockets of IF OATP1B1 and OATP1B3 by superposing on the small-molecule bound OF OATP1B1 structure and transferring the DCF to the apo system. Sidechain positions of residues close to the bound DCF were manually adjusted to match the rotamer seen in the initial 8k6l structure (for IF OATP1B1: F356, N541, L545, and S548; for IF OATP1B3: initially only N544 and F356 were adjusted; after unbinding of DCF was observed in all three repeats, further adjustments were made to F352 (Y353 in 1B1), I353, S418, N541, and S548). Coordinates for DBF were generated using CHARMM-GUI. This structure was placed in the binding site of OF OATP1B3 by alignment with DCF, using the latter sidechain-adjusted structure. Forcefield parameters for DCF and DBF were generated from CHARMM-GUI using the CHARMM General Force Field. Total atom count for these systems was between 181,687 and 187,646.

All simulations were performed using Gromacs 2023.3 (48), with the CHARMM-36m forcefield (49) and TIP3P water. Energy minimization and a multistep equilibration gradually reducing restraints on protein and lipid atoms over 5 ns were performed, following the CHARMM-GUI protocol. For simulations containing DBF/DCF, an extra equilibration step was added to allow restraints on the ligand heavy atoms to be reduced after those for the protein backbone; in addition, three separate equilibrations were carried out (for apo simulations, a single equilibration was performed). A simulation timestep of 2 fs was used. Coordinates were saved every 100 ps. Temperature and pressure were maintained at 303 K and 1 atm using v-rescale (time constant 1 ps) and a Parrinello-Rahman semiisotropic barostat (50) (time constant 5 ps, compressibility 4.5 × 10^-5^ bar^-1^), respectively. Hydrogens were constrained using the LINCS (51) algorithm. Nonbonding interactions were cut off at 1.0 nm, with force-switching from 1.2 nm. The Particle-Mesh Ewald (52) method was used to treat long-range electrostatic interactions. The nearest neighbor list was updated every 20 steps with a cutoff at 1.2 nm and a Verlet cutoff scheme.

Unbiased simulations

Repeats from each energy-minimized structure were performed following the equilibration protocol, followed by production runs of various length (see Table 3). Analysis was performed using MDAnalysis (53, 54). Protein RMSD values are computed using Cα backbone atoms; DCF/DBF RMSD values were computed for carbon atoms after an alignment on the protein helix Cα atoms. Bundle angles were calculated by first splitting each bundle into an “intracellular” (residues 26–42, 74–89, 92–104, 182–200, 204–218, 269–278 for N-bundle, 335–349, 385–401, 403–418, 546–562, 567–583, 633–650 for C-bundle) and “extracellular” (residues 42–60, 62–73, 105–118, 169–181, 219–232, 257–268 for N-bundle, 350–368, 370–384, 419–428, 527–545, 584–599, 617–632 for C-bundle) halves. Each bundle is represented by the vector from the extracellular center-of-mass to intracellular center-of-mass, and the angle between them calculated as the angle between these vectors. The angle is set positive if the extracellular halves are closer than the intracellular halves (i.e., an IF conformation), and negative if the intracellular halves are closer (i.e., OF conformation).

PC analysis

PCA was performed using the pca module in MDAnalysis. All MD trajectories were concatenated, and the protein backbone coordinates of the N- and C-bundle TM helices (excluding the loop regions) from all snapshots were aligned to the first frame. Then, we constructed a 3N × 3N covariance matrix for the N selected backbone atoms with respect to the mean protein coordinates in the trajectory. The percent variance explained by each PC was calculated by dividing its eigenvalue by the sum of all 3N eigenvalues. For visualization, we extracted the frames that represent the most extreme projections (lowest and highest weight) along the PC (see insets in Fig. S2A), and performed a linear interpolation between the two (see movies in supplement).

Data availability

Structures referenced in the manuscript that were used for the analysis and the basis of simulations are available from the Protein Data Bank. Additional details, including raw experimental measurements and molecular simulation files, are available from the corresponding author upon reasonable request.

Supporting information

This article contains supporting information.

Conflict of interest

The authors declare that they have no conflicts of interest with the contents of this article.