Investigating Lipid and Energy Dyshomeostasis Induced by PFAS Congeners in Mouse Model Using Systems Biology Approaches

1. Introduction

Per-(poly)fluoroalkyl substances (PFAS) represent a class of synthetic chemicals, widely used in various industrial products. Their structure consists of a hydrophobic fluorinated carbon backbone with hydrophilic ends. Different PFAS congeners have variance in characteristics depending on the alkyl-chain length [1]. Commonly termed as “Forever Chemicals”, PFAS are known with their low biodegradability, allowing their resistance in ecosystems. Perfluoroalkyls favor the bioaccumulation in protein and lipid-rich organs such as liver [2]. Therefore, hepatoxicity is among the primary effects associated with PFAS exposure, as extensively reported in rodent models [3,4,5,6,7,8]. Despite the lack in definite mechanism of PFAS toxicity, research proposed that PFAS-driven hepatic dysfunction could happen through perturbations in lipid metabolism[9,10], hepatomegaly [11,12], elevation in peroxisomal β-oxidation, hepatocellular hypertrophy, and peroxisome proliferation [12,13,14]. Former in vivo studies of perfluorooctanoic acid (PFOA), a PFAS congener, exposure showed lipid alterations in liver via targeting Peroxisome Proliferator-Activated Receptor-α (PPARα) which acts as a transcriptional regulator of enzymes involved in lipid and fatty acid (FA) metabolism [15,16]. This nuclear receptor is long considered a main molecular target for most PFAS congeners [17,18], and PFOA-driven transcriptional alterations in mice liver was reported to be PPARα-dependent [19]. While transcriptomic data provides a snapshot of gene-level regulation, further investigation is needed through the lens of metabolic networks, where the functional impact of gene expression changes on cellular metabolism can be more accurately addressed [20]. Given the complexity of biological systems, genome-scale metabolic models (GEMs) serve as powerful tools to dissect molecular changes. Through integrating multi-omics data with biochemical knowledge, GEMs enable a systems-level understanding of metabolic alterations that underpin molecular dysfunction associated with external exposure [21]. Central to these models is the gene–protein–reaction (GPR) rule, which maps annotated genes to their corresponding enzymes and the metabolic reactions they regulate. Consequently, alterations in gene expression can be systematically linked to changes in metabolic fluxes, enabling the prediction of functional shifts within the metabolic network [22]. In the present study, we coupled transcriptomics data with mice genome-scale metabolic reconstruction to generate context-specific models and predict mechanism of hepatotoxicity associated with three PFAS congeners, PFESA-BP2, PFOA, and GenX. Our study highlights the use of systems biology approaches to identify metabolic alterations in response to PFAS exposure.

2. Materials and Methods

2.1. Acquisition of Liver Transcriptomics Datasets

We used the NCBI’s Gene Expression Omnibus (GEO) [23] database to extract total and bulk RNA sequencing datasets of mouse liver tissue exposed to PFAS. We chose three studies with the accession IDs: GSE147331, GSE147332, and GSE212294. In the first two studies, male and female BALB/c mice of 10–12 weeks age were exposed to PFESA-BP2 (commercially known as Nafion byproduct-2) at five doses of 0.00, 0.03, 0.30, 3.0, and 6.0 mg/kg.bw (Figure 1A). The mice were exposed for seven days (acute exposure) via oral gavage with ad libtum access to food and drinking water [24]. The third study consisted of male C57BL/6J wildtype and PPARα^−/− mice of 9-10 weeks age exposed to two doses of PFOA (0.05 and 0.30 mg/kg.bw/day), and one dose of GenX (0.30 mg/kg.bw/day). The mode of exposure was through drinking water for 20 weeks (chronic exposure) [19]. The mice had ad libtum access to high-fat diet (45% kcal fat) and drinking water (Figure 2A). In all studies, mice were euthanized following the end of exposure and liver tissues were collected for performing total RNA sequencing. Transcriptomics data for these studies was deposited in GEO database.

2.2. Transcriptomics and Pathway-Based Enrichment Analysis of Metabolic Genes

We obtained the raw gene counts for GSE147331, GSE147332, and GSE212294 datasets and mapped them to the metabolic genes in mouse genome-scale metabolic model, iMM1865 [22]. Differential gene expression analysis was conducted via DEseq2 [25] package in R programming language [26]. Genes with a log FC of 0.25 and adjusted p-value <0.05 were considered as differentially expressed genes (DEGs). Pathway enrichment analysis was performed for DEGs using enrichr-KG [27] and we selected Reactome-2022 [28], GO-(Biological-Process) [29,30], MGI-(Mammalian-Phenotype) [31], and KEGG-Human [32] libraries for the analysis. The output of biological pathways and associated genes were visualized using R package of ggplot2 [33], while Venn diagrams were generated using packages ggVennDiagram [34] and VennDiagram[35]. Principle component analysis (PCA) of VST (Variance stabilizing transformation) transformed normalized gene count was conducted using plotPCA in DEseq2 package. Visualization in boxplots and violin plots was conducted through using tidyplots [36]. For easier illustration of differences between groups in GSE212294, we gave new identifications for the exposure groups: Group 1 (WT-Ctrl), Group 2 (KO-Ctrl), Group 3 (WT-GenX), Group 4 (KO-GenX), Group 5 (WT-low PFOA), Group 6 (KO-low PFOA), Group 7 (WT-high PFOA), and Group 8 (KO-high PFOA). Here, KO represents the PPARα^−/− genotype while low-PFOA and high-PFOA reflect the two PFOA doses of 0.05 and 0.3 mg/kg, respectively.

2.3. Transcriptomics Data Integration with iMM1865 Mouse GEM

To investigate the physiological changes associated with PFAS congeners’ exposure, we integrated liver transcriptomics data with the mouse genome-scale metabolic model (GEM), iMM1865 [22]. The mouse GEM contains 1,865 metabolic genes, 10,612 metabolic reactions, 93 subsystems, and 5,839 metabolites [22]. For generating condition-specific metabolic models of mice, we implemented iMAT (Integrative Metabolic Analysis Tool) [37] algorithm using COBRA toolbox v3.0 [38] in MATLAB 2019a, with the academic license of Gurobi Optimizer v12.0 to solve LP and MILP problems. This resulted in the generation of context-specific metabolic models for all samples considered for this study. We performed flux balance analysis (FBA) for all context-specific models by optimizing BIOMASS reaction as the objective function. Following the approach in Baloni et al [20], we considered reactions with flux value < 1e-06 as inactive (carrying no flux) and these were excluded from the analysis. The reactions belonging to same subsystems were grouped together by taking their average flux and subtracting it from control group to detect subsystem-level perturbations due to PFAS exposure [20]. In case of GSE212294 study, we considered the WT-control group as the comparison reference and hence, we could detect the influence of knocking-out PPARα on the liver metabolic profile along with its role in exposure effects. Projective decomposition normalization of SIGDA [39] was applied to rescale the flux matrix so that both subsystems and exposure conditions are comparable in magnitude, without distorting relative flux patterns. This enabled unbiased clustering and visualization of metabolic subsystem perturbations [39].

Changes in metabolic reaction fluxes were investigated using flux variability analysis (FVA) [40]. Minimum and maximum flux values were computed for reactions within model subsystems involved in metabolism of lipid, fatty acids, cholesterol, bile acids, triglycerides, glycolysis, TCA, oxidative phosphorylation, and ATP production. Subsystems were classified as activated if the fluxes of their associated reactions were higher than those observed in the corresponding control groups (i.e., vehicle control for the PFESA-BP2 study or Group 1 for the PFOA and GenX study). Conversely, subsystems were considered suppressed if the reaction fluxes were lower than those of the respective controls. Flux sampling was performed to explore the entire solution space of feasible flux distributions that satisfy the metabolic model's constraints, including mass balance, reaction bounds, and the objective function [41]. OptGPSampler algorithm was used to generate 1000 flux samples per reaction to evaluate the distribution and variability of possible metabolic states across mouse models. FBA, FVA, and flux sampling were carried out in COBRApy [42], in Python v3.9.6 while using Gurobi Optimizer v12.0 to solve LP and MILP problems.

2.4. Statistical Analysis

DEG analysis was done using Wald test in DEseq2, followed by multiple hypothesis correction of p-values using Benjamini-Hochberg (BH) method [43]. Genes were considered significant if they passed the thresholds of adjusted p-value (FDR) < 0.05 and an absolute log2 fold change (|LFC|) > 0.25. Considering the control group as the comparison reference in reaction fluxes, statistical analyses were performed using unpaired two-tailed t-tests. A p-value < 0.05 was considered statistically significant.

2.5. Data and Code Availability

Transcriptomic data with the accession numbers GSE147331, GSE147332, and GSE212294 were downloaded from GEO database. The normalized gene expression mapped to metabolic genes in iMM1865, the mouse genome-scale metabolic model, along with all analysis scripts, are publicly available on GitHub at: (https://github.com/BaloniLab/PFAS_Congeners_Mouse_GEMs.git).

3. Results

3.1.1. PFESA-BP2 Hepatotoxicity in BALB/c Mice Targets Lipid Metabolism in a Sex- and Dose-Dependent Pattern 

PFESA-BP2 (commonly known as Nafion byproduct-2) is a long-chain PFAS, with a structure resembling PFOS (perfluorooctane sulfonate) and is reported to have a half-life of 296 days [44,45]. We used the mouse genome-scale metabolic model, iMM1865 [22], as a template model for this study and mapped metabolic genes to liver transcriptomics data from GSE147331 (male BALB/c) and GSE147332 (female BALB/c) exposed to PFESA-BP2 [24]. The mapped genes and expression values are provided in Supplementary file 1. Sex appeared to be a major factor for variance, as shown in the principal component analysis (PCA). PCA revealed distinct clustering of male and female along PC1, which explained 90.9% of the variance, while PC2 accounted for 3.7% of the variance and was primarily driven by the dose-dependent effects of PFESA-BP2 (Figure 1B). So, we focused on sex-specific metabolic changes due to PFESA-BP2 exposure and the PCA of sex-specific analysis are provided in Supplementary file 2.

We used control group (0.0 mg/kg) as a reference for the DEG analysis in both datasets. In female mice, 359 upregulated and 384 downregulated genes were identified (Figure 1C). The number of DEGs varied across doses: 0.03 mg/kg dose yielded the fewest DEGs (28 upregulated and 42 downregulated), followed by 0.3 mg/kg (103 up, 90 down), 3.0 mg/kg (182 up, 168 down), with the highest number observed at 6.0 mg/kg (330 upregulated and 336 downregulated) (Figure 1C). Since 29 genes exhibited varied expression across different doses, we focused on the ones that showed consistent up- or down-regulation across all dose levels for the dose-comparison (Figure 1E). Similarly, male mice showed a total of 597 metabolic DEGs across doses (328 upregulated and 287 downregulated) (Figure 1D). Number of DEGs differed corresponding to the dose: 0.03 mg/kg dose got the least number of DEGs (26 up, 18 down), followed by dose 0.3 mg/kg (66 up, 49 down), dose 3.0 mg/kg (120 up and 67 down), and dose 6.0 mg/kg with (272 up, 251 down) DEGs (Figure 1F). For visualizing the dose-comparison, we excluded 18 genes with ambiguous expression across doses (Figure 1F). Together, these findings reveal a clear sex- and dose-dependent transcriptional response in metabolic genes, with higher doses (3.0 and 6.0 mg/kg) inducing the highest number of DEGs in both sexes, highlighting potential sex- and dose-specific molecular signatures.

We performed gene set enrichment analysis (GSEA) of metabolic DEGs identified across both sexes and four exposure doses to identify metabolic alterations due to PFESA-BP2 exposure. Lipid metabolism emerged as the most significantly enriched pathway linked to upregulated genes in both sexes across all doses. Interestingly, it was also the top enriched pathway among downregulated genes in both sexes, except at the lowest dose (0.03 mg/kg) (Figure 1 G, H). In both male and female mice, upregulated genes associated with lipid metabolism were mostly involved in fatty acid biosynthesis and elongation (lipogenesis) such as Elovl1, Elovl6, Fasn, Fabp4, Acly, and Acss2 [46,47], along with cytochrome P450 enzymes of Cyp1a2, Cyp2c55, and Cyp4a31 [48]. Other upregulated genes were involved in fatty acid β-oxidation such as Cpt2, Ehhadh, Abcd1, Hsd17b4, and Slc25a20 [49] (Figure 1I, J; Supplementary file 1). Downregulated genes associated with lipid metabolism in both sexes were involved in glycerophospholipid and triglycerides metabolism which included Agpat2, Dgat2, Gpcpd1, Crls1, and Enpp2. Other downregulated genes belonged to sphingolipid and ceramides metabolism including Fabp5, Ugcg, and Asah1 (Figure 1I, J; Supplementary file 1).

Cholesterol biosynthesis and steroid metabolism were the second most significant pathways in enrichment associated with upregulated genes in both sexes, under all doses except the lowest dose (0.03 mg/kg) (Figure 1G, H). The genes included Lss, Hmgcr, Sqle, Mvk, Hmgcs1, and Dhcr7. Same set of genes along with Fdps, Elovl6, Fasn, Mvd, and Fdft1 were involved in the activation of SREPF (Sterol Regulatory Element-Binding Transcription Factors) signaling at dose 3.0 mg/kg in female mice and dose 0.3 mg/kg in male mice (Figure 1I, J; Supplementary file 1). Interestingly, Cyp7a1 was upregulated in male mice while downregulated in female mice following the exposure (Figure 1I, J; Supplementary file 1). This gene encodes for cholesterol 7α-hydroxylase which is the main enzyme in bile acid biosynthesis from cholesterol [50]. High doses of PFESA-BP2 (3.0 and 6.0 mg/kg) significantly activated PPAR signaling in male mice via upregulating the expression of Fabp4, Acox1, Ehhadh, Me1, Pck1, and Cyp7a1. Activation of PPAR signaling was not observed in female mice although they had a higher number of lipid metabolic genes compared to male mice (Supplementary file 1). These results illustrate that PFESA-BP2 disrupts liver function primarily through modulation of key enzymes and signaling pathways related to lipid metabolism in a sex- and dose-dependent manner.

3.1.2. PFESA-BP2 Hepatotoxicity in BALB/c Mice Is Associated with Energy Dyshomeostasis 

Lipid metabolism plays a critical role in fulfilling body's overall energy demands, with lipids serving as a major energy source alongside glucose [51]. We investigated the influence of PFESA-BP2 on energy metabolism, hypothesizing that such perturbations in hepatic lipids will result into energy dyshomeostasis. Female mice exposed to high dose of PFESA-BP2 (6.0 mg/kg) exhibited significant enrichment in citric acid (TCA) cycle via upregulating the expression of Aco2, Mdh1/2, Idh1/2, and Suclg1/2. Same dose resulted in enrichment of pyruvate metabolism via upregulating the expression of genes including Eno1, Pkm, Pdha1, Pdhb, Pck2, and Me1/2 (Figure 1I, J; Supplementary file 1). Fatty acid β-oxidation was enriched in female mice exposed to PFESA-BP2 doses of 6.0 mg/kg and 0.3 mg/kg, for both upregulated and downregulated genes (Figure 1G, I). Similarly, oxidative phosphorylation was enriched among pathways associated with downregulated genes in female mice exposed to the 0.3 mg/kg dose of PFESA-BP2. Interestingly, purine metabolism was also consistently enriched across all doses and similarly linked to downregulated gene expression, indicating broader perturbations in nucleotide and energy-related pathways (Figure 1G).

Our analysis indicated male mice having significant enrichment of pyruvate metabolism and β-oxidation associated with upregulated genes under exposure doses 3.0 and 6.0 mg/kg, respectively (Figure 1H). Pyruvate metabolic enrichment was modulated by increasing the cellular expression of Aldh3a2, Css2, Me1, Akr1a1, Pck1, Dld, and Pck2 (Figure 1J). Both oxidative phosphorylation and respiratory electron transport were significantly enriched among downregulated genes, under exposure doses 0.3 and 6.0 mg/kg (Figure 1H). Although oxidative phosphorylation was enriched for downregulated genes in both sexes, a larger number of associated genes were affected in male mice. Among these Cox6c and Sdhd were identified as commonly downregulated in both males and females (Figure 1I, J). These findings highlight a dose- and sex-dependent disruption in hepatic energy metabolism following PFESA-BP2 exposure, marked by altered activity in key pathways such as the TCA cycle, β-oxidation, pyruvate metabolism, and oxidative phosphorylation.

3.1.3. PFOA and GenX Exposure Targets Fatty Acid and Lipid Metabolism 

Long-chain PFAS such as PFOA (perfluorooctanoic acid) and short-chain PFAS such as GenX (hexafluoropropylene oxide dimer acid) are widely used PFAS with known environmental persistence and physiological effects. GenX was developed as a substitute of PFOA, however both compounds have been reported to impact various physiological and molecular processes [11,52]. To further investigate PFAS-induced hepatotoxicity, with an emphasis on lipid metabolic alterations, we examined the impact of PFOA and GenX on hepatic transcriptome using mouse liver data from GSE212294 study (Figure 2A) [19]. Principal component analysis (PCA) showed variances across exposure groups, with PC1 and PC2 explaining 48.6 and 20.2% of the total variance, respectively. Mice group exposed to high-dose PFOA (0.3 mg/kg) was distinctly clustered from other exposure groups. As the study involved wild-type and PPARα^-/^- mice, PCA revealed a distinct separation between genotypes irrespective of the exposure (Figure 2B). DEG analysis of PPARα^-/^- mice identified 305 upregulated and 237 downregulated metabolic genes (Figure 2C). Pathway enrichment analysis of upregulated genes in PPARα^-/^- mice indicated solute carrier protein-mediated transport, lipid metabolism, biological oxidation, and amino acid metabolism in the top enriched pathways. For the downregulated genes in PPARα^-/^- mice we identified lipid metabolism, FA oxidation, and FA metabolism in the top enriched pathways (Figure 2D). PPARα is a well-established regulator of enzymes involved in lipid metabolism and FA oxidation; thus, its knockout has been consistently associated with impaired β-oxidation and the development of dyslipidemia [19,53].

We also investigated transcriptomic changes in response to GenX and PFOA exposure in both wild-type and PPARα⁻/⁻ mice. Transcriptomic profiling of GenX and PFOA exposure in wild-type and PPARα⁻/⁻ mice revealed significant differential metabolic gene expression across all groups. GenX exposure resulted in 103 upregulated and 75 downregulated genes in wild-type mice, and 118 upregulated and 124 downregulated genes in PPARα⁻/⁻ mice. Low-dose PFOA induced 87/94 (up/down) DEGs in wild-type and 158/143 in PPARα⁻/⁻ mice. High-dose PFOA triggered the most extensive response, with 287/228 (up/down) DEGs in wild-type and 248/217 in PPARα⁻/⁻ mice (Figure 2E, F; Supplementary file 1). Enrichment analysis revealed that lipid metabolism was significantly enriched among upregulated genes in all exposure groups, with its most pronounced transcriptional changes observed in wild-type mice exposed to high-dose PFOA, followed by GenX (Figure 2G). This is concordant with the findings of Attema et al. [19] where they showed both GenX and high-dose PFOA to cause lipid accumulation in liver of wild-type mice. Moreover, high-dose PFOA and GenX exposure led to significant enrichment of fatty acid metabolism, β-oxidation, and PPAR signaling pathways among upregulated genes in wild-type mice. This effect was not observed in PPARα⁻/⁻ mice and low-dose PFOA exposure (Figure 2G). Some of the targeted genes included Dgat1, Gpat4, Agpat2, Lpin2, Pgs1, and Ptdss2 which are involved in glycerolipids and triglyceride biosynthesis. Other upregulated genes included Lipe, Mgll, Pnpla2, and Abhd5 which participate in lipolysis, and FA oxidation-associated genes such as Acadvl, Acadm, Acads, Cpt2, Ehhadh, Hadha, Hadhb, and Slc25a20 (Supplementary file 1). Transcriptomic changes associated with energy metabolism, particularly oxidative phosphorylation and respiratory electron transport, were enriched only in WT mice exposed to high-dose PFOA. This indicates that elevated PFOA exposure may impair mitochondrial homeostasis and compromise energy production. Peroxisome proliferation was significantly enriched for upregulated genes following the exposure of GenX and high-dose PFOA in wild-type mice. However, this pathway was also found to be enriched for downregulated genes in PPARα^-/^- mice across all exposures (Figure 2G; Supplementary file 1). Interestingly, lipid metabolism, fatty acid metabolism and degradation, along with β-oxidation were the notable enriched pathways among downregulated genes in PPARα^-/^- mice at all exposures. The analysis suggests that both GenX and PFOA likely exert their lipid metabolic effects through PPARα activation, with a dose-dependent response particularly evident for PFOA.

3.1.4. Identifying Perturbations in Energy and Lipid Metabolism Due to PFAS Exposure Using Integrated Genome-Scale Metabolic Models 

Considering the complexity of biological systems, we wanted to understand the mechanism by which PFAS exposure exerts perturbations on energy and lipid homeostasis. We used mouse genome-scale metabolic reconstruction, iMM1865 [22], as template and expression data from GSE147331, GSE147332, and GSE212294 studies to generate context-specific GEM using iMAT algorithm [37]. iMAT allows mapping the gene expression to their corresponding metabolic reactions which can be tested for defining metabolic shifts under various conditions [37]. We obtained 20 context-specific GEMs of female BALB/c, 20 for male BALB/c, 18 for WT, and 18 for PPARα^-/^-, representing in silico models of all mice subjected to different exposure conditions. We performed flux balance analysis and flux variability analysis for these context-specific GEMs using COBRA Toolbox v3.0. We investigated the flux distribution for reactions that were part of energy and lipid metabolism. In this study, a subsystem is considered activated when the flux of its associated reactions is higher than that of the control groups (PFESA-BP2 study) or Group 1 (PFOA and GenX study) and suppressed when the flux is lower than the respective controls. Metabolic flux analysis results from both FBA and FVA are provided in Supplementary file 2.

Under PFESA-BP2 exposure, male and female mice exhibited variation in the metabolic states. In female mice, lipid and energy metabolic subsystems were mostly enriched following the exposure comparing with control group, in a dose-dependent manner. Subsystems-level analysis identified citric acid cycle, oxidative phosphorylation, fatty acid oxidation, and cholesterol metabolism which was consistent with gene enrichment analyses. From our analysis we identified phosphatidylinositol phosphate, aminosugar, and NAD metabolism, along with bile acid synthesis subsystems suggesting additional metabolic alterations. Eicosanoid and glycerophospholipid metabolism were elevated at doses 0.3 and 6.0 mg/kg while suppressed at 3.0 mg/kg dose. Sphingolipid metabolism was activated across all exposure doses, except for 0.03 mg/kg. In contrast, glycosphingolipid metabolism was significantly enriched only at the 0.3 mg/kg dose (Figure 3A, B). Notably, perturbations in metabolic fluxes of glycolysis/glycogenesis, mitochondrial transport, NAD metabolism, and folate metabolism suggest impaired mitochondrial function following exposure (Figure 3A).

In male mice, PFESA-BP2 exposure caused significant activation of pentose phosphate pathway (PPP), purine metabolism, and glycolysis/glucogenesis at dose 3.0, followed by dose 6.0 mg/kg while TCA cycle was activated at all doses. The metabolic flux of oxidative phosphorylation was elevated in all doses while suppressed at dose 6.0 mg/kg. Sphingolipid metabolism was activated in all doses but dose 0.3 mg/kg while glycerophospholipid metabolic flux was activated at all doses (Figure 3C). Glycosphingolipid metabolism was significantly enriched at dose 6.0 mg/kg, followed by other doses except 3.0 mg/kg (Figure 3D).

In PFOA and GenX case studies, subsystems associated with energy metabolism such as TCA cycle and oxidative phosphorylation were activated, particularly in wild-type mice exposed to low-PFOA and GenX, respectively (Figure 3E). Fatty acid oxidation was consistently suppressed at all groups comparing with the control (Figure 3E). Metabolic fluxes for reactions associated with cholesterol metabolism was significantly suppressed in mice exposed to low-PFOA while activated under high-PFOA exposure in wild-type mice. Glycerophospholipid metabolism was significantly suppressed in mice exposed to low-PFOA and GenX. Sphingolipid metabolism was significantly activated for high-PFOA while suppressed for GenX exposure in wild-type mice (Figure 3E, F). Pyruvate metabolism was significantly suppressed following low-PFOA exposure (Figure 3E), consistent with previous findings showing decreased levels of phosphoenolpyruvic acid in liver of PFOA-exposed mice[54]. These findings highlight the value of in silico metabolic modeling in predicting perturbations in energy and lipid metabolism, offering potential metabolic signatures for assessing PFAS-induced hepatotoxicity.

3.1.5. PFESA-BP2 Exposure Causes Activation of Cholesterol Biosynthesis in a Dose-Dependent Manner 

Flux sampling enables a more comprehensive characterization of metabolic behavior beyond the optimal solution [55]. Since we captured lipid dysbiosis associated with PFESA-BP2 exposure, we wanted to investigate specific lipid species. Emerging evidence highlighted PFAS-related cholesterol perturbations, we performed flux sampling to predict flux distribution through metabolic reactions involved in the terpenoid backbone and steroid biosynthesis which contribute to cholesterol production (Figure 4A; Supplementary file 3). We examined the following reactions: HMGCOAS (catalyzed by Hydroxymethylglutaryl CoA synthase); HMGCOARx (Hydroxymethylglutaryl CoA reductase); MEVK1 (Mevalonate kinase atp); SQLEr (Squalene epoxidase endoplasmic reticular NADP); LNSTLSr (Lanosterol synthase); LSTO2r (Lathosterol oxidase); and DHCR72r (7-dehydrocholesterol reductase).

According to our analysis, PFESA-BP2 exposure in both sexes led to activation of reactions in cholesterol biosynthesis, predicting cholesterol accumulation, particularly through elevating the metabolic fluxes of mevalonate and lanosterol biosynthesis (Figure 4A). Interestingly, sex-specific differences in fluxes were observed for key reactions such as HMGCOARx, MEVK1, and DHCR72r. These differences may be explained by the sex-dependent expression patterns of their corresponding enzymes (Figure 4B). The findings suggest that PFESA-BP2 exposure potentially causes cholesterol buildup through targeting mevalonate and lanosterol biosynthesis. However, such predictions need further investigation with experimental studies.

3.1.6. PFAS Exposure Causes Energy Dyshomeostasis via Targeting Carbon Metabolism and β-Oxidation 

Lipids contribute to membrane integrity and are key regulators of energy homeostasis. Fatty acids act as key substrates in modulating TCA cycle and β-oxidation [56]. Former reports highlighted the potential energy changes following PFAS exposure, where PFOA-exposed mice experienced it via AMPK/mTOR regulation [54]. PFOS-exposed mice had alterations in pyruvate metabolism and TCA cycle [10] while GenX-exposed mice exhibited changes in glucose and amino acid metabolism [12]. Since various energy-related pathways were enriched following PFAS exposure in the investigated studies, we wanted to spot potential energy targets of PFAS toxicity.

Therefore, we investigated the flux distribution of reactions involved in central carbon metabolism for mice exposed to PFESA-BP2 (Figure 5A; Supplementary file 3) including pentose phosphate pathway, glycolysis, and the citric acid (TCA) cycle. Metabolic fluxes showed that PFESA-BP2 exerted perturbations in all reactions, with a dose-dependent pattern. No sex-specific differences were captured except in case of PGL (catalyzed by 6- phosphogluconolactonase) and PGM (phosphoglycerate mutase) (Figure 5A; ; Supplementary table 4), potentially driven by the expression pattern of their associated enzymes. Since the reactions plotted here are bidirectional, it is difficult to determine the exact directionality of metabolic shifts. However, the predicted flux distributions, together with enrichment analysis highlight potential metabolic signatures that could guide experimental investigations into PFAS-associated energy dyshomeostasis.

In PFOA- and GenX-exposed mice study, we investigated the metabolic fluxes of mitochondrial β-oxidation (Figure 5A; ; Supplementary table 4). The metabolic fluxes showed shifts with genotype- and dose-dependence. For instance, PPARα⁻/⁻ mice (Groups 2, 4, 6, 8) often show lower flux through mitochondrial reactions ACACT1m (catalyzed by Acetyl CoA C acetyltransferase) and ECOAH9m (Enoyl-coa hydratase, 3-Hydroxy-2-methylbutyryl-CoA forming), highlighting a PPARα dependence. High-PFOA dose (Groups 7,8) showed slightly more metabolic activation of reactions ECOAH9m and BUTt (butyrate transport). Notably, metabolic fluxes of ACRNtm (carnitine/acylcarnitine translocase) were captured exclusively in wild-type mice across all exposures. This aligns with its significantly suppressed gene expression following PPARα knockout, suggesting that the loss of PPARα impairs carnitine transport into mitochondria and may underlie the observed suppression of β-oxidation in PPARα⁻/⁻ mice [19,53] (Figure 5A; ; Supplementary table 4).

4. Discussion

We carried out a systematic study to identify metabolic perturbations in mice, in response to PFAS exposure. We analyzed transcriptome data for long- and short-chain PFAS. The main findings of our study are: (1) PFESA-BP2 exposure caused hepatotoxicity in mice through targeting lipid and energy metabolism via activating SREPF and PPAR signaling. The metabolic perturbations included upregulation in fatty acid metabolism, cholesterol biosynthesis, and β-oxidation while downregulating the metabolism of triglycerides and glycerophospholipids. Such changes were dose- and sex-dependent. A broader-systemic changes targeted the blood brain barrier, causing potential neurotoxicity; (2) Flux simulation of metabolic reactions in PFESA-BP2 exposed mice showed further perturbations targeting energy metabolism such as TCA, oxidative phosphorylation, and pyruvate metabolism while potentially activating cholesterol biosynthesis via targeting mevalonate and lanosterol synthesis; (3) PFOA and GenX exposure caused perturbations in lipid and energy homeostasis via activating PPARα; and (4) Computational predictions elaborated the potential disturbances in mitochondrial β-oxidation via targeting the carnitine shuttle, following PFOA and GenX exposure in a PPARα-dependent manner.

Using transcriptomic data from Lang et al. [24], we found that PFESA-BP2 induces sex- and dose-dependent hepatic lipid dyshomeostasis, presumably regulated by SREBF and PPAR signaling (Figure 1G: J). These findings align with prior studies reporting altered serum cholesterol [59] and triglyceride levels [60] following PFAS exposure. Moreover, PFESA-BP2 induced dose-specific hepatic energy shifts—particularly in TCA and pyruvate metabolism, aligning with prior PFAS studies linking carbon metabolic disruption to TCA, glycolysis, and butanoate pathways [61], suggesting mitochondrial dysfunction and consistent with reported PFAS-related mitochondrial impairments [62]. Liver represents the hub of metabolism and former research implied the influence of its metabolic dysfunction on neuronal complications and cognition [63]. Additionally, downregulation of blood-brain barrier transport genes (SLCs, Fabp5, Abcc9) under high-dose exposure suggests liver-brain axis disruption, supported by FABP5’s known role in docosahexaenoic acid (DHA) uptake and cognitive function [64]. Collectively, our findings highlight the extensive impact of PFESA-BP2 on hepatic lipid and energy homeostasis, suggesting systemic toxic effects that may extend to the brain and potentially impair cognitive function.

High-dose PFOA induced strong hepatic lipid metabolic responses in wild-type mice—including β-oxidation, lipolysis, and triglyceride biosynthesis—which were abolished in PPARα⁻/⁻ mice, underscoring PPARα regulatory role in PFAS toxicity [19]. Our findings are consistent with earlier studies highlighting the positive correlation between serum triglycerides and PFOA exposure [65,66,67]. High-dose PFOA enriched oxidative phosphorylation and respiratory electron transport in wild-type mice, indicating mitochondrial disruption (Figure 2G; Supplementary file 1), consistent with previous reports of PFOA-induced mitochondrial abnormalities [68]. GenX and high-dose PFOA enriched peroxisome proliferation and lipid metabolic pathways in wild-type mice, while PPARα⁻/⁻ mice showed consistent downregulation across all doses (Figure 2G; Supplementary file 1). Additionally, peroxisome proliferation was enriched only in wild-type mice, consistent with classical PPARα activation effects [53]. Our results provide transcriptional evidence, conforming with previous reports of lipid accumulation and mitochondrial abnormalities induced by PFAS [54,68].

While transcriptomic analysis offered valuable insights into the molecular pathways affected by PFAS congeners, the underlying mechanisms driving such toxicity remain incompletely understood. To overcome this limitation, we integrated transcriptomic data with genome-scale metabolic modeling to predict functional metabolic shifts and uncover system-level disruptions beyond gene expression alone. Computational predictions revealed dose- and sex-specific PFESA-BP2-induced shifts, activating TCA cycle, oxidative phosphorylation, FA oxidation, cholesterol, and complex lipid metabolism pathways (Figure 3A-D). Notably, glycerophospholipid activation aligns with lipidomics data showing elevated hepatic phosphatidylglycerols and phosphatidylinositols in mice after PFESA-BP2 exposure [8]. Additionally, same analysis highlighted the potential influence of PFOA and GenX exposure on metabolic fluxes of reactions associated with energy metabolism and lipid pathways through a PPARα-independent manner (Figure 3E, F).

While FBA and FVA provide insights into flux distributions, they focus on optimizing a single objective function under the constraints of the genome-scale model (GEM). However, to capture the full range of feasible metabolic states of individual reaction while satisfying mass balance and flux bounds—we performed flux sampling. Flux sampling revealed PFESA-BP2-induced activation of cholesterol biosynthesis via increased flux through HMGCOAS, HMGCOARx, and MEVK1, with sex-specific differences in DHCR72r, suggesting cholesterol buildup through the mevalonate–lanosterol pathway (Figure 4). While this aligns with our enrichment results, it contrasts with Conley et al. [4] who observed cholesterol reduction at a much higher dose (30 mg/kg). Further simulations showed PFESA-BP2 disrupted carbon metabolism in a dose-dependent manner, targeting glycolysis, TCA cycle, and PPP via G6PDH2r, PGM, and MDHm (Figure 5A). These findings support our hypothesis that lipid dyshomeostasis may extend to energy imbalance, offering new metabolic markers for future studies. Modeling simulations depicted both PFOA and GenX to cause genotype- and dose-dependent changes in mitochondrial β-oxidation disruption, with reduced ACACT1m and ECOAH9m in PPARα⁻/⁻ mice and exclusive ACRNtm activation in wild-type mice, indicating impaired mitochondrial fatty acid transport upon PPARα loss (Figure 5B). Earlier studies have associated PFAS exposure with AC upregulation and enhanced β-oxidation [69,70], which is consistent with our computational analysis in GenX- and PFOA-exposed wild-type mice. Elevation in acylcarnitines imply an incomplete processing of long-chain fatty acids through β-oxidation and a limited availability of intermediates required for TCA modulation [71]. Altogether, this study reveals that PFAS congeners, particularly PFESA-BP2, PFOA, and GenX, induce sex- and dose-dependent disruptions in hepatic energy and lipid metabolism. Through integrative transcriptomics and metabolic modeling, we identified critical metabolic disruptions that offer mechanistic insights and reveal potential biomarkers of PFAS-induced toxicity.