Supplementary data for: Comparison of transcriptomic profiles between HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in mouse, rat, and pooled human hepatocytes
Data files
May 02, 2024 version files 1.07 GB
-
README.md
-
Supplementary_File_S3a_DEPs_Mouse_CD-1.xlsx
-
Supplementary_File_S3b_DEPs_Mouse_B6129SF2J.xlsx
-
Supplementary_File_S3c_DEPs_Rat.xlsx
-
Supplementary_File_S3d_DEPs_Human.xlsx
-
Supplementary_File_S4_Hypergeometric_Gene_Set_Enrichment.csv
-
Supplementary_File_S5_Concentration-Responsive_Genes.csv
-
Supplementary_File_S6_Concentration-Responsive_Pathways.csv
-
Supplementary_File_S7_GSEA_Gene_Set_Enrichment_Human.xlsx
-
Supplementary_Files_S1___S2_Sequencing_Quality___Cytotoxicity_Results.xlsx
Abstract
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.
README: Supplementary data for: Comparison of transcriptomic profiles between HFPO-DA and prototypical PPARa, PPARg, and cytotoxic agents in mouse, rat, and pooled human hepatocytes
https://doi.org/10.5061/dryad.stqjq2c94
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: https://bmdexpress-2.readthedocs.io/en/feature-readthedocs/ 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: https://bmdexpress-2.readthedocs.io/en/feature-readthedocs/ 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 (https://www.ncbi.nlm.nih.gov/geo/) (GEO series accession number GSE248251).
Methods
Primary Hepatocyte Isolation and Culture
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 Supplementary File 1. 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 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, CD-1[1] and B6129SF2/J mouse hepatocytes were seeded at a density of 0.6 x 106 cells/mL, and rat and human hepatocytes were seeded at a density of 1.3 x 106 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.
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 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).
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/). 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.
Differential Gene Expression Analyses
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) < 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, 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.
Identification of Pathway-Level Responses to Treatment
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 <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 <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.
[1] Experiment was repeated for CD-1 mouse hepatocytes due to a sample mishandling error.