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

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.

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.