DOI: 10.5061/dryad.hx3ffbgq5
title: 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

description: 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.
description: <em>Chemicals and Dosing Solutions</em>
 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;
 &gt;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.
 <em>Animal Husbandry and Exposure</em>
 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
 <em>ad libitum</em>. 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-Ppara<sup>tm1Gonz</sup>/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).
 <em>Clinical Observations & Tissue
 Isolation</em> 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 <em>in extremis</em> 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. <em>RNA Preparation and
 Sequencing</em> 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). <em>Sequencing Data Processing and
 Assessment of Quality</em> 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/). <em>Differential Gene
 Expression Analyses</em> 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
 &lt;10% based on <em>p</em> 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
 &lt;10%). <em>Identification of Pathway-Level
 Responses to Exposure</em> 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 &lt;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 &lt;5% were considered significantly enriched. <em>Benchmark Dose Analyses</em> 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 (<em>p</em> value &lt; 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 &lt; 10-fold below the lowest tested dose, a best BMD ≤ the highest tested dose, and a winning model fit <em>p</em> 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) &gt;20, BMD upper limit (BMDU)/BMD &gt;20, and BMDU/BMDL &gt;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 <em>p</em>-value &lt; 0.05 were considered significantly enriched. <em>Analysis Match and Upstream Regulator Predictions</em> Fold change and FDR values (&lt;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 <em>in vitro</em> 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 <em>in vivo</em> study herein, as well as the single top analysis match from previous <em>in vitro</em> 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 &gt;|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.
description: # 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
 [https://doi.org/10.5061/dryad.hx3ffbgq5](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 &lt;= Highest Dose = number of genes with BMD values less than the highest dose * Genes with BMD p-Value &gt;= = number of genes with fit P-values * Genes with BMD/BMDL &lt;= = number of genes with BMD to BMDL ratios less than 20 * Genes with BMDU/BMD &lt;= = number of genes with BMDU to BMD ratios less than 20 * Genes with BMDU/BMDL &lt;= = number of genes with BMDU to BMDL ratios less than 40 * Genes with BMD &lt;= 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 &lt;= = total number of genes that passed P-value pre-filter (0.05) * Genes with Prefilter Adjusted P-Value &lt;= = 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 &gt;60% of genes in category were up-regulated; considered “Down” if &gt;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/](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/](https://www.ncbi.nlm.nih.gov/geo/)) (GEO series accession number GSE282998).