DOI: 10.5061/dryad.stqjq2c94
title: Supplementary data for: Comparison of transcriptomic profiles between HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in mouse, rat, and pooled human hepatocytes

description: Like many per- or polyfluorinated alkyl substances (PFAS), toxicity
 studies with HFPO-DA
 (ammonium,2,3,3,3-tetrafluoro-2-(heptafluoropropoxy)-propanoate), a
 short-chain PFAS used in the manufacture of some types of fluorinated
 polymers, indicate that the liver is the primary target of toxicity in
 rodents following oral exposure. Although the current weight of evidence
 supports the PPARa mode of action (MOA) for liver effects in
 HFPO-DA-exposed mice, alternate MOAs have also been hypothesized including
 PPARg or cytotoxicity. To further evaluate the MOA for HFPO-DA in rodent
 liver, transcriptomic analyses were conducted on samples from primary
 mouse, rat and pooled human hepatocytes treated for 12, 24 or 72 hours
 with various concentrations of HFPO-DA, or agonists of PPARa (GW7647),
 PPARg (rosiglitazone), or cytotoxic agents (i.e., acetaminophen or
 d-galactosamine). Concordance analyses of enriched pathways across
 chemicals within each species demonstrated greatest concordance between
 HFPO-DA and PPARa agonist GW7647-treated hepatocytes compared to the other
 chemicals evaluated. These findings were supported by benchmark
 concentration modeling and predicted upstream regulator results. In
 addition, transcriptomic analyses across species demonstrated a greater
 transcriptomic response in rodent hepatocytes treated with HFPO-DA or
 agonists of PPARa or PPARg, indicating rodent hepatocytes are more
 sensitive to HFPO-DA or PPARa/g agonist treatment. These results are
 consistent with previously published transcriptomic analyses and further
 support that liver effects in HFPO-DA-exposed rodents are mediated through
 rodent-specific PPARa signaling mechanisms as part of the MOA for PPARa
 activator-induced rodent hepatocarcinogenesis. Thus, effects observed in
 mouse liver are not appropriate endpoints for toxicity value development
 for HFPO-DA in human health risk assessment.
description: <em>Primary Hepatocyte Isolation and
 Culture</em> Mouse hepatocytes were isolated from
 the livers of 8 – 12 week-old male CD-1 mice (CD-1 IGS, strain code: 022)
 purchased from Charles River (Raleigh, NC) and 11 week-old male B6129SF2/J
 mice (stock #101045) purchased from The Jackson Laboratory (Bar Harbor,
 ME). Rat hepatocytes were isolated from the livers of 8 – 12 week-old male
 Sprague Dawley (SD, strain code: 001) rats purchased from Charles River.
 Human hepatocytes were isolated from a pool of 10 donor livers (5 male and
 5 female). Additional donor information is available in
 <strong>Supplementary File 1</strong>. Hepatocytes from all
 species/strains 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,
 CD-1[1] and B6129SF2/J mouse hepatocytes were seeded at
 a density of 0.6 x 10<sup>6</sup> cells/mL, and rat and human
 hepatocytes were seeded at a density of 1.3 x 10<sup>6</sup>
 cells/mL. After a 24 h adaptation period, hepatocytes from each
 species/strain were treated for 12, 24 or 72 h with supplemented 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 species/strain, 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 species/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/). The number of sequenced reads per TempO-Seq probe
 were extracted from FASTQ files generated from the sequencing experiments
 for each species/strain, 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 species/strain; 2) total number of sequenced probes lower than two
 standard deviations below the mean number of probes sequenced per sample
 from the same species/strain. Count data from all samples that were not
 excluded were used for further comparative analyses.
 <em>Differential Gene Expression
 Analyses</em> Data was normalized using the
 DESeq2 R package (v1.40.2) (Love et al. 2014) to account for
 sample-to-sample variation in sequencing depth within each species/strain.
 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
 species/strain 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, 20,922 rat, or 19,683 human genes as measured by
 30,146 mouse, 22,253 rat, or 22,533 human probes, respectively, 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 hepatocytes from each species/strain
 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. DEGs for each treatment group within a species/strain and timepoint (i.e., an FDR of &lt;10% as described above) were tested for overrepresentation among the 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.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. [1] Experiment was repeated for CD-1 mouse hepatocytes due to a sample mishandling error.
description: # Supplementary data for: Comparison of transcriptomic profiles between
 HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in mouse, rat,
 and pooled human hepatocytes
 To further evaluate the
 MOA for HFPO-DA in rodent liver, transcriptomic analyses were conducted on
 samples from primary mouse, rat and pooled human hepatocytes treated for
 12, 24 or 72 hours with various concentrations of HFPO-DA, or agonists of
 PPARa (GW7647), PPARg (rosiglitazone), or cytotoxic agents (i.e.,
 acetaminophen or d-galactosamine). ## 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 from the analyses described in the methods 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 File S1: Contains
 three sub-tables. Table 1 contains human liver donor information, Table 2
 lists the sequencing quality criteria across species and strains, and
 Table 3 lists the samples that failed sequencing quality criteria.
 Supplementary File S2: Provides hepatocyte cytotoxicity results.
 Supplementary File S3a-d: Provides differential gene expression analysis
 results of CD-1 mouse (File S3a), B6129SF2/J mouse (File S3b), rat (File
 S3c) and human (File S3d) hepatocytes treated with HFPO-DA, GW7647,
 rosiglitazone, acetaminophen, or d-galactosamine, as calculated using
 DESeq2. Within each spreadsheet, each tab provides differential gene
 expression results for one of the five chemicals tested, with the
 following variables listed as column headers in each spreadsheet tab: *
 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. More information on these variables and how they are
 estimated are provided in the DESeq2 vignette (Love et al. 2014).
 Supplementary File S4: Provides hypergeometric test results for gene set
 enrichment for each species/strain and chemical treatment group. The first
 four columns within the spreadsheet, i.e., Species_Strain, Chemical,
 Concentration, and Timepoint, provide information on the species and
 strain, chemical, chemical concentration, and exposure duration tested,
 respectively, for each canonical pathway/gene set (i.e., column named
 "Pathway") evaluated. Each subsequent column provides
 information on the significance level of the enrichment of each
 pathway/gene set, number of genes from the differential gene expression
 analysis that are part of a given pathway/gene set and overall estimated
 regulation direction of that pathway/gene set. See the piano vignette for
 a detailed description of each column header (Varemo et al. 2013).
 Supplementary File S5: Provides concentration-responsive gene results for
 each species/strain and chemical treatment group. The first three columns
 within the spreadsheet, i.e., Species_Strain, Chemical and Timepoint,
 provide information on the species and strain, chemical and exposure
 duration tested, respectively. The column named "Probe ID"
 provides the gene name and TempO-Seq probe ID number. Subsequent columns
 provide modeling results for each statistical model evaluated for each
 gene in BMDExpress. See Phillips et al. (2019) and BMDExpress
 documentation:
 for a
 detailed description of each column header. Supplementary File S6:
 Provides enriched gene sets across concentration-responsive genes for each
 species/strain and chemical treatment group. The first three columns
 within the spreadsheet, i.e., Species_Strain, Chemical and Timepoint,
 provide information on the species and strain, chemical and exposure
 duration tested, respectively. Subsequent columns provide information on
 the gene set, enrichment and modeling results for each Reactome gene set
 evaluated in BMDExpress. See Phillips et al. (2019) and BMDExpress
 documentation:
 for a
 detailed description of each column header. Supplementary File S7.
 Provides human hepatocyte gene set enrichment methods and results for each
 chemical treatment group using the GSEA (gene set enrichment analysis)
 test method. The first tab in the spreadsheet provides the methods for
 this analysis. The remaining tabs provide gene set enrichment results for
 each of the five chemicals tested, with results for one chemical per tab.
 For each tab, the first column within the spreadsheet, named
 Dose_Timepoint, indicates the concentration and treatment duration of the
 chemical tested. The second column, named Pathway, indicates the name of
 each canonical pathway/gene set evaluated. Each subsequent column provides
 information on the significance level of the enrichment of each
 pathway/gene set, number of genes from the differential gene expression
 analysis that are part of a given pathway/gene set and overall estimated
 regulation direction of that pathway/gene set. See the piano vignette for
 a detailed description of each column header (Varemo et al. 2013). For all
 columns except for the columns named GeneSymbols, GeneSymbols_up and
 GeneSymbols_down, an "NA" indicates that the value for that
 variable could not be estimated because none of the underlying genes in
 the gene set were determined to be differentially expressed. An
 "NA" present in one of the three columns that begin with
 "GeneSymbols" indicates that Gene Symbols were not available for
 the respective gene set. ## Sharing/Access information RNA sequencing data
 are publicly available at NCBI’s Gene Expression Omnibus
 () (GEO series accession number
 GSE248251).