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
Data files
Aug 06, 2024 version files 575.97 MB
-
README.md
-
Supplementary_File_S3a_DEPs_Mouse_B6129SF2J.xlsx
-
Supplementary_File_S3b_DEPs_Mouse_PPARa-KO.xlsx
-
Supplementary_File_S4_Hypergeometric_Gene_Set_Enrichment.csv
-
Supplementary_File_S5_Concentration-Responsive_Genes.csv
-
Supplementary_File_S6_Concentration-Responsive_Pathways.csv
-
Supplementary_Files_S1___S2_Sequencing_Quality___Cytotoxicity_Results.xlsx
Abstract
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.
README: 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
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 <= 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 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 >60% of genes in category were up-regulated; considered “Down” if >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/) (GEO series accession number GSE248251).
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.
Methods
Primary Hepatocyte Isolation and Culture
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-Pparatm1Gonz/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 NaHCO3 (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/CO2).
Hepatocyte Treatments
Using 48-well plates, WT and PPARa KO mouse hepatocytes were seeded at densities of 0.6 ´ 106 cells/mL and 0.5 ´ 106cells/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.
Cytotoxicity Assay
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).
RNA Preparation and Sequencing
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).
Data Processing and Analysis
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.
Differential Gene Expression Analyses
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) < 10%, based on p 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.
Identification of Pathway-Level Responses to Treatment
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 <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 <5% were considered significantly enriched.
Benchmark Concentration Analyses
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 p 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 >10-fold below the lowest concentration or a best BMC > 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 p value ≥0.1), and removing genes according to the default parameters as follows: genes with BMC/BMCL >20, BMCU/BMC >20, and BMCU/BMCL >40. No filters for minimum or maximum number of genes per gene set were applied. Benchmark concentrations for the gene sets were also calculated.