Exposure to PFAS chemicals induces sex-dependent alterations in key rate-limiting steps of lipid metabolism in liver steatosis.

Toxicants with the potential to bioaccumulate in humans and animals have long been a cause for concern, particularly due to their association with multiple diseases and organ injuries. Per- and polyfluoro alkyl substances (PFAS) and polycyclic aromatic hydrocarbons (PAH) are two such classes of chemicals that bioaccumulate and have been associated with steatosis in the liver. Although PFAS and PAH are classified as chemicals of concern, their molecular mechanisms of toxicity remain to be explored in detail. In this study, we aimed to identify potential mechanisms by which an acute exposure to PFAS and PAH chemicals can induce lipid accumulation and whether the responses depend on chemical class, dose, and sex. To this end, we analyzed mechanisms beginning with the binding of the chemical to a molecular initiating event (MIE) and the consequent transcriptomic alterations. We collated potential MIEs using predictions from our previously developed ToxProfiler tool and from published steatosis adverse outcome pathways. Most of the MIEs are transcription factors, and we collected their target genes by mining the TRRUST database. To analyze the effects of PFAS and PAH on the steatosis mechanisms, we performed a computational MIE-target gene analysis on high-throughput transcriptomic measurements of liver tissue from male and female rats exposed to either a PFAS or PAH. The results showed peroxisome proliferator-activated receptor (PPAR)-α targets to be the most dysregulated, with most of the genes being upregulated. Furthermore, PFAS exposure disrupted several lipid metabolism genes, including upregulation of fatty acid oxidation genes (Acadm, Acox1, Cpt2, Cyp4a1-3) and downregulation of lipid transport genes (Apoa1, Apoa5, Pltp). We also identified multiple genes with sex-specific behavior. Notably, the rate-limiting genes of gluconeogenesis (Pck1) and bile acid synthesis (Cyp7a1) were specifically downregulated in male rats compared to female rats, while the rate-limiting gene of lipid synthesis (Scd) showed a PFAS-specific upregulation. The results suggest that the PPAR signaling pathway plays a major role in PFAS-induced lipid accumulation in rats. Together, these results show that PFAS exposure induces a sex-specific multi-factorial mechanism involving rate-limiting genes of gluconeogenesis and bile acid synthesis that could lead to activation of an adverse outcome pathway for steatosis.

Transcriptional pathways linked to fetal and maternal hepatic dysfunction caused by gestational exposure to perfluorooctanoic acid (PFOA) or hexafluoropropylene oxide-dimer acid (HFPO-DA or GenX) in CD-1 mice.

Per- and polyfluoroalkyl substances (PFAS) comprise a diverse class of chemicals used in industrial processes, consumer products, and fire-fighting foams which have become environmental pollutants of concern due to their persistence, ubiquity, and associations with adverse human health outcomes, including in pregnant persons and their offspring. Multiple PFAS are associated with adverse liver outcomes in adult humans and toxicological models, but effects on the developing liver are not fully described. Here we performed transcriptomic analyses in the mouse to investigate the molecular mechanisms of hepatic toxicity in the dam and its fetus after exposure to two different PFAS, perfluorooctanoic acid (PFOA) and its replacement, hexafluoropropylene oxide-dimer acid (HFPO-DA, known as GenX). Pregnant CD-1 mice were exposed via oral gavage from embryonic day (E) 1.5-17.5 to PFOA (0, 1, or 5 mg/kg-d) or GenX (0, 2, or 10 mg/kg-d). Maternal and fetal liver RNA was isolated (N = 5 per dose/group) and the transcriptome analyzed by Affymetrix Array. Differentially expressed genes (DEG) and differentially enriched pathways (DEP) were obtained. DEG patterns were similar in maternal liver for 5 mg/kg PFOA, 2 mg/kg GenX, and 10 mg/kg GenX (R2: 0.46-0.66). DEG patterns were similar across all 4 dose groups in fetal liver (R2: 0.59-0.81). There were more DEGs in fetal liver compared to maternal liver at the low doses for both PFOA (fetal = 69, maternal = 8) and GenX (fetal = 154, maternal = 93). Upregulated DEPs identified across all groups included Fatty Acid Metabolism, Peroxisome, Oxidative Phosphorylation, Adipogenesis, and Bile Acid Metabolism. Transcriptome-phenotype correlation analyses demonstrated > 1000 maternal liver DEGs were significantly correlated with maternal relative liver weight (R2 >0.92). These findings show shared biological pathways of liver toxicity for PFOA and GenX in maternal and fetal livers in CD-1 mice. The limited overlap in specific DEGs between the dam and fetus suggests the developing liver responds differently than the adult liver to these chemical stressors. This work helps define mechanisms of hepatic toxicity of two structurally unique PFAS and may help predict latent consequences of developmental exposure.

Impact of Aligner, Normalization Method, and Sequencing Depth on TempO-seq Accuracy.

High-throughput transcriptomics has advanced through the introduction of TempO-seq, a targeted alternative to traditional RNA-seq. TempO-seq platforms use 50 nucleotide probes, each specifically designed to target a known transcript, thus allowing for reduced sequencing depth per sample compared with RNA-seq without compromising the accuracy of results. Thus far, studies using the TempO-seq method have relied on existing tools for processing the resulting short read data. However, these tools were originally designed for other data types. While they have been used for processing of early TempO-seq data, they have not been systematically assessed for accuracy or compared to determine an optimal framework for processing and analyzing TempO-seq data. In this work, we re-analyze several publicly available TempO-seq data sets covering a range of experimental designs and use corresponding RNA-seq data sets as a gold standard to rigorously assess accuracy at multiple levels. We compare 6 aligners and 5 normalization methods across various accuracy and performance metrics. Our results demonstrate the overall robust accuracy of the TempO-seq platform, independent of data processing methods. Complex aligners and advanced normalization methods do not appear to have any general advantage over simpler methods when it comes to analyzing TempO-seq data. The reduced complexity of the sequencing space, and the fact that TempO-seq probes are all equal length, appears to reduce the need for elaborate bioinformatic or statistical methods used to address these factors in RNA-seq data.

Whole genome sequencing of low input circulating cell-free DNA obtained from normal human subjects.

Cell-free DNA circulates in plasma at low levels as a normal by-product of cellular apoptosis. Multiple clinical pathologies, as well as environmental stressors can lead to increased circulating cell-free DNA (ccfDNA) levels. Plasma DNA studies frequently employ targeted amplicon deep sequencing platforms due to limited concentrations (ng/ml) of ccfDNA in the blood. Here, we report whole genome sequencing (WGS) and read distribution across chromosomes of ccfDNA extracted from two human plasma samples from normal, healthy subjects, representative of limited clinical samples at <1 ml. Amplification was sufficiently robust with ~90% of the reference genome (GRCh38.p2) exhibiting 10X coverage. Chromosome read coverage was uniform and directly proportional to the number of reads for each chromosome across both samples. Almost 99% of the identified genomic sequence variants were known annotated dbSNP variants in the hg38 reference genome. A high prevalence of C>T and T>C mutations was present along with a strong concordance of variants shared between the germline genome databases; gnomAD (81.1%) and the 1000 Genome Project (93.6%). This study demonstrates isolation and amplification procedures from low input ccfDNA samples that can detect sequence variants across the whole genome from amplified human plasma ccfDNA that can translate to multiple clinical research disciplines.

Utility of Extrapolating Human S1500+ Genes to the Whole Transcriptome: Tunicamycin Case Study.

The TempO-Seq S1500+ platform(s), now available for human, mouse, rat, and zebrafish, measures a discrete number of genes that are representative of biological and pathway co-regulation across the entire genome in a given species. While measurement of these genes alone provides a direct assessment of gene expression activity, extrapolating expression values to the whole transcriptome (~26 000 genes in humans) can estimate measurements of non-measured genes of interest and increases the power of pathway analysis algorithms by using a larger background gene expression space. Here, we use data from primary hepatocytes of 54 donors that were treated with the endoplasmic reticulum (ER) stress inducer tunicamycin and then measured on the human S1500+ platform containing ~3000 representative genes. Measurements for the S1500+ genes were then used to extrapolate expression values for the remaining human transcriptome. As a case study of the improved downstream analysis achieved by extrapolation, the "measured only" and "whole transcriptome" (measured + extrapolated) gene sets were compared. Extrapolation increased the number of significant genes by 49%, bringing to the forefront many that are known to be associated with tunicamycin exposure. The extrapolation procedure also correctly identified established tunicamycin-related functional pathways reflected by coordinated changes in interrelated genes while maintaining the sample variability observed from the "measured only" genes. Extrapolation improved the gene- and pathway-level biological interpretations for a variety of downstream applications, including differential expression analysis, gene set enrichment pathway analysis, DAVID keyword analysis, Ingenuity Pathway Analysis, and NextBio correlated compound analysis. The extrapolated data highlight the role of metabolism/metabolic pathways, the ER, immune response, and the unfolded protein response, each of which are key activities associated with tunicamycin exposure that were unrepresented or underrepresented in one or more of the analyses of the original "measured only" dataset. Furthermore, the inclusion of the extrapolated genes raised "tunicamycin" from third to first upstream regulator in Ingenuity Pathway Analysis and from sixth to second most correlated compound in NextBio analysis. Therefore, our case study suggests an approach to extend and enhance data from the S1500+ platform for improved insight into biological mechanisms and functional outcomes of diseases, drugs, and other perturbations.

Evaluation of 5-day In Vivo Rat Liver and Kidney With High-throughput Transcriptomics for Estimating Benchmark Doses of Apical Outcomes.

A 5-day in vivo rat model was evaluated as an approach to estimate chemical exposures that may pose minimal risk by comparing benchmark dose (BMD) values for transcriptional changes in the liver and kidney to BMD values for toxicological endpoints from traditional toxicity studies. Eighteen chemicals, most having been tested by the National Toxicology Program in 2-year bioassays, were evaluated. Some of these chemicals are potent hepatotoxicants (eg, DE71, PFOA, and furan) in rodents, some exhibit toxicity but have minimal hepatic effects (eg, acrylamide and α,ß-thujone), and some exhibit little overt toxicity (eg, ginseng and milk thistle extract) based on traditional toxicological evaluations. Male Sprague Dawley rats were exposed once daily for 5 consecutive days by oral gavage to 8-10 dose levels for each chemical. Liver and kidney were collected 24 h after the final exposure and total RNA was assayed using high-throughput transcriptomics (HTT) with the rat S1500+ platform. HTT data were analyzed using BMD Express 2 to determine transcriptional gene set BMD values. BMDS was used to determine BMD values for histopathological effects from chronic or subchronic toxicity studies. For many of the chemicals, the lowest transcriptional BMDs from the 5-day assays were within a factor of 5 of the lowest histopathological BMDs from the toxicity studies. These data suggest that using HTT in a 5-day in vivo rat model provides reasonable estimates of BMD values for traditional apical endpoints. This approach may be useful to prioritize chemicals for further testing while providing actionable data in a timely and cost-effective manner.

KRAS-retroviral fusion transcripts and gene amplification in arsenic-transformed, human prostate CAsE-PE cancer cells.

CAsE-PE cells are an arsenic-transformed, human prostate epithelial line containing oncogenic mutations in KRAS compared to immortalized, normal KRAS parent cells, RWPE-1. We previously reported increased copy number of mutated KRAS in CAsE-PE cells, suggesting gene amplification. Here, KRAS flanking genomic and transcriptomic regions were sequenced in CAsE-PE cells for insight into KRAS amplification. Comparison of DNA-Seq and RNA-Seq showed increased reads from background aligning to all KRAS exons in CAsE-PE cells, while a uniform DNA-Seq read distribution occurred in RWPE-1 cells with normal transcript expression. We searched for KRAS fusions in DNA and RNA sequencing data finding a portion of reads aligning to KRAS and viral sequence. After generation of cDNA from total RNA, short and long KRAS probes were generated to hybridize cDNA and KRAS enriched fragments were PacBio sequenced. More KRAS reads were captured from CAsE-PE cDNA versus RWPE-1 by each probe set. Only CAsE-PE cDNA showed KRAS viral fusion transcripts, primarily mapping to LTR and endogenous retrovirus sequences on either 5'- or 3'-ends of KRAS. Most KRAS viral fusion transcripts contained 4 to 6 exons but some PacBio sequences were in unusual orientations, suggesting viral insertions within the gene body. Additionally, conditioned media was extracted for potential retroviral particles. RNA-Seq of culture media isolates identified KRAS retroviral fusion transcripts in CAsE-PE media only. Truncated KRAS transcripts suggested multiple retroviral integration sites occurred within the KRAS gene producing KRAS retroviral fusions of various lengths. Findings suggest activation of endogenous retroviruses in arsenic carcinogenesis should be explored.

Development of a Zebrafish S1500+ Sentinel Gene Set for High-Throughput Transcriptomics.

Sentinel gene sets have been developed with the purpose of maximizing the information from targeted transcriptomic platforms. We recently described the development of an S1500+ sentinel gene set, which was built for the human transcriptome, utilizing a data- and knowledge-driven hybrid approach to select a small subset of genes that optimally capture transcriptional diversity, correlation with other genes based on large-scale expression profiling, and known pathway annotation within the human genome. While this detailed bioinformatics approach for gene selection can in principle be applied to other species, the reliability of the resulting gene set depends on availability of a large body of transcriptomics data. For the model organism zebrafish, we aimed to create a similar sentinel gene set (Zf S1500+ gene set); however, there is insufficient standardized expression data in the public domain to train the gene correlation model. Therefore, our strategy was to use human-zebrafish ortholog mapping of the human S1500+ genes and nominations from experts in the zebrafish scientific community. In this study, we present the bioinformatics curation and refinement process to produce the final Zf S1500+ gene set, explore whole transcriptome extrapolation using this gene set, and assess pathway-level inference. This gene set will add value to targeted high-throughput transcriptomics in zebrafish for toxicogenomic screening and other research domains.

Arsenite malignantly transforms human prostate epithelial cells in vitro by gene amplification of mutated KRAS.

Inorganic arsenic is an environmental human carcinogen of several organs including the urinary tract. RWPE-1 cells are immortalized, non-tumorigenic, human prostate epithelia that become malignantly transformed into the CAsE-PE line after continuous in vitro exposure to 5µM arsenite over a period of months. For insight into in vitro arsenite transformation, we performed RNA-seq for differential gene expression and targeted sequencing of KRAS. We report >7,000 differentially expressed transcripts in CAsE-PE cells compared to RWPE-1 cells at >2-fold change, q<0.05 by RNA-seq. Notably, KRAS expression was highly elevated in CAsE-PE cells, with pathway analysis supporting increased cell proliferation, cell motility, survival and cancer pathways. Targeted DNA sequencing of KRAS revealed a mutant specific allelic imbalance, 'MASI', frequently found in primary clinical tumors. We found high expression of a mutated KRAS transcript carrying oncogenic mutations at codons 12 and 59 and many silent mutations, accompanied by lower expression of a wild-type allele. Parallel cultures of RWPE-1 cells retained a wild-type KRAS genotype. Copy number analysis and sequencing showed amplification of the mutant KRAS allele. KRAS is expressed as two splice variants, KRAS4a and KRAS4b, where variant 4b is more prevalent in normal cells compared to greater levels of variant 4a seen in tumor cells. 454 Roche sequencing measured KRAS variants in each cell type. We found KRAS4a as the predominant transcript variant in CAsE-PE cells compared to KRAS4b, the variant expressed primarily in RWPE-1 cells and in normal prostate, early passage, primary epithelial cells. Overall, gene expression data were consistent with KRAS-driven proliferation pathways found in spontaneous tumors and malignantly transformed cell lines. Arsenite is recognized as an important environmental carcinogen, but it is not a direct mutagen. Further investigations into this in vitro transformation model will focus on genomic events that cause arsenite-mediated mutation and overexpression of KRAS in CAsE-PE cells.

Whole exome sequencing in the rat.

BACKGROUND: The rat genome was sequenced in 2004 with the aim to improve human health altered by disease and environmental influences through gene discovery and animal model validation. Here, we report development and testing of a probe set for whole exome sequencing (WES) to detect sequence variants in exons and UTRs of the rat genome. Using an in-silico approach, we designed probes targeting the rat exome and compared captured mutations in cancer-related genes from four chemically induced rat tumor cell lines (C6, FAT7, DSL-6A/C1, NBTII) to validated cancer genes in the human database, Catalogue of Somatic Mutations in Cancer (COSMIC) as well as normal rat DNA. Paired, fresh frozen (FF) and formalin-fixed, paraffin-embedded (FFPE) liver tissue from naive rats were sequenced to confirm known dbSNP variants and identify any additional variants. RESULTS: Informatics analysis of available gene annotation from rat RGSC6.0/rn6 RefSeq and Ensembl transcripts provided 223,636 unique exons representing a total of 26,365 unique genes and untranslated regions. Using this annotation and the Rn6 reference genome, an in-silico probe design generated 826,878 probe sequences of which 94.2% were uniquely aligned to the rat genome without mismatches. Further informatics analysis revealed 25,249 genes (95.8%) covered by at least one probe and 23,603 genes (93.5%) had every exon covered by one or more probes. We report high performance metrics from exome sequencing of our probe set and Sanger validation of annotated, highly relevant, cancer gene mutations as cataloged in the human COSMIC database, in addition to several exonic variants in cancer-related genes. CONCLUSIONS: An in-silico probe set was designed to enrich the rat exome from isolated DNA. The platform was tested on rat tumor cell lines and normal FF and FFPE liver tissue. The method effectively captured target exome regions in the test DNA samples with exceptional sensitivity and specificity to obtain reliable sequencing data representing variants that are likely chemically induced somatic mutations. Genomic discovery conducted by means of high throughput WES queries should benefit investigators in discovering rat genomic variants in disease etiology and in furthering human translational research.

HAfTs are novel lncRNA transcripts from aflatoxin exposure.

The transcriptome can reveal insights into precancer biology. We recently conducted RNA-Seq analysis on liver RNA from male rats exposed to the carcinogen, aflatoxin B1 (AFB1), for 90 days prior to liver tumor onset. Among >1,000 differentially expressed transcripts, several novel, unannotated Cufflinks-assembled transcripts, or HAfTs (Hepatic Aflatoxin Transcripts) were found. We hypothesized PCR-cloning and RACE (rapid amplification of cDNA ends) could further HAfT identification. Sanger data was obtained for 6 transcripts by PCR and 16 transcripts by 5'- and 3'-RACE. BLAST alignments showed, with two exceptions, HAfT transcripts were lncRNAs, >200nt without apparent long open reading frames. Six rat HAfT transcripts were classified as 'novel' without RefSeq annotation. Sequence alignment and genomic synteny showed each rat lncRNA had a homologous locus in the mouse genome and over half had homologous loci in the human genome, including at least two loci (and possibly three others) that were previously unannotated. While HAfT functions are not yet clear, coregulatory roles may be possible from their adjacent orientation to known coding genes with altered expression that include 8 HAfT-gene pairs. For example, a unique rat HAfT, homologous to Pvt1, was adjacent to known genes controlling cell proliferation. Additionally, PCR and RACE Sanger sequencing showed many alternative splice variants and refinements of exon sequences compared to Cufflinks assembled transcripts and gene prediction algorithms. Presence of multiple splice variants and short tandem repeats found in some HAfTs may be consequential for secondary structure, transcriptional regulation, and function. In summary, we report novel, differentially expressed lncRNAs after exposure to the genotoxicant, AFB1, prior to neoplastic lesions. Complete cloning and sequencing of such transcripts could pave the way for a new set of sensitive and early prediction markers for chemical hepatocarcinogens.

RNA-Seq-based toxicogenomic assessment of fresh frozen and formalin-fixed tissues yields similar mechanistic insights.

Formalin-fixed, paraffin-embedded (FFPE) pathology specimens represent a potentially vast resource for transcriptomic-based biomarker discovery. We present here a comparison of results from a whole transcriptome RNA-Seq analysis of RNA extracted from fresh frozen and FFPE livers. The samples were derived from rats exposed to aflatoxin B1 (AFB1 ) and a corresponding set of control animals. Principal components analysis indicated that samples were separated in the two groups representing presence or absence of chemical exposure, both in fresh frozen and FFPE sample types. Sixty-five percent of the differentially expressed transcripts (AFB1 vs. controls) in fresh frozen samples were also differentially expressed in FFPE samples (overlap significance: P < 0.0001). Genomic signature and gene set analysis of AFB1 differentially expressed transcript lists indicated highly similar results between fresh frozen and FFPE at the level of chemogenomic signatures (i.e., single chemical/dose/duration elicited transcriptomic signatures), mechanistic and pathology signatures, biological processes, canonical pathways and transcription factor networks. Overall, our results suggest that similar hypotheses about the biological mechanism of toxicity would be formulated from fresh frozen and FFPE samples. These results indicate that phenotypically anchored archival specimens represent a potentially informative resource for signature-based biomarker discovery and mechanistic characterization of toxicity.

RNA-Seq profiling reveals novel hepatic gene expression pattern in aflatoxin B1 treated rats.

Deep sequencing was used to investigate the subchronic effects of 1 ppm aflatoxin B1 (AFB1), a potent hepatocarcinogen, on the male rat liver transcriptome prior to onset of histopathological lesions or tumors. We hypothesized RNA-Seq would reveal more differentially expressed genes (DEG) than microarray analysis, including low copy and novel transcripts related to AFB1's carcinogenic activity compared to feed controls (CTRL). Paired-end reads were mapped to the rat genome (Rn4) with TopHat and further analyzed by DESeq and Cufflinks-Cuffdiff pipelines to identify differentially expressed transcripts, new exons and unannotated transcripts. PCA and cluster analysis of DEGs showed clear separation between AFB1 and CTRL treatments and concordance among group replicates. qPCR of eight high and medium DEGs and three low DEGs showed good comparability among RNA-Seq and microarray transcripts. DESeq analysis identified 1,026 differentially expressed transcripts at greater than two-fold change (p<0.005) compared to 626 transcripts by microarray due to base pair resolution of transcripts by RNA-Seq, probe placement within transcripts or an absence of probes to detect novel transcripts, splice variants and exons. Pathway analysis among DEGs revealed signaling of Ahr, Nrf2, GSH, xenobiotic, cell cycle, extracellular matrix, and cell differentiation networks consistent with pathways leading to AFB1 carcinogenesis, including almost 200 upregulated transcripts controlled by E2f1-related pathways related to kinetochore structure, mitotic spindle assembly and tissue remodeling. We report 49 novel, differentially-expressed transcripts including confirmation by PCR-cloning of two unique, unannotated, hepatic AFB1-responsive transcripts (HAfT's) on chromosomes 1.q55 and 15.q11, overexpressed by 10 to 25-fold. Several potentially novel exons were found and exon refinements were made including AFB1 exon-specific induction of homologous family members, Ugt1a6 and Ugt1a7c. We find the rat transcriptome contains many previously unidentified, AFB1-responsive exons and transcripts supporting RNA-Seq's capabilities to provide new insights into AFB1-mediated gene expression leading to hepatocellular carcinoma.