Metabolomic and transcriptomic analysis reveals endogenous substrates and metabolic adaptation in rats lacking Abcg2 and Abcb1a transporters

Abcg2/Bcrp and Abcb1a/Pgp are xenobiotic efflux transporters limiting substrate permeability in the gastrointestinal system and brain, and increasing renal and hepatic drug clearance. The systemic impact of Bcrp and Pgp ablation on metabolic homeostasis of endogenous substrates is incompletely understood. We performed untargeted metabolomics of cerebrospinal fluid (CSF) and plasma, transcriptomics of brain, liver and kidney from male Sprague Dawley rats (WT) and Bcrp/Pgp double knock-out (dKO) rats, and integrated metabolomic/transcriptomic analysis to identify putative substrates and perturbations in canonical metabolic pathways. A predictive Bayesian machine learning model was used to predict in silico those metabolites with greater substrate-like features for either transporters. The CSF and plasma levels of 169 metabolites, nutrients, signaling molecules, antioxidants and lipids were significantly altered in dKO rats, compared to WT rats. These metabolite changes suggested alterations in histidine, branched chain amino acid, purine and pyrimidine metabolism in the dKO rats. Levels of methylated and sulfated metabolites and some primary bile acids were increased in dKO CSF or plasma. Elevated uric acid levels appeared to be a primary driver of changes in purine and pyrimidine biosynthesis. Alterations in Bcrp/Pgp dKO CSF levels of antioxidants, precursors of neurotransmitters, and uric acid suggests the transporters may contribute to the regulation of a healthy central nervous system in rats. Microbiome-generated metabolites were found to be elevated in dKO rat plasma and CSF. The altered dKO metabolome appeared to cause compensatory transcriptional change in urate biosynthesis and response to lipopolysaccharide in brain, oxidation-reduction processes and response to oxidative stress and porphyrin biosynthesis in kidney, and circadian rhythm genes in liver. These findings present insight into endogenous functions of Bcrp and Pgp, the impact that transporter substrates, inhibitors or polymorphisms may have on metabolism, how transporter inhibition could rewire drug sensitivity indirectly through metabolic changes, and identify functional Bcrp biomarkers.


Introduction
Breast cancer resistant protein (Bcrp; Abcg2) and P-glycoprotein (Pgp; Abcb1) are multidrug efflux transporters with a substantial overlap in their substrates [1,2] and a similar tissue distribution [3]. Bcrp and Pgp transport controls substrate (endogenous substances, dietary constituents, and drugs) absorption from the intestine into blood and from the blood into brain and cerebrospinal fluid (CSF). These transporters also enhance substrate clearance from the liver, kidney and other tissues, thereby regulating tissue exposure of their substrates. Thus, the absence or inhibition of these transporters may lead to accumulation of their substrates, leading to toxicity from endogenous compounds or dietary components or altered efficacy of their substrate drugs. These properties make Bcrp and Pgp among the most important FDA-regulated efflux transporters for drug-approval processes [4].
Therefore, by performing untargeted metabolomics, our goal was to understand the effect of both Bcrp's and Pgp's absence on the CSF and plasma metabolomes, which could help characterize functional biomarkers of these transporters. In addition, to understand the inter-relationship of transporter substrates with tissue gene expression, transcriptomic analysis was performed on kidney, liver and brain frontal cortex from wild type (WT) and Bcrp-Pgp double knockout (dKO) rats. We compared WT versus dKO male rat CSF and plasma because (a) BCRP shows significant sexual dimorphism (male >female) for expression in many tissues [19]; (b) CSF amount is significantly greater in rats vs. mice; (c) dKO animals have a significantly higher CSF concentration of Bcrp/Pgp dual substrates compared to single Bcrp or Pgp KO animals. This is due to upregulated expression of the alternative transporter (e.g., Bcrp) when the other is deleted (e.g., Pgp) that negates the full impact of each transporter to the disposition of dual substrates [20]; and (d) the CSF-to-plasma ratios of Bcrp and Pgp substrates in KO animals are routinely different from those in WT animals, while the plasma concentration does not always vary between WT and KO animals [21]. Thus, sampling only the plasma or using single KO animals would miss identifying some Bcrp or Pgp substrates.
Our metabolomic analysis results demonstrate that Bcrp and Pgp transporters have a significant and distinct impact on plasma and CSF metabolomes. Moreover, machine learning models identified significantly altered metabolites in dKO animals that are candidate Bcrp and/or Pgp substrates. Biological pathway analysis of the transcriptome of WT vs dKO brain, kidney and liver coupled with gene and metabolite enrichment analysis identified potential important biological functions of Bcrp and Pgp efflux transporters. To our knowledge, this is the first

Metabolomic analysis
Untargeted metabolomics of CSF and plasma samples from WT and dKO rats (n = 8 each genotype) was performed by UPLC-MS/MS and GC-MS-based methods at Metabolon 1 (Durham, NC, USA) [22]. Metabolon 1 maintains a metabolite library consisting of retention index (RI), mass/charge (m/z) and chromatographic data of authenticated standard metabolites including the MS/MS spectral information. Named metabolites in the sample were identified by comparing the mass (± 0.005 amu), RI and spectral data of the sample peaks with the library standards, the area under the peak quantified, median normalized and then compared between the WT and dKO groups using the two-sample t-test with Welch's correction. An estimate of false discovery rate (q-value) was also calculated because of multiple comparisons and q<0.1 was considered to have higher confidence in the data. Significant metabolites were identified based on p< 0.05 and highlighted dark red for increase or dark green for decrease in CSF of dKO compared to WT rats. Trending significant is considered when 0.1>p>0.05 and highlighted light red for increase or light green for decrease in CSF and plasma of dKO compared to WT rats. Principal component analysis was used to identify the data pattern [23] in the samples (WT vs. dKO plasma and WT vs. dKO CSF). Results were further analyzed by Random Forest Analysis to identify important metabolic features that contribute the most towards the group separation. The final output of fold change and associated p and q values were log 2 -transformed, and volcano plots created on Tibco 1 Spotfire 1 7.5 (TIBCO Software Inc., Palo Alto, CA, USA) for better visualization and identification of the significant metabolites. Based on the hierarchical clustering analysis (heatmap) by Metabolon 1 , 3 CSF samples were assumed to be contaminated with plasma, and this information was used to calculate 'no contamination' fold change and p-values for the CSF samples, and is listed along with the uncorrected results for all CSF metabolites generated by Metabolon 1 . The non-normalized raw data from Metabolon 1 was used to calculate the ratio of the signal intensity of metabolites that were detected both in CSF and plasma (CSF/plasma ratio). CSF/plasma ratio of the metabolites were analyzed by two sample t-test considering unequal variance to identify metabolites that have significant alterations in CSF abundance in the dKO rats compared to WT rats. Fold change of (CSF/plasma) dKO / (CSF/plasma) WT were log 2 -transformed and plotted against log 2 (t values) on Tibco 1 Spotfire 1 7.5 for further identification and better visualization of metabolites with significantly altered CSF abundance in the dKO rats.

Red blood cell protoporphyrin IX (PPIX) levels
Protoporphyrin IX levels in red blood cells were measured in female WT (n = 10) and Bcrp KO (n = 11) rats by florescence activated cell sorting as described previously [24]. Briefly, peripheral blood was diluted in phosphate-buffered saline and analyzed using a BD LSR II flow cytometer (BD Biosciences, San Jose, CA, USA) at Ex 405 nm and Em 695/40 nm for intracellular PPIX levels. Mean fluorescence intensity of PPIX from 10 to 12 animals in each genotype was analyzed using FlowJo software (Tree Star, Inc., Ashland, OR, USA).

Metabolomic analysis using MetaboAnalyst
Raw metabolite intensity data from Metabolon 1 was further analyzed separately using the MetaboAnalyst software (3.0 and 5.0). Briefly, raw data was first screened for missing values and any metabolite with more than 50% missing values was omitted from further analysis. Due to a low number of features (<2000), the data filtration step was skipped. Raw data was further log-transformed and auto-scaled to obtain a Gaussian data distribution for further statistical analysis. Significant metabolites were identified based on �1.5-fold increase or decrease in median intensity in dKO compared to WT CSF or plasma, and a p-value less than 0.05. Partial least square discriminant analysis and hierarchical clustering analysis were performed to assess the group separation (WT vs. dKO). The list of significant metabolites was compared with the Metabolon 1 data to assess the reproducibility of the significant features.

RNA extraction and microarray analysis
Brain, kidney and liver tissues were weighed, diluted 10 times with TRIzol™, homogenized and RNA isolated according to the manufacturer's instructions. RNA samples were analyzed for quality and integrity on the Agilent 2100 Bioanalyzer and the concentration of RNA was determined using NanoDrop™ 8000 spectrophotometer. Tissues from the same five rats/genotype were analyzed by microarray. 125 ng of total RNA was processed for use on the microarray using the Ambion WT kit and the Affymetrix Terminal labeling kit according to the manufacturer's instructions. The resultant single-stranded cDNA was labeled, fragmented and then hybridized to the Rat Gene 2.0 ST array for sixteen hours at 45˚C at 60 rpm. The arrays were then washed, stained, and scanned using the Affymetrix Model 450 Fluidics Station (Affymetric, Inc. Santa Clara, CA, USA) and Affymetrix Model 3000 7G scanner Station (Affymetric, Inc. Santa Clara, CA, USA) using the manufacturer's recommended protocols. Affymetrix Rat Gene 2.0 ST array was used to perform the gene expression level analysis using Affymetrix Expression Console software and results saved as RMA (Robust Multiarray Averaging) level signal. The transcriptomics data discussed in this manuscript has been deposited in NCBI's Gene Expression Omnibus (GEO) and are accessible through GEO Series accession number GSE159058.
The RMA signal file was uploaded into Nexus expression 3.0 software (BioDiscovery Inc., El Segundo, CA, USA) along with the biological pathway information for the Affymetrix Rat Gene 2.0 ST Array, pooled based on intensity, corrected for multiple test comparisons at q < 0.05 and fold change �1.2, to identify the significantly altered genes and biological pathways.
Gene Set Enrichment Analysis (GSEA) for the transcriptomic data was performed following a previously published method [25].

Integrated pathway enrichment of metabolite and transcriptome
Combined metabolite and transcript pathways of small molecules were constructed from the SMPDB [26]. Differential analysis of genes and metabolites was performed using LIMMA [27]. Pathway enrichment was performed using Fisher's Exact Test with Benjamin-Hochberg multiple testing. We used a p-value of 0.1 and FDR of 25% as cutoffs for this analysis.

In silico prediction of BCRP and/or PGP putative substrates or inhibitors
CSF metabolites identified by Metabolon 1 as significant (p< 0.05) or trending significant (0.05<p< 0.1) were analyzed for their interaction with BCRP and/or PGP by a predictive Bayesian Machine learning model implemented in Assay Central 1 at Collaborations Pharmaceuticals (Raleigh, NC, USA) according to previously published methods [28][29][30]. Models were generated for human PGP (ChEMBL TargetID 3467) using 1217 molecules and for human BCRP (ChEMBL5393) with 392 molecules. A five-fold cross-validation was performed to assess the receiver operator characteristic. The metabolites were then scored with these models to output a score as described previously [30]. Furthermore, physicochemical descriptors, e.g. AlogP, molecular polar surface area, hydrogen bond acceptor and donor count, and molecular weight for the metabolites were obtained from the public database. Regression analysis employing standard least square method was used to predict observed fold changes (dKO/ WT) in the plasma and CSF for each metabolite using the molecular descriptors, along with the ABCG2 or PGP Bayesian model score to determine a structure-property relationship and understand if the metabolites detected in CSF and plasma have any predictive property of BCRP and/ or PGP substrate.

Metabolomics results
Metabolomic analysis identifies differences in CSF and plasma metabolite levels between WT and dKO rats. The Venn diagram depicts the untargeted metabolomics analysis of metabolites in male Bcrp-Pgp dKO and WT rats that were uniquely present only in plasma (306), or in CSF (25), or in both compartments (233) (Fig 1A). Principal component analysis (PCA) (S2 Fig) and Partial Least Square Discriminant Analysis (PLS-DA) (Fig 1B and  1C) identified unique patterns of metabolites in the CSF and plasma of WT vs. dKO rats. Unsupervised hierarchical clustering analysis of all identified metabolites (Fig 1D and 1E) demonstrated that genotype (WT or dKO) segregated all samples in both tissues.
Significantly altered metabolites in CSF and plasma of dKO rats (compared to WT) display unique signatures in each tissue. Among the metabolites significantly altered only in the CSF (71), plasma (79), or both (19), amino acid related metabolites were the predominant CSF-altered metabolites (57%) (Fig 2A), while lipids comprised the foremost (49%) significantly altered metabolites in the plasma ( Fig 3A) in dKO compared to WT rats. Volcano plots showed that while the majority of metabolites in CSF were increased in dKO rats ( Fig  2B), an equal number of metabolites were increased as were decreased in the dKO plasma ( Fig  3B). The metabolite with the largest increase in dKO rat CSF was O-methyl catechol sulfate (5.5-fold) (microbiome-generated [31]) and in dKO plasma was 4-hydroxychlorothalonil (3.6-fold) (a soil degradation product of the organochlorine pesticide chlorothalonil [32]). Only one metabolite, indolin-2-one, was significantly decreased (0.29-fold) in the dKO CSF, while many lipid metabolites were significantly decreased in dKO plasma.
Random Forest Analysis (RFA) identified the metabolites with the greatest effect on group separation (WT and dKO) in CSF and plasma (S3 Fig). Analysis of the metabolites with Meta-boAnalyst, using a more stringent data filtering protocol, identified 29 metabolites altered significantly (fold change � 1.

PLOS ONE
Endogenous metabolites and metabolic pathways altered in the CSF and plasma of dKO compared to WT rats. Altered plasma lipid metabolism. The plasma of dKO rats had a significant decrease in medium and long chain fatty acids, unsaturated fatty acids, dicarboxylate, monohydroxy fatty acids, amide and amino fatty acids and polyunsaturated fatty acids (PUFAs) such as eicosenoate, eicosapentaenoate, docosapentaenoate and docosahexaenoate (Table 1A, S1 and S2 Tables), indicating a significant alteration in the fatty acid metabolism pathway. Pathway enrichment analysis of plasma metabolites confirmed that fatty acid pathways were enriched in the plasma (Fig 4).
Conversely, dKO plasma lysolipids increased (Table 1B), which may result from a decline in complex lipid hydrolysis as supported by low levels of glycerol and monoacylglycerols such as 2-palmitoylglycerol and 1-oleoylglycerol. Indeed, high levels of free carnitine in dKO plasma and modestly low levels of long chain carnitine conjugated lipids may highlight a decrease in fatty acid import into the mitochondria, while reduced levels of medium chain fatty acids may further underlie decreased β-oxidation. A similar trend in dKO CSF was suggested by lower levels of the ketone body 3-hydroxybutyrate (Table 1C), a historical marker of lipid oxidation [33,34].
Histidine metabolism. might be altered because 3-methylhistidine, N-acetyl-1-methylhistidine, trans-urocanate and cis-urocanate accumulated in dKO CSF (S3 and S4 Tables), and corresponded with an enriched histidine metabolism pathway in dKO CSF (Fig 4). 1-methylhistidine is a marker of oxidative damage in some tissues and 3-methylhistidine is associated with protein degradation [35]. Importantly, derivatives of histidine with antioxidant properties including carnosine, N-acetylcarnosine, homocarnosine and anserine were diminished in dKO CSF, and this may predispose dKO rats to increased risk of oxidative stressinduced neurological injury as supported by high levels of methionine sulfoxide (oxidation product of methionine) [36].
Altered purine and pyrimidine biosynthesis pathway metabolites. The purine catabolites xanthosine, urate, allantoin and allantoic acid accumulated in the CSF and plasma of dKO vs. WT rats (Fig 5; S1 and S3 Tables), corresponding to enrichment of the purine metabolism pathway in dKO CSF (Figs 4 and 5). Alteration in purine pathway metabolites could be the direct result of a lack of Bcrp function since Bcrp transports urate and purine analogues [38].
The pyrimidine degradation products uracil, 2-deoxyuridine, pseudouridine and N-acetylbeta-alanine accumulated in the CSF and plasma of dKO rats and may reflect that they are Bcrp substrates or that there is increased pyrimidine degradation (S1 and S3 Tables). Likewise,   Raw signal intensity of the metabolites was used to compute the fold change (dKO/ WT) and p-value was calculated using Welch's two sample t-test to identify metabolites as significant (p<0.05) (in dark green & red) or trending significant (0.1>p>0.05) (light green & red). Table 1A lists the significantly decreased overall plasma fatty acid and lipid signature. Significant increase in lysolipids are listed in Table 1B. Increase in few fatty acid and lipid pathway metabolites were observed in CSF (Table 1C). https://doi.org/10.1371/journal.pone.0253852.t001 the pyrimidine pathway was enriched in dKO CSF (Fig 4). The pyrimidine metabolites orotate and orotidine were also higher in dKO CSF and to a lesser extent in plasma compared to WT rats ( Fig 5; S1 and S3 Tables). This might reflect urate inhibition of UMP synthetase blocking the conversion of orotate and orotidine to UMP in the dKO rats [39].

Fig 4. Pathway enrichment and visualization of significant metabolic pathways from Bcrp/Pgp dKO vs WT rat CSF and Plasma metabolites.
Significance is calculated at p< 0.01; Rank order is created based on enrichment score: (k/m)/(n/N); k = number of significant metabolites in pathway; m = number of detected metabolites in pathway; n = total number of significant metabolites; N = total number of metabolites. https://doi.org/10.1371/journal.pone.0253852.g004

PLOS ONE
Altered levels of methylated nucleotides, and more generally methylated metabolites. The metabolites N1-methylguanosine, 7-methylguanosine, N1-methyladenosine and other methylated metabolites including dimethylarginine and S-methylcysteine, 3-methylhistidine, N-acetyl-1-methylhistidine, and S-methylcysteine as well as the methylated metabolite of quercetin (a known Bcrp substrate [40]) were increased in dKO CSF, suggesting Bcrp may recognize some methylated metabolites (S1 and S3 Tables). Alternatively, high levels of these methylated metabolites may suggest a difference in methylation potential between the genotypes that could have potential epigenetic implications (DNA methylation and gene expression).
Dietary metabolites increased in dKO CSF and plasma. Included riboflavin and the related biochemical metabolite flavin adenine dinucleotide (FAD), while urate was significantly higher only in dKO rat CSF (S5 Fig; S1 and S3 Tables). The dietary Bcrp substrate pheophorbide A [7] was detected only in the plasma of dKO, but not WT rats (S5 Fig). Another dietary Bcrp substrate, the soy isoflavone daidzein and its metabolite daidzein sulfate, were also significantly increased in the plasma of dKO (2.05-and 2.41-fold, respectively) compared to WT rats (S5 Fig, S1 Table). D-pinitol, a natural compound in plants [41], was significantly increased in dKO rat plasma (Table 1A).
The metabolomic analysis showed significant alteration of gut microbiota metabolites in both CSF and plasma ( Table 2). Microbiome metabolites, particularly sulfated metabolites of the amino acid tyrosine, p-cresol sulfate and phenol sulfate [42], decreased 0.45-fold and increased 1.93-fold, respectively, in the CSF of dKO rats, while O-methyl catechol sulfate increased 5.5-fold in the CSF. The microbiome metabolites p-cresol sulfate, equol sulfate, and daidzein sulfate are known Bcrp substrates [43,44], and given Bcrp's high affinity for sulfated metabolites, the other microbiota-generated sulfated metabolites might be likely Bcrp substrates as well.
Likewise, there was a significant increase in sulfated metabolites of drugs, chemicals, food components in both plasma (Table 3A) and CSF (Table 3B), that are also likely Bcrp substrates.
The concentrations of a few xenobiotics were altered in the dKO rat plasma (S1 Table). Among those, ketamine, the drug used to anesthetize the rats to obtain CSF (and that we recently published as a putative dual Pgp/Bcrp substrate [45]), was elevated 1.74-and 1.93-fold higher in dKO plasma and CSF, respectively. The amount of 4-hydroxychlorothalonil, a metabolite of the organochlorine pesticide chlorothalonil [32,46] was elevated 3.6-fold in dKO vs. WT rat plasma.
Metabolites with greater CSF-to-plasma ratios in dKO vs. WT rats. The CSF-to-plasma ratio of the raw LC-MS/MS signal intensity of metabolites identified in WT and dKO rats were compared because Bcrp and Pgp are known to have a significant impact on the brain permeability of their substrates [47] and many Bcrp and Pgp substrates show an increased CSF-tobrain or CSF-to-plasma ratio in transporter KO mice [21,48,49]. Fifty metabolites (Fig 6; S5 Table) had increased CSF-to-plasma ratios in dKO vs. WT rats. Urate, a known Bcrp substrate, had a significantly higher CSF-to-plasma ratio (3.23 fold, p<0.05) in dKO rats, while riboflavin, because it had similar fold increase in both CSF and plasma, had a similar CSF-to-plasma ratio (1.11-fold). Metabolites with significantly increased CSF-to-plasma ratios in dKO vs. WT rats (S5 Table) included catechol sulfate, O-methyl catechol sulfate, methionine sulfoxide, isovalerylcarnitine and pyridoxate. Sixty percent of these metabolites represent amino acid (AA) pathways, with some microbiome-generated AA metabolites (e.g., phenol sulfate from tyrosine and isovalerate [42,50]). Thirteen percent of the metabolites were xenobiotics, such as catechol sulfate (pyrocatechol sulfate) and O-methyl catechol sulfate, with 2.13-and 13-fold higher CSF-to-plasma ratios, respectively, in the dKO rats, and are likely Bcrp substrates given the preference of Bcrp for sulfated metabolites [51]. Pyrocatechol sulfate, a known uremic toxin, can originate from phenolic compounds in the diet, and be further modified by the gut microbiome, and subsequently sulfated and methylated in the liver [51]. 11-dehydrocorticosterone, a glucocorticoid metabolite, had a CSF-to-plasma ratio of 1.4 in dKO rats, and is a likely Pgp substrate [9].
Since it was unlikely that all 169 metabolites significantly altered in CSF and plasma of dKO rats were transporter substrates, an in-silico prediction-based approach was used to determine those metabolites that were most likely Bcrp or Pgp substrates. Bayesian machine-learning models built with compounds known to interact with Bcrp and Pgp were used to score the metabolites based on their level of transporter interaction (high score for higher interaction). Raw data for the significantly altered CSF metabolites that are used for this analysis is presented in S6 Table. Compared to the xenobiotics, many of these metabolites are smaller molecules, and might also have weaker interaction with the transporters, which is also reflected by their lower interaction score. Bcrp or Pgp Bayesian interaction score for each metabolite, along with their molecular descriptors (molecular weight, molecular polar surface area, AlogP etc), were used to predict (see Methods) the observed log 2 -fold change of the metabolites in CSF and plasma. While neither CSF nor plasma metabolites showed a significant correlation with the Pgp score (S6 Fig) CSF metabolites with > 1.5-fold increase in the dKO rats showed a good correlation with the Bcrp score (Fig 7), suggesting many of the scored metabolites are more likely to be Bcrp substrates.

Transcriptomics results
Bcrp and Pgp gene deletion results in an altered transcriptome affecting multiple biological processes in rat kidney, liver and brain. Microarray analysis of WT and dKO rat tissues revealed significant differences in gene expression in kidney (435 genes), liver (398 genes) and brain frontal cortex (240 genes) with 75, 45 and 40 genes commonly altered between kidney and liver, kidney and brain, and brain and liver, respectively. In contrast to a decrease in hepatic Cyp2b2 reported in Bcrp and Pgp single KO rats [52], hepatic Cyp2b1 increased in the dKO rats (S7 Table). Similar to findings in transporter single KO rats, expression of none of the major renal drug efflux transporters (S8 Table) was altered between dKO and WT rats. However, some Slc uptake transporters in kidney and liver had altered expression in the dKO rats (S9 Table) (no results were reported for these Slc genes in the transporter single KO rats), e.g. Slc22a5/OCTN2, a carnitine transporter, which was significantly decreased in the kidney of dKO rats. However, the loss of Bcrp/Pgp did not result in a change in expression of alternative renal transporters that might have indirectly altered metabolites such as indoxyl sulfate and p-creosol sulfate (Oat1 substrate) [14] or creatinine (Oct1 substrate) [53], strengthening the conclusion that these metabolite changes were caused directly by loss of Bcrp/Pgp. Indeed, indoxyl sulfate has recently been identified as a Bcrp substrate [54].
The altered metabolome in Bcrp/Pgp dKO rats had an apparent physiological impact promoting compensatory transcriptional changes in not only individual genes (S10 and S11 Tables) but also in biological processes in dKO vs. WT rats (Table 4A, S12 Table) including

PLOS ONE
urate biosynthesis and response to lipopolysaccharides in brain, oxidation-reduction processes and response to oxidative stress in kidney, and circadian rhythm genes in liver, with many processes altered in all three tissues in the same direction (Table 4B). Cholesterol and steroid metabolic processes were altered in the kidney and liver of dKO rats, respectively (Table 4A, S10-S12 Tables).
The unfolded protein response/ER stress response appeared to be induced in liver and kidney of dKO rats. (Fig 8; S10 and S11 Tables). In the absence of Bcrp/Pgp, hepatic genes upregulated included Caspase 12 (Casp12), an endoplasmic reticulum (ER) membrane cell death initiator [55]; proteome activator subunit 3 (Psme3), a subunit of a proteasome regulator that controls proteasome mediated death; Fam134b, involved in degradation of misfolded proteins; and Marveld1, a marvel domain containing 1 protein that may regulate decay of damaged mRNA; and in response to those changes, induction of Bri3 binding protein, Bri3bp, a protein that in the ER that can decrease apoptosis [56]. Kidneys of dKO rats also showed induction of several proteins in ER stress pathways including activating transcription factor 3, Atf3, a transcription factor involved in ER quality control/protein folding after ER stress; Psme3; Nr4a1/Nurr77, a gene induced by ER stress and that can promote kidney injury [57]; and Uncoupling protein 1, Fig 6. Identification of metabolites with greater accumulation in CSF than in plasma of dKO rats, compared to WT rats. For metabolites identified both in plasma and CSF, a ratio of CSF/plasma metabolite signal intensity was generated for WT and dKO rats. The CSF/plasma ratios of dKO was compared with WT using a two-tailed t-test with two-sample unequal variance. This identified metabolites with significantly higher CSF/plasma ratios in the dKO compared to WT rats at t< 0.05. Log 2 ratio of CSF/plasma in dKO and WT was also calculated for each metabolite and plotted against negative log 10 (of their respective t-values) to generate the volcano plot. The horizontal red line describes the t value of 0.05. The vertical red and green lines designate CSF/plasma ratio in dKO that is 1.5 fold higher or lower than WT, respectively. Metabolites that are on the right side of the vertical red line and on top of the horizontal red line, shows metabolites that have significantly higher CSF/plasma ratio in the dKO over the WT rats. Similarly, metabolites on the left of the green vertical line and above the horizontal red line have significantly lower CSF/ plasma ratios in the dKO over the WT rats. Dark green and dark red circles represent metabolites significantly high or low, respectively in the CSF/plasma of dKO over WT rats. Box plots indicate known Bcrp substrates and other metabolites that are significantly accumulated in dKO rat CSF. https://doi.org/10.1371/journal.pone.0253852.g006

PLOS ONE
Ucp1. It is possible that elevated uric acid is the metabolic cue activating ER stress [58] in Bcrp's absence, and is also consistent with Bcrp normally protecting from apoptotic injury [59].
Altered light and odorant sensing signatures in Bcrp/Pgp dKO rats. A number of genes involved in light sensing and vision were decreased in the dKO brain including retinol binding protein 3 (Rbp3) that transports retinoids between photoreceptors and retinal pigment epithelium; G Protein Subunit Alpha Transducin 2 (Gnat2) that stimulates coupling of cGMP and rhodopsin for vision; Cone-rod homeobox, (Crx) that transcriptionally regulates photoreceptor genes; Phosphodiesterase 6C, (Pde6c) important for cone phosphodiesterase activity; phosducin (Pdc) that is involved in visual phototransduction; guanylate cyclase activator 1a (Guca1a) involved in restoring cGMP for phototransduction [60]; and interphotoreceptor matrix proteoglycan (Impg1, Impg2) (S10 and S11 Tables). One possible explanation is the elevated protoporphyrin IX-conjugates in the harderian gland, in the absence of Bcrp [61], may improve light sensing/night vision [62], (this is hypothetical, but a similar chemical structure has been shown to improve night vision).
Expression of some vomeronasal receptors (Vom1r82, Vom1r83) and olfactory receptors (Olr1668, Olr1394) were altered in dKO kidney (S10 and S11 Tables). Vomeronasal receptors are expressed in human and rodent kidney [63,64] and over half of olfactory receptors have  been identified extra-nasally, including kidney [65]. These receptors can bind directly to metabolites to regulate cell signaling. For example, the olfactory receptor Olfr68 expressed in kidney, binds microbiome-generated short chain fatty acids and regulates blood pressure [66]. Top biological processes were identified by enrichment of significantly altered genes in the dKO rat brain frontal cortex, kidney and liver at p < 0.05. Top 22 common biological processes that were identified by enrichment of significantly altered genes in the dKO rat brain frontal cortex, kidney and liver at p < 0.05.
https://doi.org/10.1371/journal.pone.0253852.t004 potent pheromones known to activate olfactory receptors [67]. In total, these findings suggest that alterations in Bcrp-mediated renal secretion of sulfated chemicals may have a role in influencing signal strength of ligands for renal (and nasal) vomeronasal and olfactory type receptors and their downstream signaling, and perhaps animal behavior.
Comparison of gene expression between WT and dKO rat tissues by microarray identified 'circadian rhythm' pathway as a significantly altered biological pathway. Bcrp/Pgp dKO liver showed a down regulation of Nr1d1 (Rev-Erb alpha), a master regulator of the circadian rhythm pathway [68] that is ligand activated by heme, and by high concentrations of protoporphyrin IX (PPIX). While PPIX was below the limits of detection in the plasma of WT and dKO rats, it was elevated (6-fold) in the red blood cells (RBCs) of Bcrp KO rats (Fig 8). It is possible the chronically elevated PPIX levels and chronically activated Nr1d1 drives the altered rhythmicity of the circadian genes (Fig 8) in dKO vs WT rat transcriptome. RORalpha (Rorα), a gene repressed by ligand activated Nr1d1, and Per2, a gene downstream of Nr1d1 and RORα regulated Bmal1, as well as NFIL3, a transcription factor that helps establish circadian rhythm, are out-of-phase between the dKO and WT rats [69] (Fig 8). Interestingly, the International Mouse Phenotyping consortium (IMPC) phenome database showed that mice lacking Bcrp have abnormal sleep behavior/altered percent wake time (www.mousephenotype.org) [70]. The desynchronization of biological rhythms and altered sleeping time would be consistent with symptoms seen in people with porphyrias.
Integrated plasma metabolome and kidney transcriptome analysis identified the porphyrin pathway as altered in dKO rats. Integrated pathway enrichment analysis identified down regulation of genes associated with porphyrias and porphyrin metabolism in kidney-specifically, heme biosynthetic enzymes uroporphyrinogen III synthase (UROS), uroporphyrinogen decarboxylase (UROD) and coproporphyrinogen III oxidase (CPOX), and Cytochrome c oxidase assembly protein COX15 homolog (COX15) also known as heme A synthase, a gene known to be regulated by porphyrin levels [71] and levels of metabolites in the plasma (S13 Table) and CSF (S14 Table). Presumably, the excess PPIX (Fig 8) in the dKO rats is causing feed-back inhibition of heme biosynthesis. All perturbations (up/down) in mRNAs in liver, kidney and brain and in plasma and CSF metabolites in the same metabolic pathways can be visualized through the hyperlinks in S13 and S14 Tables.
Use of gene expression and known substrate knowledge for understanding metabolomic features. Significantly altered biological pathways and the associated gene expression changes indicated that some metabolomic changes may be due to a direct effect of compensatory pathways, or an indirect effect of many metabolites. Fig 8 and S7 Fig show a few potential mechanisms that could underlie the effect of Bcrp/Pgp on gene-metabolite interactions and biological pathways. Key metabolomic and transcriptomic findings from this study are summarized in Fig 9.

Discussion
The goal of the current study was to examine the effect of Bcrp and Pgp efflux transporter ablation on the rat metabolome and transcriptome in order to better understand their overall endogenous functions, predict drug-induced metabolic alterations/side effects secondary to transporter inhibition by xenobiotics or dietary inhibitors, and predict effects of the wellknown BCRP polymorphism [72] on population physiology and drug treatment outcome. In the absence of Bcrp/Pgp, alterations in some metabolites in plasma and CSF might have occurred for reasons such as they were transporter substrates, or represented compensatory transcriptional response to the altered metabolome, and some may result from reprogramming of the gut microbiome by transporter absence.
Collectively, altered CSF metabolites in dKO rats suggest Bcrp and Pgp may contribute to the regulation of a healthy central nervous system in rats. For example, higher CSF uric acid levels would be expected to be beneficial to the dKO rats as uric acid is thought to be protective to the brain, which is sensitive to oxidative damage. Indeed, CSF uric acid levels are inversely correlated with risk of developing Parkinson's disease, and administration of uric acid is neuroprotective [73][74][75]. Likewise, citrulline and riboflavin, whose levels were elevated in dKO CSF, are reported to prevent neuronal cell death [76], and protect tissue from oxidative injury [77], respectively. Conversely, diminished dKO CSF levels of the antioxidants carnosine, Nacetylcarnosine, homocarnosine, and anserine could predispose the dKO rats to increased risk of oxidative stress-induced neurological injury as supported by high levels of methionine sulfoxide, a marker of oxidative stress [36]. Additionally, the elevations in CSF (not plasma) BCAA could competitively inhibit transport of aromatic amino acids (tyrosine, tryptophan and phenylalanine) into the brain and consequently may alter synthesis of an array of neurotransmitters. Elevated CSF levels of xanthine, urate, hypoxanthine and glutamate may suggest differences in neuronal function in the dKO rats since their elevated CSF levels are often associated with excitotoxicity in neurological diseases [78]. In addition, CSF has been shown to be important for providing and distributing factors important for stem cell proliferation [79]. For example, purines, known Bcrp substrates, directly stimulate the differentiation of some stem cells/progenitors to differentiated cell lineages [80,81].
Elevated uric acid in dKO rats may have disparate physiological consequences in different tissues because uric acid is paradoxically both an antioxidant, but also causes oxidative stress [82,83]. In liver it can cause endoplasmic reticulum stress leading to changes in fatty acid metabolism [58]. Indeed, metabolomic analysis found a significant increase in methionine sulfoxide, a possible biomarker of oxidative stress [36], and 'oxidative stress' was a significantly altered biological pathway in the transcriptomic analysis (Table 4A). In pathological conditions, elevated serum urate is paradoxically known to be positively correlated to gout, cardiovascular and renal diseases, while in the CSF it has been shown to be negatively correlated to Alzheimer's and Parkinson's diseases [84].
Uric acid appears to be a master regulator of multiple biological systems in dKO rats. Elevated uric acid in the dKO rats is likely leading to the increased CSF pyrimidine anabolites orotidine and orotoic acid (orotate) upstream of UMP synthase, since uric acid is a potent inhibitor of UMP synthase [39] that catalyzes the final two steps in pyrimidine biosynthesis. We propose that elevated uric acid elicits a compensatory decrease in hepatic molybdenum cofactor sulfurase, the enzyme required to sulfurate/activate molybdoflavoprotein requiring enzymes such as xanthine oxidoreductase that converts xanthine to urate. Decreased molybdenum cofactor sulfurase would be expected to decrease xanthine oxidase activity and urate formation [85]. To further identify possible biological responses to elevated uric acid, we investigated for overlap in renal transcription changes in Bcrp/Pgp dKO mice with those of another mouse model with altered uric acid homeostasis-URAT1 KO mice [86]. Carbonic anhydrase 3 (Car3) and thyroid hormone responsive SPOT14 protein (Thrsp) were two of the most highly induced genes in kidney of both the URAT1 and Bcrp/Pgp dKO mice. Kidneyspecific Car3 has been shown to increase in models of proximal tubule dysfunction and oxidative stress [87], suggesting that induction of Car3, in the absence of urate transporters, is a biomarker of renal response to elevated uric acid. Additionally, it is possible that the system, sensing the decreased (absent) Bcrp/Pgp, is compensatorily inducing Car3 to try and increase transporter function because another Car isoform has been shown to physically interact with Pgp and increase Pgp ATPase activity and transport activity [88].
However, our data suggest PPARα antagonism as there was a significant increase of some fatty acids in dKO rat CSF (Table 1C). In support of this proposition, humans treated with the PPARα agonist fenofibrate show decreased urinary acylcarnitine metabolites including propionyl-, isobutyryl-, 2-methylbutyryl-and isovalerylcarnitine [97]. Importantly, some of these metabolites were increased in dKO rat CSF (Fig 2; Table 1C, S3 Table). However, incomplete PPAR-mediated fatty acid oxidation might account for this as acetylcarnitines, are, in part, intermediate products of amino acids (also elevated in dKO CSF) and fatty acids [98], and the levels of many amino acids were elevated in the dKO CSF. The decrease in PPARα target genes in dKO rat liver including interferon-induced protein with tetratricopeptide repeats 1 (Ifit1) and carnitine palmitoyltransferase 1a (Cpt1a) [97] and the carnitine transporter OCTN2 [99] (S9 Table) is consistent with PPARα antagonism. Further, nicotinamide, a metabolite increased by PPARα activation [97] was decreased in dKO CSF (S3 and S4 Tables).
Our findings suggest that the absence of Pgp/Bcrp affects the gut microbiota, that in turn, influences both the plasma and CSF metabolome of dKO rats. Pgp KO mice have been shown to have an altered microbiome [101], and we have evidence that Bcrp's absence alters the gut microbiome (manuscript in preparation). Evidence for an altered gut microbiome in the Bcrp/ Pgp dKO rats includes alterations in metabolites generated exclusively by the microbiome including: (a) decreased secondary bile acids deoxycholate (0.55-fold, 0.1>p>0.05), glycodeoxycholate (0.1 fold, p<0.01), taurodeoxycholate (0.57 fold, 0.1>p>0.05) and tauroursodexycholate (0.58 fold, 0.1>p>0.05), (Tables 1A and 2); (b) altered tyrosine metabolites produced by gut microbiota [102]: cresol sulfate was significantly decreased in the CSF while phenol sulfate was increased in the plasma of dKO rats ( Table 2). Phenol sulfate is produced by E.coli species of bacteria, which have been reported to increase in abundance in mice lacking Mdr1a (Pgp) [101,103]; and (c) increased microbiota-generated metabolites of tryptophan (indolelactate, 1.5-fold), phenylalanine and tyrosine (phenylacetylglutamine, 0.46-fold) in dKO rat plasma and CSF (Table 2). While our findings suggest an alteration in gut microbiome metabolites in the absence of Bcrp/Pgp transporters in rats (Table 2), the extent to which the altered microbiome and/or altered transport of these metabolites individually contributes to their increased levels in the dKO rat plasma and CSF is unknown. This is particularly true for sulfated metabolites whose levels were significantly increased in CSF and plasma of dKO rats ( Table 3, S1 and S3 Tables), because many arise from the intestinal microbiome, but may also be Bcrp substrates [51,104].
The combined metabolomics/transcriptomic findings suggest that microbiome metabolites are activating the lipopolysaccharide (LPS) pathway and the aryl hydrocarbon receptor (AhR) in the Bcrp/Pgp dKO rats (Fig 8). DKO CSF metabolites that are increased are early biomarkers of LPS-induced sepsis in rats [105], specifically palmitoyl carnitine, acetyl-L-carnitine, propionyl-L-carnitine, methionine, N-acetyl methionine and methionine sulfoxide, and branched chain amino acids (Table 1C, S3 Table). In addition, transcriptomic analysis found 'response to LPS' as one of the top ten significantly altered biological pathways in liver, kidney and brain of dKO rats (Table 4B). AhR is activated by microbiota-generated tryptophan metabolites/uremic toxins [106] such as indoleacetate and 3-indoxyl sulfate (a Bcrp substrate) [54], and increased levels of at least one of these metabolites in dKO rats could be contributing to ligand activation/Ahr-mediated induction of Cyp1a1, Fos and Jun in kidney and NAD(P)H dehydrogenase quinone 2 (Nqo2), an antioxidant enzyme, in brain (S10 and S11 Tables) [107].
CSF, as opposed to plasma, might be a better matrix to identify Bcrp and/or Pgp endogenous substrates because of the significant role of BBB Bcrp/Pgp on the brain penetration of their substrates [108,109], and the possibility that alternative transporters in kidney or liver offsets their systemic impact for shared substrates. Indeed, when we used in silico Pgp and Bcrp interaction scores generated by a predictive Bayesian Machine learning model, along with molecular descriptors (AlogP, molecular weight, polar surface area) to predict the fold change of metabolites observed in plasma and CSF, only significantly altered metabolites in CSF, and not plasma, displayed a significant correlation between predicted and observed data (Fig 7, S6 Fig).
The metabolomic analysis of Bcrp/Pgp CSF and plasma identified novel exogenous substrates such as ketamine, the drug used for anesthetizing the rats, that was significantly increased in both CSF and plasma of dKO rats and that we recently showed was a putative dual Pgp/Bcrp substrate in mice [45]. Interestingly, there are other studies showing that ketamine, as a receptor antagonist, has diverse effects on endogenous metabolism [110]. Hence, considering ketamine is a fast-acting mitochondrial toxicant, it is possible that the different concentrations of ketamine in dKO vs WT rats might have an impact on the metabolomic changes observed, and that could be investigated in future studies. Our observation that 4-hydroxychlorothalonil, a metabolite of the pesticide chlorothalonil [32,46] was increased in the plasma of dKO rats, indicates possible elevated exposure risk to this pesticide in people with compromised BCRP and/or PGP function.
A previous study of gene expression changes in Pgp and Bcrp WT vs. single KO rats [52], profiled changes in expression of 112 drug metabolism and drug uptake and efflux transporter genes in the liver, intestine, brain and kidney because the Bcrp and Pgp KO models have historically been used to compare the pharmacokinetics of drug transporter substrates with WT animals. The only significant changes they found [52] were increased intestinal carboxylesterase in Pgp KO and catechol-o-methyltransferase in Bcrp KO tissues, while our study did not analyze intestine gene expression and catechol-o-methyltransferase was decreased (P<0.016) in dKO rats. Similar to findings in Pgp or Bcrp single knockout rats [52], dKO rats did not show significant compensatory changes in expression of renal or hepatic uptake (S9 Table) or efflux transporters (S8 Table). There was a significant increase (2-fold) in hepatic Cyp2b1 expression in dKO rats (S7 Table), an enzyme that metabolizes glucocorticoids, and this might help explain the decrease in corticosterone in dKO rat plasma (S1 Table, S7 Fig), which might have been expected to increase in animals lacking Pgp [9]. In total, similar to single transporter KO rats [52] dKO rats showed only modest changes in ADME genes suggesting dKO rats could be used for studies of transporter-mediated pharmacokinetics, particularly for dual Pgp/ Bcrp substrates.
However, unlike our report, the previous transcriptomic analysis of Pgp and Bcrp single KO rats [52] did not comprehensively analyze changes in other genes. The transcriptomic changes in brain, liver and kidney of dKO vs WT rats not only suggested apparent PPARγ activation and PPARα antagonism, likely in reaction to altered levels of receptor ligands, but also that in response to elevated PPIX, porphyrin biosynthesis decreased and the circadian rhythm is out-of-phase. Intriguingly, there was altered expression of light and odorant sensing genes in the dKO rats perhaps in answer to increased PPIX-conjugates and altered levels of vomeronasal and odorant receptor ligands, respectively. Induction of the unfolded protein response/ ER stress response in liver and kidney of dKO rats may have resulted from elevated levels of uric acid, while perturbations in the levels of intestinal microbiota metabolites in dKO rats may have led to activation of the AhR receptor and LPS pathway. In total, these results suggest that Bcrp/Pgp mediated transport of endogenous, xenobiotic and microbiota generated chemicals impacts multiple biological pathways (Table 4A and 4B and Fig 9).
Endogenous drug transporter biomarkers are emerging tools for better predicting drug transport and disposition. Currently, no endogenous probe exclusive for Pgp efflux has been identified [111]. While PPIX is an endogenous Bcrp biomarker, it was apparently restricted to RBCs and not present in the dKO plasma. A number of endogenous metabolites associated with Bcrp function were significantly altered in the dKO rat CSF and/or plasma (e.g., urate, riboflavin). Indeed, a human metabolomics genome wide association study of 7,824 adults [112,113] found plasma urate levels correlated with BCRP rs2231142, the coding variant associated with decreased BCRP expression [72]. Also, multiple studies have reported that patients who carry ABCG2 rs2231142 are more likely to have gout, and hence, higher serum levels of uric acid [114]. A candidate Bcrp dietary biomarker is pheophorbide A (PhA), a chlorophyll metabolite previously identified as a selective Bcrp substrate [7,115] that was only detected in Bcrp/Pgp dKO, but not WT rat plasma (S5 Fig, S1 Table), and that we recently identified as a probe which can be used to test for oral Bcrp drug-drug interactions in vivo [116]. While no selective Pgp metabolite biomarkers were definitively identified in this study, notably the primary bile acid glycochenodeoxycholate was significantly decreased in the plasma of dKO rats, and its levels in human plasma were negatively associated the ABCB1 SNPs [112,113]. Although N-acyl-ethanolamine-type endocannabinoids (e.g., oleic ethanolamide), were recently identified as Pgp substrates [13], oleic ethanolamide was only weakly altered (decreased) in dKO rat plasma (Table 1A) demonstrating it would not be a suitable systemic Pgp biomarker.
This systems level understanding of the metabolomic impact and altered biological function due to lack of function of Bcrp and Pgp transporters may lead to better understanding of the desired / undesired effects of therapies that interact with Bcrp and Pgp. For example, pharmacometabolomic studies have identified changes in metabolites following drug exposure. It is possible that some direct effects and side-effects of medications on metabolism [117] result from inhibition of these transporters-i.e., transporter inhibition is causing the metabolomic changes. It will be of interest in the future to compare the metabolic signatures of some drugs profiled in pharmacometabolomics studies with metabolite changes observed in the dKO rats.
Finally, it should also be considered that the absence of ABCG2 function could rewire drug sensitivity indirectly through metabolic changes. For example, since elevated uric acid and subsequently elevated orotate can antagonize 5-fluorouracil metabolism and toxicity [39], it remains to be tested whether persons with ABCG2 genetic variation or inhibition have altered 5-fluorouracil sensitivity. Similarly, elevated levels of histidine and urocanic acid in the absence of ABCG2 function might enhance methotrexate sensitivity in these same individuals. Methotrexate kills cells by inhibiting dihydrofolate reductase and depleting tetrahydrofolate (THF) required for nucleotide synthesis. Elevated histidine and urocanic acid seen in the absence of Abcg2 function would be anticipated to increase metabolism of THF (a co-factor for formimidoyltransferase cyclodeaminase) further depleting THF and nucleotide synthesis and enhance methotrexate cytotoxic effects [118].
Overall, our study highlighted a significant impact of Bcrp and Pgp efflux drug transporters on the systemic metabolome and transcriptome. Further studies are crucial to establish their contribution on selective pathways, but at a systemic level, our study indicated an essential biological role of these transporters. Most importantly, our results suggest that long-term inhibition of these transporters might have significant consequences with respect to neurological function, gut health and immune-inflammatory status of an individual.   Table. Significantly altered SLC transporter genes in dKO rats compared to WT rats. (XLSX) S10 Table. Significantly altered genes identified by the Nexus expression analysis, at p < 0.05 and fold change > 1.2-fold. (XLSX) S11 Table. Results of analysis of microarray data by GSEA (Gene set enrichment analysis). (XLSX) S12 Table. Biological processes perturbed in Bcrp/Pgp dKO vs. WT rat brain frontal cortex, kidney and liver identified by enrichment of significantly altered genes. (XLSX) S13