R-ODAF: Omics data analysis framework for regulatory application

House dust-derived mixtures of organophosphate esters alter the phenotype, function, transcriptome, and lipidome of KGN human ovarian granulosa cells

Organophosphate esters (OPEs), used as flame retardants and plasticizers, are present ubiquitously in the environment. Previous studies suggest that exposure to OPEs is detrimental to female fertility in humans. However, no experimental information is available on the effects of OPE mixtures on ovarian granulosa cells, which play essential roles in female reproduction. We used high-content imaging to investigate the effects of environmentally relevant OPE mixtures on KGN human granulosa cell phenotypes. Perturbations to steroidogenesis were assessed using ELISA and qRT-PCR. A high-throughput transcriptomic approach, TempO-Seq, was used to identify transcriptional changes in a targeted panel of genes. Effects on lipid homeostasis were explored using a cholesterol assay and global lipidomic profiling. OPE mixtures altered multiple phenotypic features of KGN cells, with triaryl OPEs in the mixture showing higher potencies than other mixture components. The mixtures increased basal production of steroid hormones; this was mediated by significant changes in the expression of critical transcripts involved in steroidogenesis. Further, the total-OPE mixture disrupted cholesterol homeostasis and the composition of intracellular lipid droplets. Exposure to complex mixtures of OPEs, similar to those found in house dust, may adversely affect female reproductive health by altering a multitude of phenotypic and functional endpoints in granulosa cells. This study provides novel insights into the mechanisms of actions underlying the toxicity induced by OPEs and highlights the need to examine the effects of human relevant chemical mixtures.

Empirical Characterization of False Discovery Rates of Differentially Expressed Genes and Transcriptomic Benchmark Concentrations in Zebrafish Embryos

High-throughput transcriptomics (HTTr) is increasingly applied to zebrafish embryos to survey the toxicological effects of environmental chemicals. Before the adoption of this approach in regulatory testing, it is essential to characterize background noise in order to guide experimental designs. We thus empirically quantified the HTTr false discovery rate (FDR) across different embryo pool sizes, sample sizes, and concentration groups for toxicology studies. We exposed zebrafish embryos to 0.1% dimethyl sulfoxide (DMSO) for 5 days. Pools of 1, 5, 10, and 20 embryos were created (n = 24 samples for each pool size). Samples were sequenced on the TempO-Seq platform and then randomly assigned to mock treatment groups before differentially expressed gene (DEG), pathway, and benchmark concentration (BMC) analyses. Given that all samples were treated with DMSO, any significant DEGs, pathways, or BMCs are false positives. As expected, we found decreasing FDRs for DEG and pathway analyses with increasing pool and sample sizes. Similarly, FDRs for BMC analyses decreased with increasing pool size and concentration groups, with more stringent BMC premodel filtering reducing BMC FDRs. Our study provides foundational data for determining appropriate experiment designs for regulatory toxicity testing with HTTr in zebrafish embryos.

Unlocking the Power of Transcriptomic Biomarkers in Qualitative and Quantitative Genotoxicity Assessment of Chemicals

To modernize genotoxicity assessment and reduce reliance on experimental animals, new approach methodologies (NAMs) that provide human-relevant dose-response data are needed. Two transcriptomic biomarkers, GENOMARK and TGx-DDI, have shown a high classification accuracy for genotoxicity. As these biomarkers were extracted from different training sets, we investigated whether combining the two biomarkers in a human-derived metabolically competent cell line (i.e., HepaRG) provides complementary information for the classification of genotoxic hazard identification and potency ranking. First, the applicability of GENOMARK to TempO-Seq, a high-throughput transcriptomic technology, was evaluated. HepaRG cells were exposed for 72 h to increasing concentrations of 10 chemicals (i.e., eight known in vivo genotoxicants and two in vivo nongenotoxicants). Gene expression data were generated using the TempO-Seq technology. We found a prediction performance of 100%, confirming the applicability of GENOMARK to TempO-Seq. Classification using TGx-DDI was then compared to GENOMARK. For the chemicals identified as genotoxic, benchmark concentration modeling was conducted to perform potency ranking. The high concordance observed for both hazard classification and potency ranking by GENOMARK and TGx-DDI highlights the value of integrating these NAMs in a weight of evidence evaluation of genotoxicity.

Toxicological Mechanisms and Potencies of Organophosphate Esters in KGN Human Ovarian Granulosa Cells as Revealed by High-throughput Transcriptomics

Despite the growing number of studies reporting potential risks associated with exposure to organophosphate esters (OPEs), their molecular mechanisms of action remain poorly defined. We used the high-throughput TempO-Seq™ platform to investigate the effects of frequently detected OPEs on the expression of ∼3000 environmentally responsive genes in KGN human ovarian granulosa cells. Cells were exposed for 48 h to one of five OPEs (0.1 to 50 μM): tris(methylphenyl) phosphate (TMPP), isopropylated triphenyl phosphate (IPPP), tert-butylphenyl diphenyl phosphate (BPDP), triphenyl phosphate (TPHP), or tris(2-butoxyethyl) phosphate (TBOEP). The sequencing data indicate that four OPEs induced transcriptional changes, whereas TBOEP had no effect within the concentration range tested. Multiple pathway databases were used to predict alterations in biological processes based on differentially expressed genes. At lower concentrations, inhibition of the cholesterol biosynthetic pathway was the predominant effect of OPEs; this was likely a consequence of intracellular cholesterol accumulation. At higher concentrations, BPDP and TPHP had distinct effects, primarily affecting pathways involved in cell cycle progression and other stress responses. Benchmark concentration (BMC) modelling revealed that BPDP had the lowest transcriptomic point of departure. However, in vitro to in vivo extrapolation modeling indicated that TMPP was bioactive at lower concentrations than the other OPEs. We conclude that these new approach methodologies provide information on the mechanism(s) underlying the effects of data-poor compounds and assist in the derivation of protective points of departure for use in chemical read-across and decision-making.

Investigation of the effects of phthalates on in vitro thyroid models with RNA-Seq and ATAC-Seq

Introduction

Phthalates are a class of endocrine-disrupting chemicals that have been shown to negatively correlate with thyroid hormone serum levels in humans and to cause a state of hyperactivity in the thyroid. However, their mechanism of action is not well described at the molecular level.

Methods

We analyzed the response of mouse thyroid organoids to the exposure to a biologically relevant dose range of the phthalates bis(2-ethylhexyl) phthalate (DEHP), di-iso-decylphthalate (DIDP), di-iso-nonylphthalate (DINP), and di-n-octylphthalate (DnOP) for 24 h and simultaneously analyzed mRNA and miRNA expression via RNA sequencing. Using the expression data, we performed differential expression and gene set enrichment analysis. We also exposed the human thyroid follicular epithelial cell line Nthy-ori 3-1 to 1 µM of DEHP or DINP for 5 days and analyzed changes in chromatin accessibility via ATAC-Seq.

Results

Dose-series analysis showed how the expression of several genes increased or decreased at the highest dose tested. As expected with the low dosing scheme, the compounds induced a modest response on the transcriptome, as we identified changes in only mmu-miR-143-3p after DINP treatment and very few differentially expressed genes. No effect was observed on thyroid markers. Ing5, a component of histones H3 and H4 acetylation complexes, was consistently upregulated in three out of four conditions compared to control, and we observed a partial overlap among the genes differentially expressed by the treatments. Gene set enrichment analysis showed enrichment in the treatment samples of the fatty acid metabolism pathway and in the control of pathways related to the receptor signalling and extracellular matrix organization. ATAC-Seq analysis showed a general increase in accessibility compared to the control, however we did not identify significant changes in accessibility in the identified regions.

Discussion

With this work, we showed that despite having only a few differentially expressed genes, diverse analysis methods could be applied to retrieve relevant information on phthalates, showing the potential of in vitro thyroid-relevant systems for the analysis of endocrine disruptors.

Applying genomics in regulatory toxicology: a report of the ECETOC workshop on omics threshold on non-adversity

In a joint effort involving scientists from academia, industry and regulatory agencies, ECETOC's activities in Omics have led to conceptual proposals for: (1) A framework that assures data quality for reporting and inclusion of Omics data in regulatory assessments; and (2) an approach to robustly quantify these data, prior to interpretation for regulatory use. In continuation of these activities this workshop explored and identified areas of need to facilitate robust interpretation of such data in the context of deriving points of departure (POD) for risk assessment and determining an adverse change from normal variation. ECETOC was amongst the first to systematically explore the application of Omics methods, now incorporated into the group of methods known as New Approach Methodologies (NAMs), to regulatory toxicology. This support has been in the form of both projects (primarily with CEFIC/LRI) and workshops. Outputs have led to projects included in the workplan of the Extended Advisory Group on Molecular Screening and Toxicogenomics (EAGMST) group of the Organisation for Economic Co-operation and Development (OECD) and to the drafting of OECD Guidance Documents for Omics data reporting, with potentially more to follow on data transformation and interpretation. The current workshop was the last in a series of technical methods development workshops, with a sub-focus on the derivation of a POD from Omics data. Workshop presentations demonstrated that Omics data developed within robust frameworks for both scientific data generation and analysis can be used to derive a POD. The issue of noise in the data was discussed as an important consideration for identifying robust Omics changes and deriving a POD. Such variability or "noise" can comprise technical or biological variation within a dataset and should clearly be distinguished from homeostatic responses. Adverse outcome pathways (AOPs) were considered a useful framework on which to assemble Omics methods, and a number of case examples were presented in illustration of this point. What is apparent is that high dimension data will always be subject to varying processing pipelines and hence interpretation, depending on the context they are used in. Yet, they can provide valuable input for regulatory toxicology, with the pre-condition being robust methods for the collection and processing of data together with a comprehensive description how the data were interpreted, and conclusions reached.

Long-Term Use of Amoxicillin Is Associated with Changes in Gene Expression and DNA Methylation in Patients with Low Back Pain and Modic Changes

Long-term antibiotics are prescribed for a variety of medical conditions, recently including low back pain with Modic changes. The molecular impact of such treatment is unknown. We conducted longitudinal transcriptome and epigenome analyses in patients (n = 100) receiving amoxicillin treatment or placebo for 100 days in the Antibiotics in Modic Changes (AIM) study. Gene expression and DNA methylation were investigated at a genome-wide level at screening, after 100 days of treatment, and at one-year follow-up. We identified intra-individual longitudinal changes in gene expression and DNA methylation in patients receiving amoxicillin, while few changes were observed in patients receiving placebo. After 100 days of amoxicillin treatment, 28 genes were significantly differentially expressed, including the downregulation of 19 immunoglobulin genes. At one-year follow-up, the expression levels were still not completely restored. The significant changes in DNA methylation (n = 4548 CpGs) were mainly increased methylation levels between 100 days and one-year follow-up. Hence, the effects on gene expression occurred predominantly during treatment, while the effects on DNA methylation occurred after treatment. In conclusion, unrecognized side effects of long-term amoxicillin treatment were revealed, as alterations were observed in both gene expression and DNA methylation that lasted long after the end of treatment.

From vision toward best practices: Evaluating in vitro transcriptomic points of departure for application in risk assessment using a uniform workflow

The growing number of chemicals in the current consumer and industrial markets presents a major challenge for regulatory programs faced with the need to assess the potential risks they pose to human and ecological health. The increasing demand for hazard and risk assessment of chemicals currently exceeds the capacity to produce the toxicity data necessary for regulatory decision making, and the applied data is commonly generated using traditional approaches with animal models that have limited context in terms of human relevance. This scenario provides the opportunity to implement novel, more efficient strategies for risk assessment purposes. This study aims to increase confidence in the implementation of new approach methods in a risk assessment context by using a parallel analysis to identify data gaps in current experimental designs, reveal the limitations of common approaches deriving transcriptomic points of departure, and demonstrate the strengths in using high-throughput transcriptomics (HTTr) to derive practical endpoints. A uniform workflow was applied across six curated gene expression datasets from concentration-response studies containing 117 diverse chemicals, three cell types, and a range of exposure durations, to determine tPODs based on gene expression profiles. After benchmark concentration modeling, a range of approaches was used to determine consistent and reliable tPODs. High-throughput toxicokinetics were employed to translate in vitro tPODs (µM) to human-relevant administered equivalent doses (AEDs, mg/kg-bw/day). The tPODs from most chemicals had AEDs that were lower (i.e., more conservative) than apical PODs in the US EPA CompTox chemical dashboard, suggesting in vitro tPODs would be protective of potential effects on human health. An assessment of multiple data points for single chemicals revealed that longer exposure duration and varied cell culture systems (e.g., 3D vs. 2D) lead to a decreased tPOD value that indicated increased chemical potency. Seven chemicals were flagged as outliers when comparing the ratio of tPOD to traditional POD, thus indicating they require further assessment to better understand their hazard potential. Our findings build confidence in the use of tPODs but also reveal data gaps that must be addressed prior to their adoption to support risk assessment applications.

Dataset on DNA methylation and gene expression changes induced by 5-aza-2′-deoxycytidine in Syrian golden hamster fetal cell cultures

The Syrian hamster (SH) is an animal model used in virology, toxicology, and carcinogenesis, where a better understanding of epigenetic mechanisms is required. Finding genetic loci regulated by DNA methylation may assist in the development of DNA methylation-based in vitro assays for the identification of carcinogens. This dataset informs on the regulation of gene expression by DNA methylation. Primary cultures of SH male fetal cells (sex determined by differences in kdm5 loci on the X and Y chromosome) were exposed for 7 days to the carcinogen benzo[a]pyrene (20 µM) from which a morphologically transformed colony was collected and reseeded. The colony bypassed senescence and sustained growth. After 210 days of culture, the cells were collected and divided in 16 aliquots to create 4 experimental groups to test the effects of the DNA methylation inhibitor 5-aza-2'-deoxycytidine (5adC). The experiment was initiated 24 h after cell seeding in 10 cm plates. The groups are naïve cells (N), cells exposed for 48 h to either 0.05% DMSO as vehicle (V), or to 5adC at 1 µM and 5 µM. DNA and RNA libraries were sequenced on an Illumina NextSeq 500. Gene expression was analysed by RNAseq and differentially methylated DNA regions (DMRs: clusters of 200 base pairs (bp), read depth >20, q< 0.05, methylation difference >|25%|) were identified by reduce representation bisulfite sequencing (RRBS). Global genome DNA methylation was similar between the N (mean±SD, 47.3%±0.02) and V groups (47.3%±0.01). Although 5adC reduced methylation, the reduction was larger in the 1 µM (39.2%±0.002) than in the 5 µM group (44.3%±0.01). 5adC induced a total of 612 and 190 DMRs by 1 µM and 5 µM, among which 79 and 23 were in the promoter regions (±3,000 bp from the transcription start site), respectively. 5adC induced a total of 1,170 and 1,797 differentially expressed genes (DEGs) by 1 µM and 5 µM, respectively. The 5 µM treatment induced statistically significant toxicity (% cell viability: group N 97%±8, V 98.8%±1.3, 1 µM 97.3%±0.5, 5 µM 93.8%±1.5), which perhaps reduced cell division and daughter cell numbers with inherited changes in methylation, but increased number of DEGs due to both toxicity and methylation changes. As usually observed in the literature, a small portion of DEGs (4% and 4% at 1 µM and 5 µM, respectively) are associated with DMRs in their promoters. These promoter DMRs by themselves are sufficient among other epigenetic marks to induce DEGs. The dataset provides the genomic coordinates of the DMRs and an opportunity to further examine their roles in distal putative promoters or enhancers (yet to be described in the SH) in contributing to gene expression changes, senescence bypass and sustained proliferation as essential carcinogenic events (see companion paper [1]). Finally, this experiment confirms the possibility in future experiments to use 5adC as a positive control for effects on DNA methylation in cells derived from SH.

In vitro transcriptomic analyses reveal pathway perturbations, estrogenic activities, and potencies of data-poor BPA alternative chemicals

Since initial regulatory action in 2010 in Canada, bisphenol A (BPA) has been progressively replaced by structurally related alternative chemicals. Unfortunately, many of these chemicals are data-poor, limiting toxicological risk assessment. We used high-throughput transcriptomics to evaluate potential hazards and compare potencies of BPA and 15 BPA alternative chemicals in cultured breast cancer cells. MCF-7 cells were exposed to BPA and 15 alternative chemicals (0.0005-100 µM) for 48 h. TempO-Seq (BioSpyder Inc) was used to examine global transcriptomic changes and estrogen receptor alpha (ERα)-associated transcriptional changes. Benchmark concentration (BMC) analysis was conducted to identify 2 global transcriptomic points of departure: (1) the lowest pathway median gene BMC and (2) the 25th lowest rank-ordered gene BMC. ERα activation was evaluated using a published transcriptomic biomarker and an ERα-specific transcriptomic point of departure was derived. Genes fitting BMC models were subjected to upstream regulator and canonical pathway analysis in Ingenuity Pathway Analysis. Biomarker analysis identified BPA and 8 alternative chemicals as ERα active. Global and ERα transcriptomic points of departure produced highly similar potency rankings with bisphenol AF as the most potent chemical tested, followed by BPA and bisphenol C. Further, BPA and transcriptionally active alternative chemicals enriched similar gene sets associated with increased cell division and cancer-related processes. These data provide support for future read-across applications of transcriptomic profiling for risk assessment of data-poor chemicals and suggest that several BPA alternative chemicals may cause hazards at similar concentrations to BPA.

CODA: a combo-Seq data analysis workflow

The analysis of the combined mRNA and miRNA content of a biological sample can be of interest for answering several research questions, like biomarkers discovery, or mRNA-miRNA interactions. However, the process is costly and time-consuming, separate libraries need to be prepared and sequenced on different flowcells. Combo-Seq is a library prep kit that allows us to prepare combined mRNA-miRNA libraries starting from very low total RNA. To date, no dedicated bioinformatics method exists for the processing of Combo-Seq data. In this paper, we describe CODA (Combo-seq Data Analysis), a workflow specifically developed for the processing of Combo-Seq data that employs existing free-to-use tools. We compare CODA with exceRpt, the pipeline suggested by the kit manufacturer for this purpose. We also evaluate how Combo-Seq libraries analysed with CODA perform compared with conventional poly(A) and small RNA libraries prepared from the same samples. We show that using CODA more successfully trimmed reads are recovered compared with exceRpt, and the difference is more dramatic with short sequencing reads. We demonstrate how Combo-Seq identifies as many genes and fewer miRNAs compared to the standard libraries, and how miRNA validation favours conventional small RNA libraries over Combo-Seq. The CODA code is available at https://github.com/marta-nazzari/CODA.

Multi-omics HeCaToS dataset of repeated dose toxicity for cardiotoxic & hepatotoxic compounds

The data currently described was generated within the EU/FP7 HeCaToS project (Hepatic and Cardiac Toxicity Systems modeling). The project aimed to develop an in silico prediction system to contribute to drug safety assessment for humans. For this purpose, multi-omics data of repeated dose toxicity were obtained for 10 hepatotoxic and 10 cardiotoxic compounds. Most data were gained from in vitro experiments in which 3D microtissues (either hepatic or cardiac) were exposed to a therapeutic (physiologically relevant concentrations calculated through PBPK-modeling) or a toxic dosing profile (IC20 after 7 days). Exposures lasted for 14 days and samples were obtained at 7 time points (therapeutic doses: 2-8-24-72-168-240-336 h; toxic doses 0-2-8-24-72-168-240 h). Transcriptomics (RNA sequencing & microRNA sequencing), proteomics (LC-MS), epigenomics (MeDIP sequencing) and metabolomics (LC-MS & NMR) data were obtained from these samples. Furthermore, functional endpoints (ATP content, Caspase3/7 and O2 consumption) were measured in exposed microtissues. Additionally, multi-omics data from human biopsies from patients are available. This data is now being released to the scientific community through the BioStudies data repository (https://www.ebi.ac.uk/biostudies/).

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.