Per- and polyfluoroalkyl substances (PFAS) are a wide range of chemicals that are used in a variety of consumer and industrial products leading to direct human exposure. Many PFAS are chemically non-reactive and persistent in the environment, resulting in additional exposure from water, soil, and dietary intake. While some PFAS have documented negative health effects, data on simultaneous exposures to multiple PFAS (PFAS mixtures) are inadequate for making informed decisions for risk assessment. The current study leverages data from previous work in our group using Templated Oligo-Sequencing (TempO-Seq­™) for high-throughput transcriptomic analysis of PFAS-exposed primary human liver cell spheroids; herein, we determine the transcriptomic potency of PFAS in mixtures. Gene expression data from single PFAS and mixture exposures of liver cell spheroids were subject to benchmark concentration (BMC) analysis. We used the 25^th lowest gene BMC as the point of departure to compare the potencies of single PFAS to PFAS mixtures of varying complexity and composure. Specifically, the empirical potency of eight PFAS mixtures were compared to predicted mixture potencies calculated using the principal of concentration addition (i.e., dose addition) in which mixture component potencies are summed by proportion to predict mixture potency. In this study, for most mixtures, empirical mixture potencies were comparable to potencies calculated through concentration addition. This work supports that the effects of PFAS mixtures on gene expression largely follow the concentration addition predicted response and suggests that effects of these individual PFAS in mixtures are not strongly synergistic or antagonistic.

Data was derived as follows, 

Text copied with permission from:

Addicks GC, Rowan-Carroll A, Reardon AJF, Leingartner K, Williams A, Meier MJ, Moffat I, Carrier R, Lorusso L, Wetmore BA, et al. 2023. Per- and polyfluoroalkyl substances (PFAS) in mixtures show additive effects on transcriptomic points of departure in human liver spheroids. Toxicol Sci 194: 38–52.

 DOI: 10.1093/toxsci/kfad044

https://academic.oup.com/toxsci/article/194/1/38/7169149

Cell Culture 

3D InSightTM Human Liver Microtissues were purchased from InSphero (Brunswick, ME) in a 96 well format, with a single spheroid per well. These spheroids are a co-culture model from 10 different human liver donors including males and females and are a metabolically active system of hepatocytes and Kupffer cells. (Proctor et al. 2017; Rowan-Carroll et al. 2021). Upon arrival, culture media were replaced with InSphero Human Liver Maintenance Medium–Tox (InSphero (Brunswick, ME), and spheroids were acclimated at 37°C and 5% CO2 for 24 hours prior to PFAS exposures. PFAS were added to the media at the indicated concentrations and cells were exposed for either 24 hours or 10 days at 37°C and 5% CO2. For 10-day exposures, spent media were replaced every three days with new PFAS containing media. At the end of the exposures, media were collected and diluted 1:10 in storage buffer (200 mM Tris-HCl pH 7.3 containing 10% glycerol and 1% bovine serum albumin) then frozen at -80C to test the cytotoxicity. Spheroids were then washed once with Dulbecco’s phosphate buffered saline (DPBS) (Thermo Fisher Scientific, Franklin, MA) and lysed with 5–7 µl of TempO-Seq™ lysis buffer (BioSpyder Technologies Inc, Carlsbad, CA). Samples were triturated, incubated for 10 min at Room Temperature and then stored at -80C.

Chemical Preparation and Exposure Conditions

Generation of raw transcriptomic data for single PFAS exposures that were used for comparison of mixture exposure data were the subject of our previous investigations (Reardon et al. 2021; Rowan-Carroll et al. 2021), with the details of PFAS purchase, preparation and concentration selection for individual PFAS exposures discussed therein. Details of individual PFAS exposures are briefly described herein for clarity. For individual PFAS exposures, the highest concentration of 100 µM was based on the EPA’s ToxCast program’s highest concentration. Additional exposures were selected to capture the response over three orders of magnitude and were based on the results obtained for PFOS in our first experiment (Rowan-Carroll et al. 2021). Single PFAS exposure levels were: 0.2, 2, 10, 20, 50, 100 µM unless otherwise stated (Table 1, bottom).

Specific mixture concentrations were selected with considerations of reducing operational complexity of exposure experiments while using exposures at a comparable range to single PFAS exposures for computational modeling of concentration-response. All mixtures contained equal molarities of each PFAS in the mixture. Mixtures, exposure concentrations (the total molarity of all combined PFAS in the mixtures) were of a similar range as the single PFAS exposures in our previous studies (Reardon et al. 2021). For each mixture, with the exception of Mixture 7, there were six exposure concentrations ranging from less than 1 µM up to 100 µM. Mixture 7 had five exposure concentrations with a top exposure of 40 µM. Concentration levels were 0.02, 0.2, 1, 2, 5 µM of each PFAS for high complexity mixtures (mixtures with more than three PFAS, (Mixture 2 – 5)) and 0.2, 1, 2, 10, 20 µM for each PFAS for low complexity mixtures (mixtures with two or three PFAS, (Mixture 1, 6 and 7)). All mixtures, with the exception of Mixture 7, also had a highest concentration level with the total combined molarities of component PFAS summing to 100 µM (Table 1). To allow comparisons of PFAS mixture potencies to single PFAS potencies, all mixture concentrations are reported as total molarity of all PFAS within the mixture. Thus, total combined molarity of PFAS for each mixture exposure level was dependent on the number of PFAS in each mixture; i.e., mixture 6, with three component PFAS (PFOA, PFOS, PFNA), had exposures of 0.6, 3, 6, 30, 60 and 100 µM total PFAS (0.2, 1, 2, 10, 20 and 33 µM of each PFAS). See also Table 1 for detailed summary of mixtures and the individual and total molarities of PFAS within the mixtures. See Figure 1 for a graphical depiction of the PFAS included in each mixture.

PFOS (95% purity CAS 1763-23-1), PFOA (95% purity CAS 335-67-1) and perfluorobutanesulfonic acid (PFBS) (95% purity CAS 375-73-5) were purchased from Sigma-Aldrich (Oakville, Ontario). The remainder of PFAS were obtained through a collaboration with the U.S. Environmental Protection Agency (EPA). These PFAS were procured under EPA contract (#EP-D-12-034) by Evotec Inc (Bradford, CT, USA) with a minimum target purity concentration of 95%. These were solubilized in dimethyl sulfoxide (DMSO) and monitored to ensure no precipitation was evident. If solubility was not an issue, solutions were prepared at 30 mM; otherwise, stocks were prepared at 10 or 20 mM. EPA PFAS stocks used passed an analytical quality evaluation and were deemed to be stable and free from contaminants (Smeltz et al. 2023) The full list of PFAS used in this study and the abbreviations used throughout the manuscript are summarized in Supplemental Table 1. PFAS were dissolved in DMSO (Sigma-Aldrich, Oakville, Ontario) to prepare working stock solutions up to 30 µM. Final DMSO concentrations in cell culture media were 0.1% for PFAS concentrations below 50 µM, 0.17% for 50 µM and 0.3% for 100 µM PFAS concentrations. PFAS exposures were matched to vehicle only (DMSO) time-matched and concentration-matched controls.

The overall design included twenty 96-well plates (10 plates per time- point, (24-hour and 10-day of exposures)) of liver cell spheroids. For each time point, four separate plates contained the complete range of PFAS exposure concentrations (i.e., each exposure had one replicate on each of four separate plates in the same experiment, resulting in a total of four replicates) alongside at least two matched DMSO controls (described below). Mixture exposures occurred in the same experiment along with all other PFAS except for PFOS, PFOA, PFBS, and PFDS from Rowan-Carroll et al. (2021), which were run separately and had duplicate exposures on each of two plates (i.e., each exposure had two replicates on two separate plates in the same experiment, with a total of four replicates). Within each plate there were 10 DMSO controls: two at 0.3% for 100 µM PFAS concentrations, two at 0.17% for 50 µM PFAS concentrations, and eight at 0.1% DMSO for all other PFAS concentrations at each timepoint (totaling eight, eight and 24, respectively, across all plates). DMSO controls for plates with PFOS, PFOA, PFBS, and PFDS exposures totaled 16 at 0.1%, 8 at 0.17% and 8 at 0.3% for each timepoint. All mention of replicates within the current manuscript and supplementary information refer to replicates of PFAS-exposed microtissues (as described above) with independent downstream processing (described below).

Cytotoxicity Assessment and Cytotoxicity Exclusion Criteria

Cytotoxicity was determined using a lactate dehydrogenase (LDH) assay (LDH-Glo, Promega J2380, Madison,WI), as per the manufacturer’s instructions. Briefly, at the end of PFAS exposures cell culture media were collected and diluted 1:10 in storage buffer (200 mM Tris-HCl pH 7.3 containing 10% glycerol and 1% bovine serum albumin) and frozen at -80C. For analysis, samples were diluted 1:1 with LDH Detection Reagent and equilibrated for one hour at RT before reading relative luminescence units (RLU) (GloMax 96 Microplate Luminometer, Promega Corp, Madison,WI). Prior to analysis, RLU values from blank controls were subtracted from spheroid media RLU values. LDH ratios for each exposure were calculated by dividing experimental RLUs by the averaged RLU of respective DMSO controls. Means and standard deviations for each set of exposures (chemical, concentration, and timepoint) were then calculated from LDH ratios and expressed as fold change. Sample exclusion due to cytotoxicity was set at a 10-fold increase in LDH over controls (Rowan-Carroll et al. 2021). PFAS exposures at concentrations above the first cytotoxic concentration often resulted in decreased LDH signal, presumably due to cell death and spheroid degradation; therefore, all concentrations higher than the lowest concentrations yielding a 10-fold increase in LDH were also deemed cytotoxic. Samples that were determined to be cytotoxic based on the above metrics were excluded from all downstream analysis in order to ensure the fidelity of subsequent computational quality control of sequencing data (described below).

TempO-Seq™ Library Building and Next Generation Sequencing

Gene expression was measured using the human TempO-Seq™ S1500+ panel (House et al., 2017; Mav et al., 2018) (BioSpyder Technologies Inc, Carlsbad, CA). This panel of approximately 3000 genes was selected by the National Institute of Environmental Health Sciences to cover a diverse set of biological pathways and provide a representative list of genes that captures transcriptomic variation and toxicological response with a high level of representation when compared to RNA-seq (Bushel et al. 2018) (National Toxicology Program 2015) ( https://www.federalregister.gov/documents/2015/04/15/2015-08529/list-of-environmentally-responsive-human-genes-selected-for-use-in-screening-large-numbers-of ). For TempO-Seq™ analysis, liver microtissues were lysed as described above using a volume of 2 x TempO-Seq™ lysis buffer equal to the residual volume of DPBS. Lysates and positive controls for sequencing reactions (1 x Human Universal Reference RNA— uhrRNA (Agilent Cat # 740000) and 1 x Human Brain Total RNA brRNA (ThermoFisher AM7962)), as well as 1 x negative controls for sequencing reactions (1 x TempO-Seq™ lysis buffer alone) were hybridized with detector oligo mix according to the manufacturer’s protocol (Tempo-Seq™ Human Tox +Surrogate with a Standard Attenuation Transcriptome Kit (96 Samples)) (BioSpyder Technologies, Inc. Carlsbad, CA). Positive and negative controls were included on each plate. Hybridization was followed by nuclease digestion of excess oligos, detector oligo ligation, and amplification with tagged primers according to manufacturer’s instructions. During amplification, each sample was ligated to sample-specific barcodes to allow identification of sample sequences after pooled sequencing reactions. Labeled and pooled amplicons were column purified using NucleoSpin Gel and PCR Clean-up kits, (Takara Bio USA, Inc, Mountain View, CA). Libraries were sequenced in-house, at Health Canada, using a NextSeq 500 High-Throughput Sequencing System (Illumina, San Diego, CA) using 50 cycles from a 75-cycle high throughput flow cell.

Data Processing : Generation of Gene Expression Data

Data processing was done with R v.3.6.1 (R. Core Team 2022).

To process TempO-Seq™ data, FASTQ files were generated from the BCL files using bcl2fastq v. 2.20.0.422 (Illumina, San Diego, CA). FASTQ files were processed using the TempO-SeqR script v 3.0 provided by BioSpyder (BioSpyder Technologies, Inc. Carlsbad, CA), as implemented within our transcriptomics data processing pipeline (https://github.com/R-ODAF/R-ODAF_Health_Canada). Briefly, reads from the FASTQ files were aligned to the TempO-Seq™ Human Surrogate+Tox Panel (S1500+) v2.0  probes reference sequences (Mav et al. 2018) using STAR 2.7.8a. The qCount function from QuasR (Gaidatzis et al., 2015) was used to extract feature counts specified in a GTF file (provided by BioSpyder) from the aligned reads. The result of this workflow is a table of counts per probe per sample. The gene expression data set is available through the NCBI Gene Expression Omnibus (series numbers GSE145239 and GSE144775).

Sequencing data for individual samples underwent several rounds of quality control for inclusion or exclusion from use in downstream toxicological modeling. The study-wide quality control workflow used in this study was adapted from the recommendations made by Harrill et al. (2021) and was used to assess the quality of the alignments and exclude samples if necessary. This data analysis pipeline was an improved implementation of that used in our prior articles (Reardon et al. 2021; Rowan-Carroll et al. 2021) with higher quality control (QC) stringency. Specifically, samples with read counts below 10% of our target depth of 1 M aligned reads (i.e., 100,000 reads) were removed; this quality control step flagged both PFAS exposed and control samples. Hierarchical clustering plots were generated (hclust function: default linkage function of hclust function in R; complete-linkage) for all the samples per time point using a distance metric defined as 1-Spearman correlation in order to identify potential outliers. Samples that clustered as singletons when cutting the dendrograms at the 0.1 dissimilarity were removed from the study. We also excluded samples based on the fraction of mapped reads (samples were excluded if the alignment rate was <40%). Additionally, for each sample, we calculated the number of probes with at least 5 uniquely mapped reads; the number of probes required to capture 80% of the signal in a given sample; and the Gini coefficient. For those three metrics, any samples identified as outliers based on Tukey’s outer fence (3X interquartile range) were excluded (Harrill et al. 2021). Some samples were excluded based on multiple criteria. Further details on pipeline analysis, QA/QC, and exclusion criteria for all excluded samples are detailed in the Supplemental Rmd File 1 (provided here). Sample exclusion due to cytotoxicity and QA/QC is summarized in Figure 3. A minimum of two samples per chemical/mixture for each exposure concentration were required for benchmark concentration modeling (BMC), therefore any exposures with only one sample after QC exclusion were also excluded from BMC analysis.

For downstream analysis (e.g., BMD modeling), the count matrix (genes X samples) of all samples passing filters were log₂-transformed and normalized by their library size scaling factor derived from the median-of-ratios method in DESeq2 to account for differences in the number of reads per sample (Love et al. 2014). To account for differences in the percentage of DMSO used for the dissolution of different chemical concentrations, the log₂-transformed and DESeq2 normalized data for exposures to higher concentration PFAS that were dissolved in 0.17% and 0.30% DMSO were further normalized to their corresponding DMSO-matched controls and rescaled to the average of the DMSO 0.1% dose group. After normalizing to the controls, the DMSO 0.17% and the DMSO 0.30% controls were then removed.

Log₂-transformed DESeq2 normalized gene counts in a format appropriate for analysis with BMDExpress v3 software package (Phillips et al. 2019), that is freely available at https://github.com/auerbachs/BMDExpress-3/releases are provided here.

Buick JK, Williams A, Meier MJ, Swartz CD, Recio L, Gagné R, Ferguson SS, Engelward BP, Yauk CL. 2021. A Modern Genotoxicity Testing Paradigm: Integration of the High-Throughput CometChip® and the TGx-DDI Transcriptomic Biomarker in Human HepaRG^TM Cell Cultures. Front Public Health 9: 1144.

Bushel PR, Paules RS, Auerbach SS. 2018. A Comparison of the TempO-Seq S1500+ Platform to RNA-Seq and Microarray Using Rat Liver Mode of Action Samples. Front Genet 9: 485.

Harrill JA, Everett LJ, Haggard DE, Sheffield T, Bundy JL, Willis CM, Thomas RS, Shah I, Judson RS. 2021. High-Throughput Transcriptomics Platform for Screening Environmental Chemicals. Toxicol Sci 181: 68–89.

Jassal B, Matthews L, Viteri G, Gong C, Lorente P, Fabregat A, Sidiropoulos K, Cook J, Gillespie M, Haw R, et al. 2020. The reactome pathway knowledgebase. Nucleic Acids Res 48: D498–D503.

Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15: 550.

Mav D, Shah RR, Howard BE, Auerbach SS, Bushel PR, Collins JB, Gerhold DL, Judson RS, Karmaus AL, Maull EA, et al. 2018. A hybrid gene selection approach to create the S1500+ targeted gene sets for use in high-throughput transcriptomics. PLOS ONE 13: e0191105.

National Toxicology Program. 2015. Federal Register :: List of Environmentally Responsive Human Genes Selected for Use In Screening Large Numbers of Substances Using Toxicogenomic Approaches. https://www.federalregister.gov/documents/2015/04/15/2015-08529/list-of-environmentally-responsive-human-genes-selected-for-use-in-screening-large-numbers-of (Accessed December 12, 2022).

Phillips JR, Svoboda DL, Tandon A, Patel S, Sedykh A, Mav D, Kuo B, Yauk CL, Yang L, Thomas RS, et al. 2019. BMDExpress 2: enhanced transcriptomic dose-response analysis workflow ed. A. Valencia. Bioinformatics 35: 1780–1782.

Proctor WR, Foster AJ, Vogt J, Summers C, Middleton B, Pilling MA, Shienson D, Kijanska M, Ströbel S, Kelm JM, et al. 2017. Utility of spherical human liver microtissues for prediction of clinical drug-induced liver injury. Arch Toxicol 91: 2849–2863.

R. Core Team. 2022. R: A Language and Environment for Statistical Computing. https://www.R-project.org/.

Reardon AJF, Rowan-Carroll A, Ferguson SS, Leingartner K, Gagne R, Kuo B, Williams A, Lorusso L, Bourdon-Lacombe JA, Carrier R, et al. 2021. Potency Ranking of Per- and Polyfluoroalkyl Substances Using High-Throughput Transcriptomic Analysis of Human Liver Spheroids. Toxicol Sci. https://doi.org/10.1093/toxsci/kfab102 (Accessed October 19, 2021).

Rowan-Carroll A, Reardon A, Leingartner K, Gagné R, Williams A, Meier MJ, Kuo B, Bourdon-Lacombe J, Moffat I, Carrier R, et al. 2021. High-Throughput Transcriptomic Analysis of Human Primary Hepatocyte Spheroids Exposed to Per- and Polyfluoroalkyl Substances as a Platform for Relative Potency Characterization. Toxicol Sci 181: 199–214.

Smeltz MG, Clifton MS, Henderson WM, McMillan L, Wetmore BA. 2023. Targeted Per- and Polyfluoroalkyl substances (PFAS) assessments for high throughput screening: Analytical and testing considerations to inform a PFAS stock quality evaluation framework. Toxicol Appl Pharmacol 459: 116355.

United States Environmental Protection Agency. 2022. Benchmark Dose Software (BMDS). https://www.epa.gov/bmds/download-bmds.

Yang L, Allen BC, Thomas RS. 2007. BMDExpress: a software tool for the benchmark dose analyses of genomic data. BMC Genomics 8: 387.