Supplementary data for: Comparison of phenotypic and transcriptomic profiles between HFPO-DA and prototypical PPARα, PPARγ, and cytotoxic agents in wild-type and Ppara-null mouse livers

Recent in vitro transcriptomic analyses for short-chain per- and polyfluoroalkyl substances (PFAS), HFPO-DA (ammonium, 2,3,3,3-tetrafluoro-2-(heptafluoropropoxy)-propanoate) added to the weight of evidence supporting the peroxisome proliferator-activated receptor alpha (PPARα) activator-induced hepatocarcinogenesis mode of action (MOA) for HFPO-DA-mediated liver effects in rodents. Importantly, PPARα-mediated key events (KEs), including hepatocellular hypertrophy and proliferation that have been shown to occur prior to tumor development in this MOA, are rodent-specific and likely not human-relevant. To further inform the MOA of HFPO-DA and evaluate other hypothesized MOAs, phenotypic and transcriptomic responses in wild-type (WT) and Ppara-null mice were investigated following short-term exposure to HFPO-DA or prototypical agonists of PPARα (GW7647), PPARγ (rosiglitazone), or cytotoxicity (acetaminophen). Phenotypic and transcriptomic assessment of mouse livers demonstrated a general lack of response to HFPO-DA or GW7647 exposure in Ppara-null but not WT mice. Conversely, rosiglitazone or acetaminophen elicited similar phenotypic and transcriptomic responses between genotypes, demonstrating a lack of PPARα-dependence. In WT mice, HFPO-DA-mediated responses were similar to GW7647 but different from rosiglitazone or acetaminophen. Dose-dependent increases in liver weight, karyomegaly, and mitosis, as well as increased transcriptomic signaling related to PPARα activation and cell proliferation, were observed in HFPO-DA and GW7647-exposed WT mice. The consistent phenotypic and transcriptomic signaling patterns between HFPO-DA and GW7647 in WT mice, and the lack of changes in Ppara-null mice, provide further support that HFPO-DA-mediated early KEs in mouse liver are PPARα-dependent, occur via the rodent-specific PPARα MOA, and therefore are not appropriate for use in human health risk assessment.

Methods
Chemicals and Dosing Solutions

Ammonium perfluoro(2-methyl-3-oxahexanoate) (HFPO-DA; CASRN 62037-80-3; 95% purity) was purchased from Manchester Organics Ltd. (Runcorn, Cheshire, United Kingdom). GW7647 (CASRN 265129-71-3; ≥98% purity) was purchased from Cayman Chemical Company (Ann Arbor, MI). Rosiglitazone (CASRN 122320-73-4; 98.9% purity), acetaminophen (APAP; CASRN 103-90-2; ≥99% purity), dimethyl sulfoxide (DMSO; CASRN 67-68-5; >99.9% purity), and carboxylmethylcellulose sodium salt (CASRN 9004-32-4; Product No. C4888) were purchased from Sigma Aldrich (St. Louis, MO). HFPO-DA dosing solutions were prepared in deionized water; GW7647 and rosiglitazone dosing solutions were prepared in 5% DMSO suspended in 0.1% carboxymethylcellulose in deionized water. APAP dosing solutions were prepared in warm 0.9% saline.

Animal Husbandry and Exposure

The study was conducted at The Jackson Laboratory (JAX) (Sacramento, CA) under the supervision and approval of the Institutional Animal Care and Use Committee (IACUC study No. 20110-SR) and under Good Laboratory Practices (GLP)-like conditions. Male mice, aged 9 to 12 weeks, were acclimated for at least 2 days, and were housed in individually ventilated polysulfonate cages with HEPA filtered air at a density of up to 5 mice per cage at a temperature of 20°C - 26°C and a relative humidity of 30% - 70%, under an approximately 12-hour light-dark cycle. Male mice were chosen because of their greater sensitivity than female mice to HFPO-DA-mediated liver effects as observed in previous studies (Chappell et al. 2020; Heintz et al. 2022; Thompson et al. 2019). Mice were provided filtered tap water, acidified to a pH of 2.5-3.0, and standard rodent chow ad libitum.

Three strains of mice were used, including two wild-type (WT) strains, J:ARC(S) (stock no. 034608) and B6129SF2/J (stock no. 101045), and the Ppara-null strain B6;129S4-Pparatm1Gonz/J (stock no. 008154). J:ARC(S) mice, also referred to as JAX Swiss Outbred mice, are an outbred stock from Charles River Breeding Laboratories’ outbred CD-1(ICR) mouse colony.[1] Ppara-null mice were generated by deletion of 83 base pairs in exon 8 of the ligand binding domain of the mPPARa gene, resulting in a nonfunctional gene (Lee et al. 1995). The lack of functional PPARa protein activity was verified through the lack of activation of downstream PPARa target genes (Lee et al. 1995). No single genetic background strain of the Ppara-null mouse exists, but either the B6129SF2/J or C57BL/6J mouse strains are considered the most appropriate background strains.

This study consisted of four arms based on chemical treatment: HFPO-DA (chemical of interest), GW7647 (prototypical PPARa agonist), rosiglitazone (prototypical PPARg agonist), and APAP (prototypical cytotoxic agent). The exposure route, dose level, and duration for HFPO-DA and each of the three prototypical agonists for the relevant hypothesized MOAs were designed to capture the transcriptomic signature and relevant phenotypic effects of each chemical’s MOA. Male mice from each of the three strains were randomly allocated based on body weight to each dose group (vehicle control, low, medium, or high) within each study arm (n = 4/dose/strain/arm).

Dosing solutions of HFPO-DA, GW7647, and rosiglitazone were administered via oral gavage at a dose volume of 10 mL/kg body weight once per day for five days. The doses were as follows: 0, 3, 15, or 30 mg/kg-d HFPO-DA; 0, 5, 10, or 20 mg/kg-d GW7647; and 0, 5, 10, or 20 mg/kg-d rosiglitazone. Dose levels for HFPO-DA were selected based on previous short-term and subchronic OECD test guideline toxicity studies in mice (DuPont 2008a,b; DuPont 2010; Chappell et al. 2020). The highest dose, 30 mg/kg-d HFPO-DA, was selected based on previous 7- and 28-day toxicity studies in mice; the lowest dose, 3 mg/kg-d HFPO-DA, was selected based on minimal to no transcriptomic or histopathological changes following subchronic exposure to 0.1 or 0.5 mg/kg-d HFPO-DA. Dose levels for GW7647 and rosiglitazone were selected based on the mechanistic and physiological effects reported in previous mouse toxicity studies following short-term (5 to 9 days) exposure (Foreman et al. 2021; Tao et al. 2010; Edvardsson et al. 1999; Shen et al. 2007). The 5-day exposure duration used for HFPO-DA, GW7647, and rosiglitazone was selected to capture initial hepatic transcriptomic responses following chemical exposure. This duration is also consistent with the methods described in the EPA Transcriptomic Assessment Product (ETAP) (EPA 2024).

APAP was administered to 12-hour fasted mice as a single intraperitoneal (IP) injection at 0, 150, 300, or 600 mg/kg. Feed was restored immediately after APAP dosing, and mice were sacrificed 6 hours post-dose. Induction of acute APAP liver toxicity is both dose- and duration-sensitive (Bhushan et al. 2014); thus, in order to capture the transcriptomic signature for this mechanism, the experimental design employed was consistent with previous IP studies examining APAP-induced acute liver injury in mice (Bhushan et al. 2014; Kotulkar et al. 2023).

Clinical Observations & Tissue Isolation

Mice in the HFPO-DA, GW7647, and rosiglitazone study arms were monitored daily for clinical observations. Mice in the APAP arm were monitored closely after the single IP injection for severe adverse effects, including bleeding, severe lethargy, and blood in the urine: low dose – single observation post-dosing, medium dose – observed once/hour for up to 6 hours, and high dose – observations at 15, 30, 60, 90, and 120 minutes post-dose. If any severe adverse effects were observed in one or more animals within a dose group, all mice in that group were sacrificed in extremis at the same time. Mice that survived to scheduled study termination were euthanized via carbon dioxide asphyxiation either 24 hours after receiving the last dose via oral gavage (HFPO-DA, GW7647, and rosiglitazone) or 6 hours post-IP injection (APAP). Terminal body weights were taken prior to necropsy.

Immediately following the sacrifice, blood and livers were collected. Blood was collected via cardiocentesis and processed to serum and shipped to IDEXX (West Sacramento, CA) for measurement of clinical chemistry parameters including alkaline phosphatase (ALP), aspartate aminotransferase (AST), alanine transaminase (ALT), total cholesterol, triglycerides, high density lipoprotein (HDL), low density lipoprotein (LDL), and total bile acids. Livers were weighed and fixed in 10% formalin prior to shipment to Experimental Pathology Laboratories, Inc. (EPL; Durham, NC), where the livers were embedded in paraffin (FFPE), sectioned sequentially at 4-6 mM, and mounted on glass slides.

RNA Preparation and Sequencing

Mounted and unstained FFPE liver sections sequential to those used for histopathological examination were scraped from the slides and processed according to the TempO-Seq® protocol by BioSpyder Technologies, Inc. (Carlsbad, CA) to yield libraries of tagged (by sample) detector oligos ligated to RNA targets, as previously described (Yeakley et al. 2017). These libraries were sequenced using a NovaSeq X Ultra-High-Throughput Sequencing System (Illumina, San Diego, California).

Sequencing Data Processing and Assessment of Quality

Raw sequencing data for each sample (i.e., the FASTQ files) were analyzed using the TempO-Seq® data analysis pipeline, as previously described (Yeakley et al. 2017). For each study arm, the output of this pipeline was a table containing the number of sequenced reads for each TempO-Seq® probe per sample. Consistent with our previous studies, samples were reviewed for inclusion in the subsequent analyses based on two criteria: (1) overall sequencing depth ≥2 standard deviations below the mean across all samples within a study arm, and (2) number of sequenced probes ≥2 standard deviations below the mean across all samples within a study arm (Chappell et al. 2020; Heintz et al. 2022; Heintz et al. 2024a,b). Samples that did not meet one or both criteria were excluded from subsequent analyses, which were conducted using packages in the R software environment (version 4.4.1; cran.r-project.org/).

Differential Gene Expression Analyses

For each experimental arm, count data were normalized with the DESeq2 R package (v1.44.0) (Love et al. 2014) to account for sample-to-sample variation in sequencing depth across samples. Statistical comparisons between treatment and control groups from the same mouse strain were performed in DESeq2 to identify differentially expressed probes (DEPs) and determine fold change of DEPs. DEPs were defined as probes with a false discovery rate (FDR) of <10% based on p values adjusted for multiple comparisons using the Benjamini and Hochberg (BH) procedure (Love et al. 2014). In the TempO-Seq assay, some (but not all) genes are represented by multiple probes, such that the 30,146 mouse probes correspond to 21,398 mouse genes. Therefore, differentially expressed genes (DEGs) were identified from the corresponding DEPs (thus maintaining FDR <10%).

Identification of Pathway-Level Responses to Exposure

Gene set enrichment analysis was used to identify biological pathways associated with transcriptomic responses in the livers of mice exposed to HFPO-DA or one of the positive control chemicals. For genes with multiple corresponding probes, the probe with the highest sequencing count across all samples within a study arm was used in the gene set enrichment analysis. First, mouse gene identifiers were converted into human identifiers (when available) with the R package biomaRt (v2.60.1) based on the Ensembl genome database (http://uswest.ensembl.org/index.html). Next, the gene expression data with the human identifiers were queried for enrichment of gene sets within the canonical pathway (CP) subcollection (c2.cp.v2023.2). This collection includes gene sets from several pathway databases and is available through the Molecular Signatures Database (MSigDB; http://software.broadinstitute.org/gsea/msigdb/index.jsp). The hypergeometric test method for overrepresentation was used to identify enrichment of gene sets. DEGs (FDR <10%) for each treatment group within a study arm and mouse strain were tested for overrepresentation among gene sets in the CP subcollection using the Fisher combined probability test function within the Platform for Integrative Analysis of Omics data (PIANO) R package (v2.20.0) (Varemo et al. 2013). Gene sets with an FDR <5% were considered significantly enriched.

Benchmark Dose Analyses

BMDExpress software (v2.3; Phillips et al. 2019) was used to conduct dose–response modeling from the normalized gene expression data from DESeq2. Data were loaded into BMDExpress without transformation, using TempO-Seq probe IDs as the gene identifiers. Genes altered by chemical exposure for each mouse strain were identified with a Williams trend test (p value < 0.05; absolute log fold change ≥ 1.5 in at least one dose); corrections for multiple tests were not applied. Benchmark dose (BMD) analysis was conducted with the linear, power, hill, 2° and 3° polynomial, and exponential 2-5 models, assuming constant variance and a benchmark response (BMR) of 1 standard deviation. For each chemical and mouse strain, significant dose-responsive genes (DRGs) met the following criteria: a best BMD < 10-fold below the lowest tested dose, a best BMD ≤ the highest tested dose, and a winning model fit p value ≥ 0.1 (NTP 2018). The Reactome gene set collections available within the BMDExpress software were used for functional classification of significant DRGs and calculation of BMDs for enriched gene sets. Genes with BMD/BMD lower limit (BMDL) >20, BMD upper limit (BMDU)/BMD >20, and BMDU/BMDL >40 were removed from functional classification analyses (NTP 2018). There were no filters applied for minimum or maximum number of genes per gene set. Gene sets with a Fisher's Exact Right p-value < 0.05 were considered significantly enriched.

Analysis Match and Upstream Regulator Predictions

Fold change and FDR values (<10%) determined by DESeq2 for DEGs were analyzed using Qiagen Ingenuity Pathway Analysis (IPA, v. 01-22-01; Qiagen Bioinformatics, Redwood City, California) software. For genes with multiple corresponding probes, the probe with the highest sequencing count across all dose groups within a study arm was used. Gene expression profiles for HFPO-DA were compared and ranked against gene expression profiles of positive control chemical dose groups (excluding HFPO-DA) using IPA Analysis Match. The mid-dose group, 15 mg/kg-d HFPO-DA for each mouse strain, was selected for analysis to match comparisons to ensure that a maximum transcriptomic response was evaluated and avoid the potential for overt toxicity responses at the highest dose. In addition, gene expression profiles from the in vitro transcriptomic study in primary hepatocytes (Heintz et al. 2024a,b) were also included in the analysis.

Analysis matches were ranked using the overall experimental z-score. The overall z-score is a combined similarity score of the selected gene expression profile in question (i.e., mid-dose HFPO-DA for each mouse strain) compared to other transcriptomic profiles included in the analysis. The overall z-score is calculated using the average scores from IPA enrichment analysis of canonical pathway signatures, predicted upstream regulators, causal networks, and downstream effects. The top three analysis matches from positive control chemical dose groups from the 5-day in vivo study herein, as well as the single top analysis match from previous in vitro transcriptomic studies (Heintz et al. 2024a,b), were reported. For these top-ranked groups, activation/inhibition patterns of the top 20 predicted upstream regulators based on z-score (significance threshold >|2|) were examined.

[1] In 1981, mice from the CRL CD-1(ICR) colony were transferred to the Animal Resources Centre (ARC) in Australia. In 2020, this stock was transferred to JAX.