An Organic Anion Transporter 1 (OAT1)-centered Metabolic Network

There has been a recent interest in the broader physiological importance of multispecific “drug” transporters of the SLC and ABC transporter families. Here, a novel multi-tiered systems biology approach was used to predict metabolites and signaling molecules potentially affected by the in vivo deletion of organic anion transporter 1 (Oat1, Slc22a6, originally NKT), a major kidney-ex-pressed drug transporter. Validation of some predictions in wet-lab assays, together with re-evaluation of existing transport and knock-out metabolomics data, generated an experimentally validated, confidence ranked set of OAT1-interacting endogenous compounds enabling construction of an “OAT1-centered meta-bolicinteractionnetwork.”Pathwayandenrichmentanalysisindi- cated an important role for OAT1 in metabolism involving: the TCA cycle, tryptophan and other amino acids, fatty acids, prosta-glandins, cyclic nucleotides, odorants, polyamines, and vitamins. The partly validated reconstructed network is also consistent with a major role for OAT1 in modulating

A great deal of recent evidence suggests that solute carriers (SLC) 2 play a much broader role in physiology, including sig-naling and metabolism, than has been previously appreciated (1,2). Indeed, there has been a recent call for more systematic analysis of the roles of SLCs in metabolism and signaling (1). Perhaps due to heavy emphasis on their key role in pharmacokinetics, members of the SLC (and ABC as well) "drug" transporter families are not generally depicted in biochemical pathways involving the endogenous metabolites they transport (2). Such an omission could have clinical consequences, because drugs directly or indirectly affect pathways normally involved in the movement of key metabolites, pathway intermediates, and signaling molecules, thereby fundamentally affecting cell and organ physiology in normal, pathophysiological, and developmental situations (1,2).
According to the remote sensing and signaling hypothesis, SLC and ABC drug transporters are important in regulating the movement of small endogenous molecules such as key metabolites (e.g. ␣-ketoglutarate, tryptophan metabolites), signaling molecules (e.g. cAMP, prostaglandins, polyamines), vitamins, antioxidants (e.g. uric acid), and certain hormones (e.g. thyroxine) between tissues, organs, and even organisms (2)(3)(4)(5)(6)(7)(8)(9)(10)(11). According to this theory, the SLC and ABC transporters form an integrated network allowing remote communication between different tissues via small endogenous molecules. This integrated network functions in a manner similar to the neuroendocrine system and is, in fact, interlinked with it. The ability of SLC and ABC drug and other transporters to regulate or modulate broad aspects of systemic physiology suggests that drug-metabolite interactions might extend well beyond simple competition for transport at the binding site(s) and provide an explanation for aspects of certain drug-induced metabolic syndromes (e.g. those seen with diuretic use or chronic HIV antiviral treatment (12,13)). Furthermore, elucidation of their physiological role is likely to be useful for further defining the roles of drug transporters in modulating common metabolic diseases, such as diabetes, gout, and chronic kidney disease (2,4,14,15).
Among the aforementioned transporters, organic anion transporter 1 (OAT1/SLC22A6, originally identified as NKT (16,17)), likely the main basolateral probenecid-sensitive organic anion drug transporter of the kidney, mediates ratelimiting steps in the renal elimination of organic anionic drugs and a few cationic drugs (7, 18 -22). This drug transporter has also long been suggested to play a role in key endogenous functions (21,(23)(24)(25), whereas analysis of Oat1 knock-out mice have provided critical information about its potential role in basal physiology (23,26,27). For example, targeted and untargeted metabolomics analyses have revealed significant changes in the concentrations of a number of endogenous metabolites in the knock-out animal, including key biochemical pathway intermediates and signaling molecules (23,27). Moreover, OAT1 (as well as other SLC22 transporters) has also been associated with metabolic abnormalities and disease (2,3,7,11,25,26,28,29). Taken together, this suggests an important, if underappreciated, role for SLC drug transporters in metabolic processes and signaling.
Taking a cue from our previous systems biology efforts to analyze the physiological role of OAT1 (11), we sought to build a detailed map of metabolic and signaling pathways modulated by "drug" transporters such as OAT1. A systems biology approach involving integration of OAT1 knock-out and wild type gene expression data into a genome-scale metabolic reconstruction (GEM) together with constraint-based modeling (i.e. flux variability analysis (FVA)) was used to predict metabolites affected by the absence of OAT1. Pharmacophore-based virtual screening, re-evaluation of existing transport, and knock-out metabolomics data, as well as wet-lab validation were then used to constrain and rank the predicted compounds based on their potential to directly interact with OAT1. Input of this data into the Cytoscape plug-in, MetScape, enabled the generation of a largely experimentally validated, confidence-ranked OAT1centered metabolic interaction network. Pathway and enrichment analysis supported an important role for OAT1 in several key metabolic and signaling pathways. Moreover, the results indicate the feasibility of this novel hierarchical, integrative approach in the context of generally understanding drug transporter-related metabolism, as well as other biological processes involving SLC and ABC transporters. The results appear consistent with the Remote Sensing and Signaling Hypothesis (2-11).

Results
We have previously used GEMs to identify non-obvious, novel OAT1 substrates that were experimentally validated (3,11). However, it was also noted that the GEMs failed to detect some of the known OAT1 substrates. Thus we hypothesized that by using constraint-based modeling (to capture systems level interactions) in conjunction with pharmacophore modeling (to capture molecular substrate-receptor (transporter) binding interactions), the strengths of each methodology could be leveraged to overcome their respective weakness. We performed a loosely constrained, context-specific analysis of transcriptomic data of a global metabolic model to catch all possible metabolites potentially interacting with OAT1 (either directly or indirectly) followed by pharmacophore-based virtual screening, along with comparison to in vivo/in vitro databases and wet-lab assays to provide the specificity filters to differentiate those metabolites that would directly interact with OAT1 from those likely to indirectly be part of the pathway (but not an OAT1 substrate). Ultimately, this multi-tiered systems level analysis, together with metabolomics (from knock-out and wild type animals) and transport data from OAT1-expressing cells was used to build a detailed confidence-ranked OAT1-centered metabolic interaction network that includes many small mol-ecules with well established functions in metabolic and signaling pathways (Fig. 1). We arrived at the final network by working through 5 stages: Stage 1: Prediction of Metabolites Affected by the Deletion of OAT1-The first phase involved a systems biology approach in which transcriptomic data derived from the Oat1-KO mouse was integrated into the mouse GEM, iMM1415 (30), to predict metabolites potentially affected by the absence of this transporter (supplemental Fig. S1 and Table S1). Context-specific wild type (WT) and knock-out (KO) models were constructed from kidney gene expression data using the Gene Inactivity Moderated by Metabolism and Expression (GIMME) algorithm (31) with iMM1415. This systems biology approach simulates the metabolic state of a tissue/organ based on transcriptomic data and the execution of specific biological processes by converting the metabolic network into mathematical equations followed by the application of linear programming to calculate a solution of fluxes (i.e. measurement of rate of production or depletion of metabolites) for each of the reactions. We confirmed that the generated models could produce many of the small molecules requisite for the normal physiologic function of kidneys.
FVA (32), which calculates the boundary points of the steady state solution space, is one of most common constraint-based approaches used in investigations of the effects of gene deletion on a system (33); however, it requires an appropriate objective function (e.g. biomass production, ATP production, substrate uptake) that yields realistic flux distributions (34). However, the current metabolic models do not adequately encapsulate the in vivo metabolic alterations likely resulting from the loss of the drug transporter (OAT1) in a complex organ, and we were forced to apply approximations to the data to gain some understanding of the in vivo role of OAT1 in metabolism. Thus, whereas several objective functions were considered, including the biomass objective function, the ATP objective function and the renal objective function (35), we ultimately decided to use the biomass objective function because it is: 1) comprised of many metabolites involved in core cellular processes including those involved in regulating/modulating cell maintenance; and 2) applicable to non-proliferative cells, which either produce or consume these metabolites. Furthermore, many of the essential metabolites involved in regulating/modulating cell maintenance (e.g. Krebs cycle intermediates, prostaglandins, vitamins, uric acid, and polyamines) are transported by OAT1 and/or other closely related SLC22 family members (2,36,37). (Note that, in the simulations and comparative analyses between the wild type and knock-out models, we are not optimizing for growth, thus we are not assuming that biomass is being maximized, but rather only that the biomass components can be produced by the cells, which we know to be the case, biologically.) To increase the sensitivity of model predictions, we used broad constraints at this initial stage of investigation allowing the inclusion of all possible metabolites that were then ranked in confidence based on the application of pharmacophore screens, as well as available data on their ability to directly interact with OAT1 (from either in vivo metabolomics data or pharmacokinetic data). The resulting WT and KO models consisted of 2233 and 2143 reactions, respectively (supplemental Fig. S1). FVA was used to compare metabolic differences and to find the maximum and minimum flux values for each reaction in the network, enabling the calculation of flux span ratios of the KO over WT for each reaction. Reactions with flux span ratios equal to 1 indicated no change in reaction activities due to the deletion of Oat1; flux span ratios less than 1 implied decreased reaction activities, whereas ratios greater than 1 indicated FIGURE 1. Overall strategy employed. A diagram of the overall multi-tiered hierarchal approach employed in this study. As described under "Results," the approach is comprised of 5 stages: 1) systems biology analysis of transcriptomic data for the prediction of metabolites altered by the absence of OAT1 (supplemental Fig. S1 and Table S1); 2) validation of predictions using metabolomics and kinetic data (supplemental Fig. S1); 3) computational chemistry analysis, generation of pharmacophore hypotheses, and the application of pharmacophore filtering (supplemental Fig. S2 and Table S2); 4) wet-lab validation and identification of new OAT1 metabolites/signaling molecules (Fig. 3); and 5) construction of a substantially validated OAT1-centered metabolic network (Fig. 4, supplemental Fig. S3). increased reaction activities. It was assumed that reactions that require transport of metabolites would change in the OAT1 KO, and 1026 reactions with altered activities were identified, including 321 exchange/transport reactions, which were responsible for the handling of 177 metabolites. After excluding water and other uninformative molecules (see "Experimental Procedures"), 146 metabolites remained (supplemental Table S1), which were predicted to be linked to OAT1-mediated transport.
Stage 2: Validation of the Models Using in Vivo Data from the Oat1-KO, as Well as with in Vitro/ex Vivo Transport Data-Multiple approaches were used to validate the metabolic reconstruction, as well as to provide a measure of confidence in the ability of the predicted metabolites to interact with OAT1. Initially, metabolomics data from the Oat1-KO were interrogated to determine whether experimental observations corresponded with computational predictions. Previous metabolomics profiles identified 36 metabolites (some of which are also signaling molecules), with significantly altered plasma and/or urine concentrations between the WT and Oat1-KO ( Fig. 1, supplemental Fig. S1 and Table S1) (23,27). Among these 36 metabolites, only 19 were actually included in the iMM1415 GEM (supplemental Table S1), and direct comparison with this metabolomics data revealed that 10 of these 19 metabolites were predicted by the GEM analysis (supplemental Table S1). Applying a hypergeometric test to this sample clearly indicates statistically significant enrichment in such metabolites in our population of 146 predicted metabolites (p Ͻ 0.01), which suggests that the systems biology approach integrating transcriptomics data together with loosely constrained FVA modeling makes reasonable predictions of the in vivo metabolic differences between WT and KO.
In addition, because metabolites that interact with OAT1 would be expected to be affected by the absence of this transporter, a search of the literature was performed and a list of 108 metabolites for which kinetic data exists (i.e. K m and/or K i ) indicating interaction with OAT1 was identified. Among these 108 metabolites, 56 are present in the iMM1415 GEM, of which 21 were found within the 146 GEM-predicted metabolites (supplemental Table S1). Once again, applying a hypergeometric test to this sample reveals statistically significant enrichment in OAT1-interacting metabolites in our population of 146 predicted metabolites (p Ͻ 0.01). Moreover, combining the metabolomics results with the kinetic data generated a list of 65 non-overlapping metabolites with wet-lab support (either in vivo metabolomics or in vitro/ex vivo kinetic data) for interaction with OAT1 in the iMM1415 GEM and out of these 65 metabolites 24 were predicted by the systems biology analysis and applying the hypergeometric test indicates significant enrichment in OAT1-interacting metabolites among the 146 GEM-predicted metabolites (p Ͻ 0.01). Taken together, the significant enrichment of the 146 GEM-predicted metabolites in OAT1-interacting metabolites indicates the utility of the loosely constrained systems biology approach utilized in our study.
Stage 3: Pharmacophore Analysis of Metabolites' Potential to Interact with OAT1-As described above, the broad constraints were applied in Stage 1 to identify all possible reactions and maximize the prediction of potential endogenous metabolites and signaling molecules likely affected by the absence of this drug transporter. To generate an OAT1-centered interaction network, this broad list of metabolites was then filtered and ranked by their potential to interact with OAT1. QSAR and pharmacophore modeling have been used to analyze limited sets of OAT drugs/substrates (6,9,18,38,39). Although many drugs appear to be related to metabolites and signaling molecules (40,41), the availability of chemical libraries and computational tools have led to more systematic comparisons of metabolites, natural compounds, and drugs (42)(43)(44)(45). Indeed, pharmacophores based on OAT1 metabolites have previously been used to virtually screen chemical libraries and identify potential inhibitors that have been experimentally validated (27,46). We thus reasoned that the chemical features of known OAT1-transported drugs might be used to rank the predicted metabolites for their potential to interact with OAT1. In addition, some of these metabolites could then be prioritized for later wet-lab validation to assess direct interaction with OAT1. Therefore, OAT1 pharmacophore models based on a large set of well established drug ligands were built.
To construct pharmacophore models for OAT1, 61 drugs having a published K m or K i less than 100 M for OAT1 were selected as "actives" for model building and model validation; two-thirds of the drugs in this group (41 drugs) were used to build the pharmacophore models (training set), and onethird (20 drugs) of the actives was used as a validating set (supplemental Fig. S2A and Table S2). Because the drugs possess diverse chemical structures (consistent with the known multispecific nature of OAT1), they were first clustered into groups using their atomic property fields (APF) (e.g. hydrogen bond donors, hydrogen bond acceptors, SP2 hybridization, lipophilicity, sizes of large atoms, positive and negative charges, etc.) as discriminators (47)(48). Thus, the training actives were grouped into 7 distinct clusters ( Fig. 2A), and pharmacophore models were then built for each cluster based upon the alignment of its members (Fig. 2, B and C). Fig. 2B demonstrates how members of cluster 1 were first aligned, and "pharmacophore model 1" was built to represent the three-dimensional atomic properties shared among the members of that cluster. Then, the pharmacophore models were validated based on the validating set (known positives) and drugs from Drugbank database (serving as true negatives), and a ROC curve was generated (supplemental Fig. S2B). The calculated area under curve was 80.58, which supports the utility of these pharmacophores to identify OAT1-interacting compounds.
The 7 pharmacophores (Fig. 2C) were then used as threedimensional chemical space constraints for virtual screening of the predicted molecules that revealed that 74 of the 146 predicted metabolites satisfied the constraints (supplemental Table S1) and were therefore predicted to have direct interaction with OAT1. Of these 74 metabolites, 18 are known to have direct interaction with OAT1 based on previous experimental in vitro observations (supplemental Table S1). Thus, compared with the original list of 146 metabolites that had 21 metabolites with kinetic data indicating interaction with OAT1 (prior to pharmacophore filtering), the percentage of metabolites known to have direct interactions was enriched about 2-fold (from 14.4% (21 of 146) to 24.3% (18 of 74) after filtering) and this enrichment in OAT1-interacting metabolites was found to be statistically significant by the hypergeometric test (p Ͻ 0.01).
Stage 4: Wet-lab Validation and Identification of Novel OAT1 Ligands-Based on their ability to fit the pharmacophore models, a subset of 8 commercially available metabolites out of the remaining 56 metabolites predicted to directly interact with OAT1 in the pharmacophore screen, but for which there is no wet-lab kinetic data indicating actual interaction with this transporter, were then randomly selected for validation in wetlab transport assays. Of these 8 metabolites/signaling molecules, four were found to interact with OAT1 in transfected cells (Fig. 3, Supplemental Table S1). These metabolites/signal-ing molecules were dihydrofolic acid, palmitoleic acid, 16-hydroxy-hexadecanoic acid, and prostaglandin E 1 with calculated K i values of 93, 200, 13, and 12 M, respectively (Fig. 3); values that are well within the documented range for many compounds shown to interact with OAT1 (supplemental Table S2). These metabolites are important in whole-body physiology, signaling, and cellular metabolism. For example, prostaglandin E 1 , an endogenous vasodilator, serves to increase peripheral blood flow (48), whereas dihydrofolic acid is required to synthesize both purines and pyrimidines. Palmitoleic acid, a longchained fatty acid serving as a potential lipokine, is important in the regulation of lipid metabolism (49).  Table S2 and Fig.  S2). Seven clusters of 3 or more drugs were created and each cluster was used for the generation of a single pharmacophore hypothesis. B, the alignment of drugs from Cluster 1 and the creation of the pharmacophore model: (a) the chemical structures of the drugs in the cluster were aligned; (b) chemical determinants are superimposed on the cluster alignment as three-dimensional pharmacophore features; (c) three-dimensional pharmacophore model of the aligned drugs created and used for virtual screening. C, pharmacophore models that were generated from each of the 7 drug clusters; blue, hydrogen bond donor; red, hydrogen bond acceptor; white, aromaticity; yellow, hydrophobicity; light red, negative charges; light blue, positive charges.
Taken together with the metabolomics and existing kinetic data described above, the number of wet-lab supported metabolites was increased and the application of a hypergeometric statistical test to the overall method for prediction and identification of novel OAT1 metabolites (the combined in silico and in vitro approach), indicates significant enrichment in the number of OAT1-interacting metabolites (p Ͻ 0.01).
Stage 5: Construction of an OAT1-centered Metabolic Interaction Network-The next phase of this study involved construction and analysis of a substantially validated, confidence ranked OAT1-centered metabolite interaction network (Figs. 4, 5, supplemental Fig. S3 and Table S3). To link OAT1 to multiple metabolic pathways, an interaction network was built based on the results of the aforementioned systems biology/ pharmacophore approach and the wet-lab data (validated here or published previously), using Metscape, a Cytoscape plug-in used to construct and visualize metabolic networks based on the KEGG database (50). The resulting broader OAT1-centered metabolic network (see "Experimental Procedures") consisted of a total of 253 metabolites, including 176 experimentally validated and/or computationally predicted and 77 "plus-one" (directly connected) metabolites (Fig. 4). Of these 176 metabolites, 73 had wet-lab support for interactions with OAT1 either by in vivo metabolomics from the knock-out or in vitro assays that were performed in this study or published. Moreover, ϳ3% of the 2272 (i.e. 78/2272) metabolites comprising the MetScape database at the time of analysis were among those for which OAT1-interaction data were available; however, within the 253 metabolites comprising the OAT1-interaction network more than 28% (i.e. 73/253) have been shown to interact with OAT1. This represented a significant enrichment for OAT1-interacting metabolites in the OAT1-centered interaction network (p Ͻ 0.01). These metabolites were thus placed in the group termed "wet-lab" support (supplemental Table S3) and had the highest level of confidence for being part of an OAT1-centered metabolic network based on their ability to interact with OAT1.
Three other groups of metabolites were included, in order of level of confidence (supplemental Table S3). The metabolites with the next level of confidence were those first predicted by GEM and which also passed pharmacophore filters; these were termed "metabolites with high confidence of interacting with OAT1" (supplemental Table S3). Metabolites only predicted by GEM (but having structures such that they did not pass the drug-based pharmacophore filters with high confidence) were classified as "metabolites likely to be affected indirectly." Finally, the remaining plus-one metabolites were termed "OAT1-first neighbor compounds" (Fig. 4). The network revealed that, in the revised OAT1-centered metabolite pathways, the majority of metabolites were interconnected to constitute a main component, and there were also a number of small self-connected components (Fig. 4); network parameters were measured and some are shown in supplemental Fig. S4. The metabolites within the network participated in more than 20 different canonical metabolic pathways, suggesting the broad importance of OAT1 in metabolism. Pathway enrichment analysis of these 253 metabolites was performed using the online bioinformatics resource, Metaboanalyst (51). This overrepresentation analysis provides statistical information on the impact of the metabolites on various pathways (Table 1 and Fig. 5). The affected metabolite/signaling pathways included carbohydrate (e.g. Kreb's cycle, galactose metabolism, etc.), lipid (glycosphingolipid metabolism, bile acid biosynthesis, etc.), amino acid (alanine, aspartate and glutamate, etc.), nucleotide (purine and pyrimidine), and cofactor and vitamin (vitamin A, B2, B3, B5, B6, and B9) metabolism (Table 1; Fig. 5).
The most highly represented metabolic pathways having at least 7 metabolites are shown and ranked according to the validation percentage in Table 1 (validation percentage is equal to the number of "wet-lab supported" metabolites of the pathway divided by the total number of metabolites in the pathway). Among these pathways, the two with the highest "hits with wet-lab support" were tyrosine metabolism and TCA cycle (58.8 and 71.4%, respectively). The TCA cycle is noteworthy, because it includes metabolites known to be classical substrates of OAT1, such as ␣-ketoglutarate, citrate, fumarate, and succinate (2,9,22,52). Also included among these top-ranked pathways was tryptophan metabolism, which had considerable wet-lab support for 8 of 15 metabolites (53.3%), including anthranilate, xanthurenic acid, kynurenine, and indole-acetic acid, which are also putative uremic toxins associated with chronic kidney disease (CKD) (27,53,54).

Discussion
As previously described, systemic deletion of OAT1 was accomplished using a standard homologous recombination approach and thus resulted in the generalized loss of the transporter (23). OAT1 is largely responsible for the uptake of anionic substrates from the blood, and from the viewpoint of organ physiology, any change in the function of OAT1 (e.g. Oat1 knock-out) would concomitantly alter the concentration of . An OAT1-centered metabolic interaction network contains several essential biochemical pathways. Interaction network consisting of metabolites classified according to the level of confidence for OAT1 interaction: 1) wet-lab supported; 2) predicted to interact with high confidence; 3) predicted to be affected indirectly; and 4) first neighbor (plus-one) (supplemental Table S3). The level of confidence of metabolites is reflected by the size and color of nodes (larger and darker nodes have a higher ranking). In the network, some well represented metabolic pathways are shown. For example, purine metabolism is shown as 25 nodes: 5 wet-lab supported (urate, GMP, hypoxanthine, inosine, 3Ј,5Ј-cyclic GMP); 7 predicted to interact with high confidence (guanosine, deoxyguanosine, deoxyinosine, dGDP, deoxyadenosine, dADP, L-glutamine); 2 predicted to be affected indirectly (3Ј, 5Ј-cyclic AMP, L-aspartate); and 10 first neighbors (guanine, dGMP, reduced glutaredoxin, dAMP, adenine, AMP, IMP, xanthosine 5Ј-phosphate, GMP, 5-phospho-alpha-D-ribose 1-diphosphate). metabolites cleared by this transporter not only in the bloodstream, but also in the cells of the tissue in which it is expressed.
In the mouse, the kidney is the predominant site of OAT1 expression and Northern blot and immunohistochemical analyses of kidneys from OAT1 knock-out animals revealed undetectable levels of gene products (i.e. RNA and protein) (23). OAT1 is also found, albeit at much lower levels of expression, in some other mouse tissues, including the choroid plexus, inner ear, eye, brain, and spleen. Thus, whereas we cannot completely eliminate the possibility that these other tissues contribute to the observed metabolic alterations seen in OAT1-KO animals, based on the relatively minor contributions expected from these other tissues, it would seem likely that the observed metabolic alterations in the blood are driven mainly by the lack of kidney-specific expression of the transporter. Moreover, despite the generalized deletion of this important drug transporter, mutant animals are born at expected sex ratios and both male and female mice are viable, appear healthy, and have a normal life expectancy (23).
The OAT1 drug transporter is a focus of regulatory agencies concerned about side effects of drugs due to interaction at the level of the transporter (55,56). For some time, it has been clear that OAT1 and other drug transporters play key roles in regulating levels of endogenous metabolites and signaling molecules (11,21,27,57). For example, the high but shifting embryonic expression of OAT1 and other SLC22 transporters (OAT3, OCT1, URAT1) in the developing nervous system and other developing tissues led to the hypothesis in 2000 that these drug FIGURE 5. Pathway enrichment analysis of the OAT1-centered metabolic interaction network identifies several essential biochemical pathways affected by the absence of OAT1. A, graph of pathway enrichment analysis comparing the -log(p) to the impact on the various pathways for the network metabolites. Some of the affected pathways are indicated with circles colored based on the p value (darker red indicates a lower p value, whereas yellow indicates less significance) and sized based on the impact on the pathway (larger circles have a greater impact). Pathway impact accounts for both the number of affected nodes and its importance with the maximum importance of each pathway being 1. B, the TCA cycle pathway is shown as an example of an affected pathway as it has the highest percentage of hits with wet-lab support (see Table 1). Affected pathway metabolites are highlighted in red and those with wet-lab support for OAT1 interaction have a green border. The number within the rectangle is the KEGG ID for each metabolite, the seven metabolites shown in red are: C00022-pyruvate, C00149-malate, C00042-succinate, C00122-fumarate, C00158-citrate, C00311-isocitrate, C00026-2-oxoglutarate.

TABLE 1 Top pathways affected by OAT1 deficiency as determined by pathway analysis and ranked by the number of hits (%) with wet-lab support
The represented metabolic pathways are ranked according to their percentage of hits with wet-lab support. TCA cycle, tyrosine, and tryptophan metabolism are some of the most well validated pathways represented in the network. transporters transported endogenous small molecules that could affect morphogenesis (21). There has recently been a call for more systematic analyses of the roles of these transporters, particularly the SLCs, in metabolism and signaling (1). There is also growing evidence showing relevance of OAT1 to metabolic disease, including chronic kidney disease (15,25,58). Taken together with the finding that these highly conserved drug transporters are differentially and highly expressed in various epithelial tissues lining body fluid compartments, it has been hypothesized that they potentially participate in remote communication ("remote sensing and signaling") between organs and organisms (2-4, 7, 9, 10, 27). This may apply to other SLC families as well (1). The remote sensing and signaling hypothesis for SLC and ABC drug transporters, which are mainly found on the apical and basolateral surfaces of epithelial cells lining body fluid compartments (e.g. blood, CSF, bile, amniotic fluid), theorizes that, together, these transporters function analogously to the endocrine system (2,4,7,10). As we show here, they play a key role in regulating or modulating biochemical pathways involving a wide variety of small endogenous molecules with high informational content from the perspective of systemic or local metabolism and signaling. Our analysis thus builds upon prior work (7,11,23,27) implicating OAT1 in the in vivo regulation of pathways involving essential metabolites (e.g. TCA cycle intermediates, tryptophan metabolites derived from the gut microbiome), key signaling molecules (e.g. prostaglandins, polyamines, cyclic nucleotides), molecules with antioxidant activity (e.g. ascorbic acid, urate), hormones (e.g. thyroxine), vitamins, and cofactors (e.g. panthothenic acid) (Fig. 5, Table 1). Crucially, from our analyses we are able to construct a partly validated OAT1-centered metabolic network (Figs. 4, 5; Table 1).
The molecules in the network are central to classical metabolism and signaling pathways, as well as organ and systemic physiology. Moreover, they are important in pathophysiological states like hyperuricemia, metabolic syndrome, diabetes, and chronic kidney disease. For example, some metabolites known to be transported by OAT1 are potential classical uremic toxins (e.g. indoxyl sulfate, kynurenate, polyamines, and uric acid), which accumulate as renal function declines in CKD (53,54). Moreover, many of these pathophysiological states have been linked to SNPs or altered expression and/or function of members of the OAT subfamily of SLC22 (2, 4, 7, 10, 14, 15, 25, 58 -61).
An important corollary has to do with the role of OATs in drug-metabolite interactions due to competition for transport. This may also be important in chronic kidney disease, where uremic toxins such as those mentioned above are transported by OAT1. For example, polyamines can inhibit transport of the drug methotrexate by OAT1 (11). If the competition is more substantial, as might occur with the drug probenecid (which binds OAT1 with high affinity), our data suggests that the metabolic changes could be quite significant, potentially resulting in many metabolic alterations that substantially overlap with the OAT1 knock-out noted here. In this regard, it is worth noting that there are a number of drug-induced metabolic syndromes that are associated with OAT1-transported drugs (12,13).
Furthermore, the complexity of the OAT1-centered network points to the possibility of unexpected metabolic changes that could go well beyond the relatively straightforward concept of transporter-level competition for the ligand binding site. For example, in this study, we were able to separate likely direct versus indirect interactions of metabolites with OAT1; thus, a drug that tightly binds OAT1 may not only alter metabolites that directly compete for transport but also others in the OAT1-centered network that are not directly transported by OAT1. Because thiazide diuretics and HIV antivirals are transported by OAT1 (6, 62, 63), the OAT1-centered network may help in understanding the drug-induced metabolic syndromes associated with chronic treatment with these drugs (64 -67).
In summary, we have described and validated a novel approach for systems driven discovery through integration of distinct and complementary systems biology, computational and wet-lab approaches to gain further insight into the substrate specificity and function of the OAT1 transporter in basal physiology (Fig. 1). As more wet lab transport data accumulates, and with further improvement of pathway analysis and other computational tools-including metabolic reconstruction approaches (particularly with respect to mammalian transporters), it should be possible to further refine and validate the OAT1-centered network proposed here. Furthermore, based on our results, we suggest that the types of multifaceted analyses described here for OAT1, enabling the construction of a "drug transporter"-centered metabolic network (Figs. 4, 5; Table 1), can be applied to other SLC and ABC drug transporters to generate a more comprehensive picture of the role that these transporters play in metabolism. Eventually, this approach could potentially connect cellular metabolism in different organs via molecules transported by multispecific drug transporters (as well as other transport systems such as those involving other types of transporters or channels) as proposed in the remote sensing and signaling hypothesis (2,4,7,9,10).
Importantly, once drug-transporter metabolic networks are created for other SLC and ABC multispecific transporters, the systems biology approach employed here may be useful for explicitly predicting the metabolic alterations expected for new drugs in healthy or diseased populations with globally altered metabolism (e.g. CKD, liver disease, metabolic syndrome, diabetes, multiorgan injury).

Experimental Procedures
The overall approach taken in the study is depicted in schematic flow charts (Fig. 1), which consists of the following stages: 1) systems biology analysis to predict metabolites affected by the deletion of OAT1 (supplemental Fig. S1); 2) validation of models using in vivo and in vitro data (Fig. 1); 3) computational chemistry analysis (Fig. 2, supplemental Fig. S2); 4) wet-lab validation (Fig. 3); and 5) construction of an OAT1-centered metabolic interaction network (Fig. 4, supplemental Fig. S3).
Analysis of Transcriptomic Data with Mouse GEM Network-The transcriptomic and metabolomics data from previous studies (23, 27, 63, 68 -70) were used for the contextspecific analysis. The mouse GEM, iMM1415 (30), which contains the biochemical transformations for numerous tissues and cells in mice, was utilized (see supplemental Fig. S1 for more details). To create context-specific models, gene expression data from wild type (WT) and Oat1 knock-out (KO) mice were analyzed by Microarray Suite version 5.0 to assign present/absent calls. A minimum of 3 sets of microarray data for each condition (WT and KO) was analyzed separately, and for a gene to be considered present, it had to be present in at least 2 of 3 sets of data. In this way, the gene expression data were converted into a "binary" classification (i.e. either absent or present). This binary data without regard to actual expression values was then integrated into the genome-scale metabolic reconstruction model using the GIMME algorithm using the default setting of 90% for the percentage of the objective function needed to be met to generate the model (31).
GEM reconstructions and the application of constraintbased analysis of metabolic networks are widely used in systems biology approaches. The initial planned approach for this study described here was to apply statistical analysis or sampling of the networks with FVA. Unfortunately, the use of strictly constrained FVA models not only resulted in the prediction of a limited number of metabolites with which to work, but they also failed to capture some important interactions known to involve the transporter. Although this is not an uncommon occurrence, it is particularly problematic for this study when one takes into account the fact that our current understanding and the scope of the current metabolic models do not adequately encapsulate the metabolic interactions/alterations resulting from the deletion of SLC drug transporter. In other words, the initial approach was limited by a metabolic model that was not specifically designed to address the questions that were being asked and thus required a modified approach described below.

Sensitivity Filter Through Comparison of the Metabolic Differences between KO and WT with FVA and Model
Validation-FVA was used to compare functional differences between the WT and KO mice based on maximal achievable reaction flux ranges (32). A reduced model was created by removing those reactions classified as absent based on gene expression data from the WT and KO conditions. When establishing the lower and upper bound for exchange reactions, uptake constraints for WT and KO models were set to be the same (10 mol/h) for most metabolites, and to make the model more renal-specific, a few metabolites that were listed for the previously published renal objective function, which was used to analyze blood pressure regulation (35), were set to be either secreted or absorbed accordingly. In our analysis, a biomass maintenance function, which was defined to represent the metabolic composition necessary for the maintenance of mouse tissues, was used as the constraint to generate the context-specific models using the GIMME algorithm. We also considered the ATP objective function, as well as the so-called "renal objec-tive function" (35). As with the biomass objective function, application of these other objective functions to the data requires some degree of approximation. Although the ATP objective function produced similar results to the biomass objective function, the application of the later produced a broader range of results. Interestingly, despite its name, the application of the renal objective function resulted in a more limited set of metabolites for subsequent analysis. This is likely due to the fact that it was designed for the analysis of the role of the kidney in regulating/modulating blood pressure.
The biomass pseudoreaction was ultimately selected to ensure that the cells could produce the complement of small molecules requisite for its function and cellular maintenance as well as metabolites/signaling molecules known to be transported by OAT1 have roles in cell maintenance and growth. Flux spans were then calculated as the differences between the maximum and minimum of reaction fluxes and the pairwise ratio between the KO and WT flux spans were calculated (i.e. , resulting in a set of flux span ratios for the set of reactions shared between the two models. Flux span ratios with a value of 1 were considered to be unchanged; in this loosely constrained model, any variation from 1 was considered to be affected by the absence of OAT1 and a successful model was expected to generate a list of "GEM-predicted OAT1 metabolites" that could be further evaluated (supplemental Table S4).
Specificity Filter through Pharmacophore Model Building and Validation Based on Drugs Known to Interact with OAT1-Computational chemistry analysis was performed with ICM software developed by Molsoft L.L.C (San Diego, CA). The software was used to perform clustering, alignment, and pharmacophore building based on the APF (47) of OAT1-interacting drugs or tracers (pharmaceuticals) with known K m (substrate affinity) or K i (inhibitory affinity) of less than 100 M were selected. A total of 61 pharmaceuticals were selected; among which two-thirds (41 drugs) were selected for a training set for the model generation, and one-third (20 drugs) were placed in a validating set. The structure-data files (pubchem.ncbi.nlm.nih. gov) for the OAT1-interacting drugs were input into ICM, which superimposed and aligned the drugs based on using the APF superimposition method. This method takes into consideration a number of three-dimensional structural parameters, hydrogen bond donors, hydrogen bond acceptors, SP2 hybridization, lipophilicity, size of large atoms, and positive and negative charges (71,72). Based upon these alignments, 7 pharmacophore hypotheses were generated. After the model generation, the validating set (serving as true positives) and all the drugs from the Drugbank database (serving as true negatives) were screened against the pharmacophores. Using these screening results, a ROC curve was generated (supplemental Fig. S2).
Screening the List of Metabolites Predicted by the Metabolic Networks with Pharmacophore Models-The GEM-predicted OAT1 metabolites were compared with each of the 7 pharmacophore models and ranked by how well they fit with the threedimensional molecular space defined by the drugs known to interact with OAT1. The 30 metabolites that best fit each of the 7 pharmacophores were selected as having potential to directly interact with OAT1. Many metabolites fit more than one pharmacophore model. After eliminating overlaps, 74 metabolites remained, and these were termed "top-ranked GEM-predicted OAT1 metabolites" (passing OAT1 pharmacophore filter).
Uptake Inhibition Assay-Oat1-transfected CHO cells cultured on 96-well plates were incubated in the presence of 10 M 6-carboxyfluorescein (6,9,62) with or without individual metabolites (with controls treated with the OAT1 inhibitor probenecid). The IC 50 curves for novel ligands were plotted in Prism Software (GraphPad Inc., San Diego, CA), and the IC 50 values were converted to K i (inhibition affinity) using the Cheng-Prusoff equation.
Construction of an OAT1-centered Metabolic Interaction Network-The OAT1-centered metabolic interaction network was generated based on the categorization and ranking of metabolites according to the confidence for their interaction with OAT1. Those metabolites for which kinetic data exists indicating the ability to interact with OAT1 were merged together and labeled as "wet-lab supported" metabolites. Metabolites predicted by the GEM analysis, which also fit the pharmacophore filtering were designated as "metabolites predicted to interact with high confidence." Metabolites that were predicted by the GEM analysis but which did not fit the pharmacophore screen were designated as "metabolites predicted to be affected indirectly." The KEGG IDs for all of these metabolites were input then into Metscape, a Cytoscape plug-in (50). Metabolites present in Metscape were used as "input metabolites" to build a metabolic network. The construction of the network also introduced many "plus-one" or "first-neighbor" metabolites, and because there is currently no evidence supporting their interactions with OAT1, they were deemed to have the lowest confidence for interacting with OAT1. The network was then trimmed to eliminate uninformative nodes using the following criteria: small molecules (such as water, carbon dioxide, etc), energy-related molecules (NADH, ATP, etc.), and large peptides not known to interact with OAT1 or related transporters (somatostatin, kinetensin, etc.) were removed; unnecessary "dead-ended" and "inter-connecting" plus-ones were removed to create a more concise network (in other words, dead-ended plus-one nodes, which connected only to one node, were removed, as were inter-connecting plus-ones that did not affect connections between wet-lab validated or predicted nodes). The final network thus consisted of metabolites that fell into four categories, which in the order of level of confidence of their potential to interact with OAT1, the categories were: 1) wet-lab supported; 2) predicted to interact with high confidence; 3) predicted to be affected indirectly; and 4) plus-one after trimming.
Metabolic Pathway Analysis-Pathway and enrichment analyses were performed using Metaboanalyst 3.0 for pathway analysis and visualization (51). Lists of the KEGG IDs for the metabolites were input into this online bioinformatics resource to either the Pathway Analysis or Enrichment Analysis func-tionalities on the MetaboAnalyst 3.0 website. For the pathway analysis, the Homo sapiens (human) pathway library was selected and all compounds in the selected pathway were used. The algorithms specified were the hypergeometric test for the over-representation analysis and the relative betweeness centrality for the pathway topology analysis. For enrichment analysis, the pathway-associated metabolite set was selected as the library and all compounds in the metabolite set library were used.
Statistics-To determine whether the overall in silico approach results in significant enrichment of metabolites known to have direct interaction with OAT1, a hypergeometric-based test was performed to calculate the various p values. The hypergeometric calculation, which is based on certain assumptions, has been used in systems biological analyses for determining the probability of a result occurring just by chance (73)(74)(75). The hypergeometric test can be used as a measure of over-representation and takes into account the overall population size, the number of successes within this population, the sample size, and the number of successes within the sample population.