DOI: 10.5061/dryad.pc866t1wp
title: Supplementary data for: Comparison of transcriptomic profiles between HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in wild-type and PPARa knockout mouse hepatocytes

description: Recent in vitro transcriptomic analyses for the short-chain
 polyfluoroalkyl substance (PFAS), HFPO-DA (ammonium,
 2,3,3,3-tetrafluoro-2-(heptafluoropropoxy)-propanoate), support
 conclusions from in vivo data that HFPO-DA-mediated liver effects in mice
 are part of the early key events of the peroxisome proliferator-activated
 receptor alpha (PPARa)activator-induced rodent hepatocarcinogenesis mode
 of action (MOA). Transcriptomic responses in HFPO-DA-treated rodent
 hepatocytes have high concordance with those treated with a PPARa agonist
 and lack concordance with those treated with PPARg agonists or cytotoxic
 agents. To elucidate whether HFPO-DA-mediated transcriptomic responses in
 mouse liver are PPARa-dependent, additional transcriptomic analyses were
 conducted on samples from primary PPARaknockout (KO) and wild-type (WT)
 mouse hepatocytes exposed for 12, 24 or 72 hours with various
 concentrations of HFPO-DA, or well-established agonists of PPARa (GW7647)
 and PPARg (rosiglitazone), or cytotoxic agents (acetaminophen or
 d-galactosamine). Pathway and predicted upstream regulator-level responses
 were highly concordant between HFPO-DA and GW7647 in WT hepatocytes. A
 similar pattern was observed in PPARa KO hepatocytes, albeit with a
 distinct temporal and concentration-dependent delay potentially mediated
 by compensatory responses. This delay was not observed in PPARa KO
 hepatocytes exposed to rosiglitazone, acetaminophen, d-galactosamine. The
 similarity in transcriptomic signaling between HFPO-DA and GW7647 in both
 the presence and absence of PPARa in vitro indicates these compounds share
 a common mechanism of action and supports the PPARa-dependence of
 HFPO-DA-mediated effects in mouse liver.
description: <em>Primary Hepatocyte Isolation and
 Culture</em> Mouse hepatocytes were isolated from
 the livers of 11 week-old male B6129SF2/J mice (stock #101045) and 11
 week-old male PPARa-null mice
 (B6;129S4-Ppara<sup>tm1Gonz</sup>/J, stock #008154) purchased
 from The Jackson Laboratory (Bar Harbor, ME). As described in Lee et al.
 (1995), PPARa-null mice were generated by a targeted disruption of the
 ligand-binding domain (i.e., deletion of 83 base pairs in exon 8, see Lee
 et al. 1995 for details) of the mouse PPARa (mPPARa) gene, rendering the
 mPPARa gene nonfunctional. Although nonfunctional, mPPARa RNA was detected
 in these mice at very low expression levels; however, mice completely lack
 expression of mPPARa protein as measured by Western blotting, and lack
 functional protein activity as indicated by the inability to activate
 downstream PPARa target genes (Lee et al. 1995). PPARa-null mice are
 considered constitutive knockout mice (i.e., PPARa is nonfunctional in the
 entire animal), thus primary hepatocytes from PPARa-null mice will be
 referred to as PPARa knockout (KO) hepatocytes. B6129SF2/J mice were used
 as the genetic background strain for PPARa KO mice and are referred to as
 wild-type (WT) hepatocytes. Hepatocytes from both mouse
 genotypes were isolated using a two-step enzymatic digestion of liver
 tissue as described in Mudra and Parkinson (2001). Hepatocyte viability
 was determined by trypan blue (0.04%; Millipore Sigma, St. Louis, MO)
 exclusion and was ≥ 79%.  Primary hepatocytes were
 plated in a collagen-sandwich configuration on 48-well plates. Hepatocytes
 were maintained in Modified Eagle's medium, Dr. Chee’s modification
 (MCM; Gibco, Grand Island, NY) supplemented with 25.9 mM
 NaHCO<sub>3</sub> (Sigma Aldrich), 953.56 μM L-arginine (Sigma
 Aldrich), 3.95 mM L-glutamine (Sigma Aldrich), 40.51 μM thymidine (Sigma
 Aldrich), 0.988% (v/v) ITS+ (6.18 μg/mL insulin, 6.18 μg/mL transferrin,
 6.18 ng/mL selenium, 5.29 μg/mL linoleic acid, 1.24 mg/mL BSA; BD
 Biosciences, Bedford, MA), 98.8 μg/mL primocin (InvivoGen, San Diego, CA),
 and  0.099 μM dexamethasone (Sigma Aldrich). Cell cultures were incubated
 in a humidified culture chamber (37 ± 2 °C at 95% relative humidity, 95/5%
 air/CO<sub>2</sub>). <em>Hepatocyte
 Treatments</em> Using 48-well plates, WT and
 PPARa KO mouse hepatocytes were seeded at densities of 0.6 ´
 10<sup>6</sup> cells/mL and 0.5 ´
 10<sup>6</sup>cells/mL, respectively. After a 24 h adaptation
 period, hepatocytes from each genotype were treated for 12, 24 or 72 h
 with supplemented Modified Eagle's medium (MCM) media containing
 solvent control in the presence or absence of HFPO-DA (0.1, 5, 50 or 500
 μM) or one of the following positive controls: GW7647 (0.01, 0.1, 1 or 10
 μM), rosiglitazone (0.01, 0.1, 1 or 10 μM), acetaminophen (0.3, 1, 3 or 10
 mM) or d-galactosamine (0.3, 1, 3 or 10 mM). Deionized water (1%) served
 as the solvent control for HFPO-DA and dimethylsulfoxide (DMSO, 0.1%; cell
 culture grade; Sigma Aldrich) served as the solvent control for the
 remaining test chemicals. Treatment solutions were replaced every 24 h.
 For each genotype, treatment groups were performed in triplicate wells for
 12 and 72 h treatment durations, and quadruplicate wells for the 24 h
 treatment duration. At 24, 48 and 72 h following
 treatment, hepatocyte cultures were visualized with a Nikon TMS Microscope
 (Nikon Corporation) or Accu-Scope 3032 Inverted Microscope (Accu Scope
 Inc.), and representative hepatocytes from each species/strain treatment
 group were photographed with a PAXcam5 digital camera (MIS Inc.) to
 document morphological integrity.
 <em>Cytotoxicity Assay</em> The
 release of lactate dehydrogenase (LDH) into culture medium is an indicator
 of loss of cell membrane integrity and was used to estimate cytotoxicity.
 LDH release was measured using a commercial kit (Roche Diagnostics GmbH,
 Mannheim, Germany) according to the manufacturer’s directions. Briefly, at
 12, 24 and 72 h, medium samples were collected from hepatocytes treated
 with solvent controls (DMSO or deionized water) only, supplemented MCM
 only (negative control for LDH assay) or test chemicals. Three untreated
 hepatocyte samples per strain and timepoint were treated with 1% Triton
 X‑100 solution (positive control for LDH assay) and incubated for a period
 of 30 to 120 min. Aliquots of each medium sample were transferred to a
 96‑well plate and mixed with aliquots of the LDH assay working solution to
 begin the reaction. LDH activity for each medium sample was measured using
 a spectrophotometer at 490 nm (BioTek Instruments, Inc.). Cytotoxicity in
 a treatment group was determined based on measurements of percent LDH
 release ³ 25% in combination with changes in hepatocyte morphology
 indicative of cytotoxicity. A preliminary cytotoxicity assay was performed
 to select treatment concentrations used in the present study (data not
 shown). <em>RNA Preparation and
 Sequencing</em> Following treatment with HFPO-DA
 or positive controls (i.e., GW7647, rosiglitazone, acetaminophen or
 d-galactosamine) for 12, 24 or 72 h, primary hepatocytes for each
 species/strain were lysed using TempO-Seq® Enhanced Lysis Buffer and
 processed according to the TempO-Seq® protocol by BioSypder Technologies
 (Carlsbad, California), as previously described (Yeakley et al. 2017).
 Resultant DNA libraries were sequenced using a HiSeq 2500
 Ultra-High-Throughput Sequencing System (Illumina, San Diego,
 California). <em>Data Processing and
 Analysis</em> Sequencing data were analyzed using
 packages in the R software environment, version 4.3.1
 (cran.r-project.org/). FASTQ files generated from the sequencing
 experiments for each mouse genotype provided the number of sequenced reads
 per TempO-Seq probe, with each probe representing a gene-specific
 sequence. Samples were excluded from the downstream analyses if either or
 both of the following exclusion criteria were met: 1) overall sequencing
 depth (total number of reads across all probes) lower than two standard
 deviations below the mean sequencing depth across all samples from the
 same genotype; 2) total number of sequenced probes lower than two standard
 deviations below the mean number of probes sequenced per sample from the
 same genotype. Count data from all samples that were not excluded were
 used for further comparative analyses.
 <em>Differential Gene Expression
 Analyses</em> The DESeq2 R package (v1.40.2)
 (Love et al. 2014) was used to normalize data and account for
 sample-to-sample variation in sequencing depth within each mouse genotype.
 Fold-change and differentially expressed probes (DEPs) associated with
 chemical treatment were determined within DESeq2 by conducting statistical
 comparisons between treatment groups and controls from the same mouse
 genotype and treatment duration. DEPs were defined as those with a false
 discovery rate (FDR) &lt; 10%, based on <em>p </em>values
 adjusted for multiple testing using the Benjamini and Hochberg (BH)
 procedure (Love et al. 2014); differentially expressed genes (DEGs) were
 identified from respective DEPs, as some genes (but not all) are
 represented by multiple probes in the TempO-Seq assay. The expression
 levels of 21,398 mouse genes as measured by 30,146 mouse probes were
 reported from the TempO-Seq assay for each sample.
 <em>Identification of Pathway-Level Responses to
 Treatment</em> Biological pathways associated
 with transcriptomic responses in mouse hepatocytes following treatment
 with HFPO-DA or positive controls were identified by gene set enrichment
 analysis. For genes for which multiple probes were used to measure
 expression, the probe with the highest sequencing count across all samples
 was used in the pathway analyses. Mouse or rat gene identifiers were
 converted into human identifiers using the R package biomaRt (v2.56.1)
 based on the Ensembl genome database (http://uswest.ensembl.org/index.html). Gene expression data were then queried for enrichment of gene sets within the canonical pathway (CP) subcollection (c2.cp.v2022.1) available through the Molecular Signatures Database (MSigDB; http://software.broadinstitute.org/gsea/msigdb/index.jsp), which includes gene sets from several pathway databases.  Enrichment of sets of genes (i.e., the constituents of a molecular signaling pathway) was evaluated using the hypergeometric test method for overrepresentation. Significant DEGs (i.e., an FDR of &lt;10% as described above) for each treatment group, timepoint, and mouse genotype were tested for overrepresentation among the gene sets in the canonical pathway subcollection using the Fisher combined probability test function within the Platform for Integrative Analysis of Omics data (PIANO) R package (v2.16.0) (Varemo et al. 2013). Gene sets with an FDR &lt;5% were considered significantly enriched. <em>Benchmark Concentration Analyses</em> Concentration–response modeling was conducted using the BMDExpress software (v2.3) (Phillips et al. 2019). Normalized expression data for all samples as generated using DESeq2 were loaded into BMDExpress without transformation, using probe IDs from the TempO-Seq experiment as gene identifiers. A Williams trend test (with <em>p</em> value cutoff = 0.05) was used to identify genes altered by chemical treatment for each species/strain and timepoint. No fold-change filters or correction for multiple tests were applied. Benchmark concentration (BMC) analysis was conducted using the following models: linear, power, hill, 2° and 3° polynomial, and exponential models 2–5. The models were run assuming constant variance and a benchmark response (BMR) of 1 standard deviation. Concentration-responsive genes with a best BMC &gt;10-fold below the lowest concentration or a best BMC &gt; the highest concentration were removed. Functional classification was conducted using the Reactome gene set collections available within the BMDExpress software, based on significantly concentration-responsive genes (i.e., genes with a winning model fit <em>p</em> value ≥0.1), and removing genes according to the default parameters as follows: genes with BMC/BMCL &gt;20, BMCU/BMC &gt;20, and BMCU/BMCL &gt;40. No filters for minimum or maximum number of genes per gene set were applied. Benchmark concentrations for the gene sets were also calculated.
description: # Supplementary data for: **Comparison of transcriptomic profiles between
 HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in wild-type
 and PPARa knockout mouse hepatocytes**
 [https://doi.org/10.5061/dryad.pc866t1wp](https://doi.org/10.5061/dryad.pc866t1wp) To elucidate whether HFPO-DA-mediated transcriptomic responses in mouse liver are PPARa-dependent, transcriptomic analyses were conducted on samples from primary PPARa knockout (KO) and wild-type (WT) mouse hepatocytes exposed for 12, 24 or 72 hours with various concentrations of HFPO-DA, or established agonists of PPARa (GW7647) and PPARg (rosiglitazone), or cytotoxic agents (acetaminophen or d-galactosamine). ## **List of files** A list of each file included in this dataset: ·       **Supplementary Files S1 & S2_Sequencing Quality & Cytotoxicity Results.xlsx** – Summary of sequencing quality data and observations of cytotoxicity. ·       **Supplementary File S3a_DEPs_Mouse_B6129SF2J.xlsx** – Differential probe expression data for WT mouse. ·       **Supplementary File S3b_DEPs_Mouse_PPARa-KO.xlsx** – Differential probe expression data for PPARa-KO mouse. ·       **Supplementary File S4_Hypergeometric Gene Set Enrichment.csv** – Gene set enrichment results for both mouse strains. ·       **Supplementary File S5_Concentration-Responsive Genes.csv** – Gene dose-response analysis results for both mouse strains. ·       **Supplementary File S6_Concentration-Responsive Pathways.csv** – Pathway dose-response analysis results (i.e., enriched gene sets) for both mouse strains. ## Description of the data and file structure Each file (Excel spreadsheet) contains supplementary data related to sequencing quality results, cytotoxicity results, or processed RNA sequencing data including differentially expressed genes/probes, gene set enrichment, concentration-responsive genes, and enriched gene sets among concentration-responsive genes. Below is a list of each supplementary file along with a description of its contents: ### Supplementary Files S1 & S2\_Sequencing Quality & Cytotoxicity Results.xlsx Supplementary File S1 contains two sub-tables: Table 1 lists the mean sequencing quality criteria across mouse strains (WT and PPARa-KO). Table 2 lists the individual samples that failed sequencing quality criteria. The following variables are included in both Tables 1 (means) and 2 (by individual samples): ·      Sequencing depth = total number of reads sequenced and aligned across all probes ·      Number of probes = total number of probes sequenced Supplementary File S2: Provides hepatocyte cytotoxicity results across strains, chemical doses, and timepoints. ### Supplementary File S3a\_DEPs\_Mouse\_B6129SF2J.xlsx Supplementary File 3a provides differential gene expression analysis results of WT (B6129SF2/J) mouse hepatocytes treated with HFPO-DA, GW7647, rosiglitazone, acetaminophen, or d-galactosamine, as calculated using DESeq2. Each tab provides differential gene expression results for one of the five chemicals tested, with the following variables listed as column headers in each of the five spreadsheet tabs (as defined by Love et al. 2014): * Dose_Timepoint = concentration and treatment duration 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. ### Supplementary File S3b\_DEPs\_Mouse\_PPARa-KO.xls Supplementary File 3b provides differential gene expression analysis results of PPARa-KO mouse hepatocytes treated with HFPO-DA, GW7647, rosiglitazone, acetaminophen, or d-galactosamine, as calculated using DESeq2. Each tab provides differential gene expression results for one of the five chemicals tested, with the same variables listed as column headers in each of the five spreadsheet tabs as were listed in **Supplementary File 3a** (Supplementary File S3a_DEPs_Mouse_B6129SF2J.xlsx). ### Supplementary File S4\_Hypergeometric Gene Set Enrichment.csv Supplementary File S4 provides hypergeometric test results for gene set enrichment for each mouse strain and chemical treatment group. The following variables are listed as column headers in this spreadsheet (as defined by Varemo et al. 2013): ·      Species_Strain = species and mouse strain ·      Chemical = chemical treatment ·      Concentration = treatment concentration ·      Timepoint = exposure duration ·      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 = identify of differentially expressed genes (non-directional) ·      GeneSymbols_up = identify of differentially expressed upregulated genes ·      GeneSymbols_down = identify of differentially expressed downregulated genes ### Supplementary File S5\_Concentration-Responsive Genes.csv Supplementary File S5 provides concentration-responsive gene results for each mouse strain and chemical treatment group. 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 identical across models and therefore only listed once (as defined by Philips et al. 2019; BMDExpress, Undated; EPA, 2022): ·      Species_Strain = species and mouse strain ·      Chemical = chemical treatment ·      Timepoint = exposure duration ·      Probe ID = gene name and Tempo-Seq probe ID number ·      Entrez Gene IDs = universal gene identifier ·      Genes 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 program for dataset. Parameters listed in column headers following “Best Model” 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?” indicate calculations that could not be completed due to missing parameters that could not be calculated for a model. ### Supplementary File S6\_Concentration-Responsive Pathways Supplementary File S6 provides enriched gene sets across concentration-responsive genes for each mouse strain and chemical treatment group. 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): ·      Species_Strain = species and mouse strain ·      Chemical = chemical treatment ·      Timepoint = exposure duration ·      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 gene set’s responsive genes ·      Fisher's Exact Right P-Value = P-value of enrichment analysis of gene sets that indicates significant overrepresentation of gene set’s responsive genes ·      Fisher's Exact Two Tail = P-value of enrichment analysis of gene sets that from 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 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 designated threshold in Functional Classification setup ·      Genes Conflict Probes List = probes mapping to same gene not correlated above designated threshold in 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 mode, 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 criteria were met ·      Percent Genes With Overall Direction Up = fraction of genes up-regulated based on adverse direction up of genes ·      Percent Genes With Overall Direction Down = fraction of genes down-regulated based on adverse direction down 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. ## Sharing/Access information 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 GSE248251). ## 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.