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
Data files
Apr 21, 2025 version files 172.60 MB
-
README.md
20.52 KB
-
Supplementary_File_S4_Sequencing_Quality.csv
13.01 KB
-
Supplementary_File_S5a_HFPO-DA_Differentially_expressed_genes.csv
27.57 MB
-
Supplementary_File_S5b_GW7647_Differentially_expressed_genes.csv
26.97 MB
-
Supplementary_File_S5c_Rosiglitazone_Differentially_expressed_genes.csv
28.96 MB
-
Supplementary_File_S5d_APAP_Differentially_expressed_genes.csv
29.87 MB
-
Supplementary_File_S6_Gene_Set_Enrichment.csv
27 MB
-
Supplementary_File_S8._Gene_set_enrichment_of_dose-responsive_genes.csv
10.66 MB
-
Supplementray_File_S7_Dose-responsive_genes.csv
21.55 MB
Abstract
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.
https://doi.org/10.5061/dryad.hx3ffbgq5
To further inform the mode of action (MOA) of HFPO-DA (ammonium, 2,3,3,3-tetrafluoro-2-(heptafluoropropoxy)-propanoate) and evaluate the PPARa-dependence of HFPO-DA-mediated liver effects in vivo, 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 PPARa (GW7647) and PPARg (rosiglitazone), or cytotoxic agent (acetaminophen [APAP]). This Dryad dataset contains the supplementary data pertaining to the transcriptomic results of this study; supplementary data for phenotypic findings (i.e., clinical observations) are provided in Supplementary Files S1-S3 available with the main publication.
List of files
A list of each file included in this dataset:
- Supplementary File S4_Sequencing Quality.csv– Summary of sequencing quality data for HFPO-DA, GW7647, rosiglitazone, and APAP in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S5a_HFPO-DA_Differentially expressed genes.csv – Differential probe expression data for HFPO-DA in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S5b_GW7647_Differentially expressed genes.csv – Differential probe expression data for GW7647 in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S5c_Rosiglitazone_Differentially expressed genes.csv – Differential probe expression data for rosiglitazone in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S5d_APAP_Differentially expressed genes.csv – Differential probe expression data for APAP in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S6_Gene Set Enrichment.csv – Gene set enrichment results for HFPO-DA, GW7647, rosiglitazone, and APAP in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S7_Dose-responsive genes.csv – Gene dose-response analysis results for HFPO-DA, GW7647, and rosiglitazone in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains.
- Supplementary File S8. Gene set enrichment of dose-responsive genes.csv – Pathway dose-response analysis results (i.e., enriched gene sets) for HFPO-DA and GW7647 in WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice strains
Description of the data and file structure
Each file (Excel spreadsheet) contains supplementary data related to sequencing quality results or processed RNA sequencing data, including differentially expressed genes/probes, gene set enrichment, dose-responsive genes, and enriched gene sets among dose-responsive genes.
Below is a list of each supplementary file along with a description of its contents:
File: Supplementary_File_S4_Sequencing_Quality.csv
Description: Supplementary File S4 contains two sub-tables: Table 1 lists the mean sequencing quality criteria across mouse strains (J:ARC(S), B6129SF2/J, and Ppara-null) for each chemical (HFPO-DA, GW7647, rosiglitazone, and APAP). Table 2 lists the individual samples and indicates whether each sample failed sequencing quality criteria. The following variables are included in both Tables 1 (means) and 2 (by individual samples):
Variables
- Sequencing depth = total number of reads sequenced and aligned across all probes
- Number of probes = total number of probes sequenced
Files:
Supplementary File_S5a_HFPO-DA_Differentially_expressed_genes.csv
Supplementary File_S5b_GW7647_Differentially_expressed_genes.csv
Supplementary File_S5c_Rosiglitazone_Differentially_expressed_genes.csv
Supplementary File_S5d_APAP_Differentially_expressed_genes.csv
Description: Supplementary Files S5a-d provide differential gene expression analysis results from the livers of WT (J:ARC(S) and B6129SF2/J) and Ppara-null mice treated with HFPO-DA, GW7647, rosiglitazone, or APAP, as calculated using DESeq2. Each tab provides differential gene expression results for one of the four chemicals tested, with the following variables listed as column headers in each of the four spreadsheet tabs (as defined by Love et al. 2014):
Variables
- Chemical = chemical treatment
- Strain_Dose = mouse strain and dose of the chemical tested
- gene_source = gene name and TempO-Seq probe ID number
- baseMean = the average of the normalized gene count values, divided by size factors, taken over all samples
- log2FoldChange = the effect size estimate. This value indicates how much the gene's expression has changed between the treatment and control groups. This value is reported on a logarithmic scale to base 2. If the baseMean value = 0 (i.e., no counts for that gene), the log2 fold change will be equal to "NA", meaning a log2 fold change could not be estimated.
- lfcSE = the standard error estimate for the log2 fold change estimate. If the log2 fold change estimate = NA, the lfcSE estimate will also be equal to "NA", meaning that the standard error estimate of the log2 fold change could not be estimated.
- stat = the value of the Wald test statistic for the gene. Wald test statistic is the log2FoldChange divided by lfcSE, which is compared to a standard Normal distribution to generate a two-tailed P-value. If the log2 fold change and standard error estimate are equal to NA, the Wald test statistic with also be equal to "NA", meaning that the Wald test statistic could not be estimated.
- pvalue = P-value of the Wald test statistic for the gene. If the P-value = NA, either the Wald test statistic could not be estimated, or the gene contained an outlier count as identified using Cook's distance.
- padj = adjusted P-value for multiple testing for the gene. If the adjusted P-value = NA, either the P-value was also equal to "NA" or the gene did not pass the independent filtering threshold for the selection of genes for multiple test correction.
File: Supplementary_File_S6_Gene_Set_Enrichment.csv
Description: Supplementary File S6 provides hypergeometric test results for gene set enrichment for each mouse strain (J:ARC(S), B6129SF2/J, and Ppara-null) and chemical treatment group (HFPO-DA, GW7647, rosiglitazone, and APAP). The following variables are listed as column headers in this spreadsheet (as defined by Väremo et al. 2013):
Variables
- Chemical = chemical treatment
- Strain_dose = mouse strain and dose of the chemical tested
- Pathways = canonical pathway/gene set
- p-value – non-dir = non-directional (direction of differential expression omitted) p-value
- Adjusted p-value – non-dir = corrected non-directional p-value
- Significant (in gene set) – non-dir = number of differentially expressed genes in gene set
- Non-significant (in gene set) – non-dir = number of non-significant genes in gene set
- Significant (not in gene set) – non-dir = number of differentially expressed genes not in gene set
- Non-significant (not in gene set) – non-dir = number of non-significant genes not in gene set
- p-value – down = p-value for downregulation of gene set
- Adjusted p-value – down = corrected downregulation p-value
- Significant (in gene set) – down = number of differentially expressed downregulated genes in gene set
- Non-significant (in gene set) – down = number of non-differentially expressed genes in gene set (downregulation)
- Significant (not in gene set) – down = number of differentially expressed downregulated genes not in gene set
- Non-significant (not in gene set) – down = number of non-differentially expressed genes not in gene set (downregulation)
- p-value - up = p-value for upregulation of gene set
- Adjusted p-value - up = corrected upregulation p-value
- Significant (in gene set) - up = number of differentially expressed upregulated genes in gene set
- Non-significant (in gene set) - up = number of non-differentially expressed genes in gene set (upregulation)
- Significant (not in gene set) - up = number of differentially expressed upregulated genes not in gene set
- Non-significant (not in gene set) - up = number of non-differentially expressed genes not in gene set (upregulation)
- GeneSymbols_all = identification of differentially expressed genes (non-directional)
- GeneSymbols_up = identification of differentially expressed upregulated genes
- GeneSymbols_down = identification of differentially expressed downregulated genes
File: Supplementary_File_S7_Dose-responsive_genes.csv
Description: Supplementary File S7 provides dose-responsive gene results for each mouse strain (J:ARC(S), B6129SF2/J, and Ppara-null) and three chemical treatment groups (HFPO-DA, GW7647, rosiglitazone, and APAP). BMDExpress evaluates fit across the following models: Hill, Power, Linear, Polynomial 2, Polynomial 3, Exponential 2, Exponential 3, Exponential 4, and Exponential 5. The following variables are listed as column headers in this spreadsheet, noting that column headers regarding statistical modeling results are identical across models and therefore only listed once (as defined by Philips et al. 2019; BMDExpress, Undated; EPA, 2022):
Variables
-
Chemical = chemical treatment
-
Mouse strain = mouse strain
-
Probe ID = gene name and Tempo-Seq probe ID number
-
Entrez Gene IDs = universal gene identifier
-
Gene Symbols = gene symbols that are included in the probe
-
[Model] BMD = benchmark dose; estimated dose level associated with specified benchmark dose response
-
[Model] BMDL = benchmark dose lower bound of 95% confidence interval
-
[Model] BMDU = benchmark dose upper bound of 95% confidence interval
-
[Model] fitPValue = P-value for global goodness of fit of model to data
-
[Model] fitLogLikelihood = measure of model fit to compare between models
-
[Model] AIC = Akaike Information Criterion; estimates quality of model relative to the other models
-
[Model] adverseDirection = direction of the response (up or down)
-
[Model] BMD/BMDL = ratio of model’s BMD to BMDL for gene
-
Flagged [Model] = indicates model flagged for gene based on settings from BMD analysis (1 = flagged; 0 = not flagged)
-
[Model] Parameters Intercept, alpha, rho, control, slope, power, beta 0-18, v, n, k = various parameters for each different model
-
[Model] Execution Complete = indicates whether model calculation could be completed
-
Best Model = most appropriate model selected by the program for the dataset
- Parameters following “Best” indicate the value of the parameter from the best model
-
wAUC = weighted area under the curve
-
Prefilter P-Value = prefilter P-value from analysis of variance (ANOVA)
-
Prefilter Adjusted P-Value = prefilter P-value after Benjamini-Hochberg correction
-
Max Fold Change = maximum fold change in response across all dose levels
-
Max Fold Change Absolute Value = absolute value of the maximum fold change in response across all dose levels
-
NOTEL = no-observable transcriptional effect level
-
LOTEL = lowest observable transcriptional effect level
-
FC Dose Level n: fold change at dose n
Blank cells indicate parameters that could not be calculated for a model. “#NAME?” indicates calculations that could not be completed due to missing parameters that could not be calculated for a model.
File: Supplementary_File_S8._Gene_set_enrichment_of_dose-responsive_genes.csv
Description: Supplementary File S8 provides enriched gene sets across dose-responsive genes for each mouse strain (J:ARC(S), B6129SF2/J, and PPpara-null) and three chemical treatment groups (HFPO-DA, GW7647, rosiglitazone, and APAP). The following variables (or sets of variables) are listed as column headers in this spreadsheet that provide information on the gene set, enrichment, and modeling results for each Reactome gene set evaluated in BMDExpress (as defined by Philips et al. 2019; BMDExpress, Undated):
Variables
-
Chemical = chemical treatment
-
Mouse strain = mouse strain
-
GO/Pathway/Gene Set/Gene ID = category identifier
-
GO/Pathway/Gene Set/Gene Name = GO category, reactome pathway, or user-defined gene set names
-
All Genes (Expression Data) = total number of genes loaded that were assigned to the pathway/category
-
All Genes (Platform) = total number of genes assigned to the pathway/category from the platform utilized (i.e., Reactome)
-
Input Genes = total number of genes assigned to the pathway/category that passed ANOVA pre-filtering and BMD analyses
-
Genes with BMD <= Highest Dose = number of genes with BMD values less than the highest dose
-
Genes with BMD p-Value >= = number of genes with fit P-values
-
Genes with BMD/BMDL <= = number of genes with BMD to BMDL ratios less than 20
-
Genes with BMDU/BMD <= = number of genes with BMDU to BMD ratios less than 20
-
Genes with BMDU/BMDL <= = number of genes with BMDU to BMDL ratios less than 40
-
Genes with BMD <= N-Fold Lowest Positive Dose = number of genes with BMD values greater than 10-fold below the lowest positive dose
-
Genes with Prefilter P-Value <= = total number of genes that passed P-value pre-filter (0.05)
-
Genes with Prefilter Adjusted P-Value <= = number of genes that passed adjusted P-value filter (0.5)
-
Genes That Passed All Filters = total number of genes that passed all prefilters
-
Fisher's A Parameter = number of genes related to category term with a BMD
-
Fisher’s B Parameter = number of genes not related to category term with a BMD
-
Fisher’s C Parameter = number of genes related to category term without a BMD
-
Fisher’s D Parameter = number of genes not related to category term without a BMD
-
Fisher's Exact Left P-Value = P-value of enrichment analysis of gene sets that indicates significant underrepresentation of the gene set’s responsive genes
-
Fisher's Exact Right P-Value = P-value of enrichment analysis of gene sets that indicates significant overrepresentation of the gene set's responsive genes
-
Fisher's Exact Two-Tailed = P-value of enrichment analysis of gene sets from the two-tailed test that can indicate either significant over- or underrepresentation of gene sets
-
Percentage = percent (0-1) of total number of genes on array input into BMD analysis that had reportable BMDs based on filter criteria in Functional Classification setup
-
Entrez Gene IDs = universal gene identifiers from probes/probe sets used as input in BMD analysis
-
Gene Symbols = unique gene symbols corresponding to Gene IDs input into BMD analysis
-
Probe IDs = probe identifiers input into BMD analysis
-
Genes with Conflicting Probesets = Entrez Gene IDs with conflicting probe sets or probes (with correlation coefficient)
-
BMD/BMDL/BMDU Mean = mean BMD/BMDL/BMDU for genes in category that had a calculated BMD/BMDL/BMDU
-
BMD/BMDL/BMDU Median = median BMD/BMDL/BMDU for genes in category that had a calculated BMD/BMDL/BMDU
-
BMD/BMDL/BMDU Minimum = minimum BMD/BMDL/BMDU for genes in category that had a calculated BMD/BMDL/BMDU BMD/BMDL/BMDU
-
BMD/BMDL/BMDU Standard Deviation = standard deviation of BMD/BMDL/BMDU for genes in the category that had a calculated BMD/BMDL/BMDU
-
BMD/BMDL/BMDU wMean = weighted mean BMD/BMDL/BMDU for genes in category that had a calculated BMD/BMDL/BMDU
-
BMD/BMDL/BMDU wSD = weighted standard deviation BMD/BMDL/BMDU for genes in category that had a calculated BMD/BMDL/BMDU
-
5th/10th Percentile Index = Nth gene number representing 5th/10th percentile for all genes in category
-
BMD at 5th/10th Percentile of Total Genes = BMD at 5th/10th percentile for all genes in Reactome category, including genes without significant dose-response
-
BMD/BMDL/BMDU List = list of BMD/BMDL/BMDU values for genes in this set
-
Probes with Adverse Direction Up = number of probes/probe sets in GO category with adverse direction up (increased)
-
Probes with Adverse Direction Down = number of probes/probe sets in GO category with adverse direction down (decreased)
-
Genes with Adverse Direction Up = number of genes in set with adverse direction up
-
Genes Up List = Entrez Gene IDs with adverse direction up
-
Genes Up Probes List = Probes/probe sets with adverse direction up
-
Genes Up BMD/BMDL/BMDU Mean = mean of BMD/BMDL/BMDU with adverse direction up
-
Genes Up BMD/BMDL/BMDU Median = median of BMD/BMDL/BMDU with adverse direction up
-
Genes Up BMD/BMDL/BMDU SD = standard deviation of BMD/BMDL/BMDU with adverse direction up
-
BMD/BMDL/BMDU list (up) = list of BMD/BMDL/BMDU values for genes with adverse direction up
-
Genes with Adverse Direction Down = number of genes in set with adverse direction down
-
Genes Down List = Entrez Gene IDs with adverse direction down
-
Genes Down Probes List = Probes/probe sets with adverse direction down
-
Genes Down BMD/BMDL/BMDU Mean = mean of BMD/BMDL/BMDU with adverse direction down
-
Genes Down BMD/BMDL/BMDU Median = median of BMD/BMDL/BMDU with adverse direction down
-
Genes Down BMD/BMDL/BMDU SD = standard deviation of BMD/BMDL/BMDU with adverse direction down = list of BMD/BMDL/BMDU values for genes with adverse direction down
-
BMD/BMDL/BMDU list (down) = list of BMD/BMDL/BMDU values for genes with adverse direction down
-
Genes with Adverse Conflict Count = number of genes with probes in category that have conflicting adverse directions
-
Genes Conflict List = Entrez Gene IDs with probes not correlated above the designated threshold in the Functional Classification setup
-
Genes Conflict Probes List = probes mapping to the same gene not correlated above the designated threshold in the Functional Classification setup
-
BMD/BMDL/BMDU list (Conflict) = BMD/BMDL/BMDU values from conflicted probes
-
Model Counts = model counts denoting model, total number of genes in that model, and percentage of total model fits the model accounts for
-
Mean/Total/Min/Max/Standard Deviation/Median Fold Change = mean, summed, minimum, maximum, standard deviation, and median of absolute value of fold change values for genes in category
-
Lower bound of the 95% confidence interval - BMD/BMDL/BMDU = lower bound of 95% confidence interval of BMD/BMDL/BMDU values (assumes normal distribution)
-
Upper bound of the 95% confidence interval - BMD/BMDL/BMDU = upper bound of 95% confidence interval of BMD/BMDL/BMDU values (assumes normal distribution)
-
Overall Direction = considered “Up” if >60% of genes in category were up-regulated; considered “Down” if >60% of genes in category were down-regulated; considered “Conflict” if neither criterion were met
-
Percent Genes With Overall Direction Up = fraction of genes up-regulated based on the adverse direction upwards of genes
-
Percent Genes With Overall Direction Down = fraction of genes down-regulated based on the adverse direction downwards of genes
-
Percent Genes With Overall Direction Conflict = fraction of genes with probes that are both up and down-regulated
Blank cells indicate parameters that were not applicable or could not be calculated.
References
BMDExpress. (Undated). BMDExpress-2 Documentation. Available from https://bmdexpress-2.readthedocs.io/en/feature-readthedocs/. Accessed 2AUG24.
EPA. (2022). Benchmarck Dose Software. Version 3.3. User Guide. EPA/600/R-21/245. July 2022.
Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 15, 1-21.
Phillips, J. R., Svoboda, D. L., Tandon, A., Patel, S., Sedykh, A., Mav, D., ... & Auerbach, S. S. (2019). BMDExpress 2: enhanced transcriptomic dose-response analysis workflow. Bioinformatics, 35(10), 1780-1782.
Väremo, L., Nielsen, J., & Nookaew, I. (2013). Enriching the gene set analysis of genome-wide data by incorporating directionality of gene expression and combining statistical hypotheses and methods. Nucleic acids research, 41(8), 4378-4391.
Code/software
These files can be reviewed with Microsoft Excel or similar.
Access information
Other publicly accessible locations of the data:
- RNA sequencing data are publicly available at NCBI’s Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/) (GEO series accession number GSE282998).
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.