Integrating gene expression and splicing dynamics across dose-response oxidative modulators

Toxicological risk assessment increasingly utilizes transcriptomics to derive point of departure (POD) and modes of action (MOA) for chemicals. One essential biological process that allows a single gene to generate several different RNA isoforms is called alternative splicing. To comprehensively assess the role of splicing dysregulation in toxicological evaluation and elucidate its potential as a complementary endpoint, we performed RNA-seq on A549 cells treated with five oxidative stress modulators across a wide dose range. Differential gene expression (DGE) showed limited pathway enrichment except at high concentrations. However, alternative splicing analysis revealed variable intron retention events affecting diverse pathways for all chemicals in the absence of significant expression changes. For instance, diazinon elicited negligible gene expression changes but progressive increase in the number of intron retention events, suggesting splicing alterations precede expression responses. Benchmark dose modeling of intron retention data highlighted relevant pathways overlooked by expression analysis. Systematic integration of splicing datasets should be a useful addition to the toxicogenomic toolkit. Combining both modalities paint a more complete picture of transcriptomic dose-responses. Overall, evaluating intron retention dynamics afforded by toxicogenomics may provide biomarkers that can enhance chemical risk assessment and regulatory decision making. This work highlights splicing-aware toxicogenomics as a possible additional tool for examining cellular responses.

Introduction

Genomic and transcriptomic approaches in toxicology provide unprecedented molecular insight into the mode-of-action (MOA) and derivation of points-of-departure (PODs) while reducing time, increasing throughput, and lowering the cost (Chepelev et al., 2015; Harrill et al., 2019; 2021; Baltazar et al., 2020; Franzosa et al., 2021). Thus, these approaches are becoming increasingly well-accepted as alternatives to conventional in vivo studies to derive PODs for assessing chemical safety in screening studies (Pagé-Larivière et al., 2019; Baltazar et al., 2020; Johnson et al., 2020; Pouzou et al., 2020). Recently, a focus of transcriptomics in toxicology has been the dose-response modeling of gene expression. This application couples benchmark dose (BMD) modeling to derive transcriptomic POD values determined by the gene expression changes induced by the chemical and gene ontology pathway analysis (EFSA Scientific Committee et al., 2017; Farmahin et al., 2017; Jensen et al., 2019). This approach identifies specific POD levels that are altered due to chemical exposure. The National Toxicology Program (NTP) report outlines an approach for genomic dose-response that has been effective for in vivo studies, where the biological effect POD closely approximates apical toxicology endpoints. For in vitro studies, additional steps like in vitro to in vivo extrapolation (IVIVE) modeling may be needed before the transcriptomic results can be effectively used for quantitative risk assessment (National Toxicology Program, 2018). Despite the advent of bioinformatic approaches and pipelines to determine the MOA and POD values of a chemical by transcriptomics, there is not yet any set consensus on data exploration. Therefore, it is crucial to conduct a systematic analysis of the bioinformatic approaches employed for deriving BMD and POD values.

The commonalities and differences observed in gene ontology pathways between chemicals can be attributed to their shared ability to induce oxidative stress, a well-documented mechanism in toxicology. For instance, pathways related to mitochondrial function, cell respiration, and DNA damage repair are often affected by oxidative stress as the cell attempts to counteract the damage caused by reactive oxygen species (ROS). Additionally, chemicals may target specific pathways or processes depending on their unique chemical structures and MOA. To further elucidate these relationships, additional studies and experiments are needed to provide a more detailed understanding of how oxidative stress mediates the observed gene ontology pathway changes in response to specific chemicals.

Alternative splicing is a fundamental cellular mechanism that enables a single gene to produce multiple distinct mRNA isoforms, thereby expanding proteomic diversity and regulating gene expression (Baralle and Giudice, 2017). Among the various types of alternative splicing events, intron retention, where introns are retained within the mature RNA, has garnered increasing attention due to its potential role in cellular processes and disease. While the significance of alternative splicing is well-established, the functional implications of intron retention remain relatively unexplored. Recent studies have begun shedding light on the regulatory roles of intron retention in diverse biological contexts, including neurodevelopment, immune responses, and cancer progression (Braunschweig et al., 2014; Wong et al., 2016; Jacob and Smith, 2017). The involvement of intron retention in toxicology and its potential as a molecular mechanism remain largely untapped.

In this study, we examine transcriptomics for short-term in vitro exposures with chemicals modulating cellular oxidative stress (Kappus, 1987). Oxidative stress is caused by an imbalance between the production and accumulation of ROS in cells and tissues and the ability of a biological system to detoxify these reactive products. ROS can be generated as by-products from the metabolism of environmental pollutants or xenobiotics that lead to imbalance and eventually tissue damage. Physiological and biological stress responses such as immune cell activation, inflammation, infection, cancer, or aging are all also responsible for the endogenous production of oxidative species. Excess ROS can adversely affect several cellular structures, such as membranes, lipids, proteins, lipoproteins, and DNA (Mittal et al., 2014). Chemical-induced oxidative stress represents an external influence on normal homeostatic processes and may exert effects at the transcriptome level differently for different chemicals that impact oxidative stress response.

Transcriptomic events occurring at doses prior to cell death can provide a clear insight into whether toxicity is an accumulation of cell damage or an acute response to high-dose exposure to a chemical. In an effort to investigate these possibilities, we examined five effectors of oxidative stress - diazinon, paraquat, prochloraz, sulforaphane, and 2,4-dinitrochlorobenzene (DNCB)—in a dose-response exposure using the A549 lung adenocarcinoma cancer cell line. The A549 lung adenocarcinoma cell line was selected because it is a well-characterized model that has been extensively used to study oxidative stress and chemical toxicity (Lukaszewicz et al., 2019; Upadhyay et al., 2019; Niechoda et al., 2023), and harbor functional cytochrome P450 enzymes and antioxidant response pathways that allow them to respond to and metabolize xenobiotics that modulate oxidative stress (Hukkanen et al., 2000). The chemicals were chosen for their varied mechanisms of impacting oxidative stress (Supplementaryu Table S1). Cell viability and RNA-seq experiments were performed using a nine-concentration dose-response experimental setup. The RNA-seq data were analyzed for differential gene expression, alternative splicing events, BMD modeling, and subsequent POD derivation following ontology pathway enrichment analyses. Taken together, these results presented in this study underscore the underutilized potential of intron retention in the context of toxicogenomics and its possible in illuminating MOAs for various chemicals. By comparing gene expression and alternative splicing events, this study examines the value of combining these tools by evaluating compounds with different characteristics that affect oxidative stress.

Materials and methods

Cell culture

A549 cells were obtained from the American Type Culture Collection (ATCC). Proliferation media consisted of Dulbecco’s Modified Essential Medium (Gibco) with 10% fetal bovine serum (FBS) (Atlanta Biologicals) and 1% penicillin-streptomycin (Gibco). The plating medium consisted of Dulbecco’s Modified Essential Medium (Gibco) with 10% Heat-Inactivated FBS (HI FBS; Atlanta Biologicals) and 1% penicillin-streptomycin (Gibco). For passaging, cells were washed with sterile Dulbecco’s Phosphate Buffered Saline (DPBS) and enzymatically removed using 0.25% trypsin-EDTA (Gibco). In preparation for the proliferation assay, A549 cells were thawed (2.0 x 106, passage 3) and placed into a T25 flask (United States Scientific). After 5 days in culture, cells were passaged into a T75 flask (2.0 × 106, passage 4) using a plating medium (10% HI FBS). After 5 days, cells were plated at a density of 35,000 cells/well in black-walled, clear-bottom, tissue-culture-treated 96-well plates (Greiner Bio-One) using a plating medium (10% HI FBS). Cells were cultured for 48 h in 96-well plates prior to treatment. Cell viability was measured using a bioluminescent assay for ATP in lysed cells (Promega G7571, CellTiter-Glo assay) with 3 biological and 3 technical replicates for each condition.

Chemicals and reagents

Test chemicals, sulforaphane (Chemical Abstracts Service Registry number (CASRN) 4478-93-7), prochloraz (CASRN 67747-09-5), paraquat (CASRN 1910-42-5), diazinon (CASRN 333-41-5), and 2,4-dinitrochlorobenzene (DNCB - CASRN 97-00-7) were procured from Sigma. All chemicals had >98% purity, except sulforaphane, which had >90% purity. Stock solutions were made up of molecular biology grade (>99.9% purity) dimethyl sulfoxide (DMSO) (Sigma) or deionized water (paraquat and sulforaphane). The stock solutions were diluted into media prior to adding to cells, keeping the solvent concentration the same in all wells. Each compound was tested in eight concentrations in half-log increments from the highest dose tested. Diazinon, prochloraz, and DNCB were dissolved in DMSO, while paraquat and sulforaphane were dissolved in sterile deionized water. DMSO vehicle control wells contained 0.1% DMSO in culture medium, matching the final DMSO concentration in chemical treatment wells. Deionized water vehicle controls were used for paraquat and sulforaphane treatments. For serial dilutions, the chemicals were diluted using the same respective solvents to maintain consistent vehicle concentrations across all doses. All treatment wells contained 500 μL of cell culture medium. Chemicals or vehicle compounds were added in a 50 μL volume to achieve the desired final concentrations. Treatment plates were prepared fresh before each experiment to avoid evaporation issues.

RNA-Seq data generation and analysis

Using 4 biological replicates for each condition, whole cell RNA-seq was performed by Novogene Corp. using polyA enrichment. The RNA isolation was performed by utilizing the RNeasy Micro Kit (Qiagen Cat. No. 74004) with polyA enrichment. The retained RNA underwent quality assessment with a RIN >9 threshold, and purity was verified based on absorbance ratio measurements (260:280 nm ratio of −2). Illumina 150bp paired-end read kits (TruSeq Stranded mRNA Library Prep Kit (Cat #20020594) as per the instructions supplied with the kits. Reads were mapped to the Ensembl human reference genome (hg38) using HISAT2 (version 2.05) (Zhang et al., 2021). Each biological replicate of each condition was sequenced to an average depth of −40 million reads (standard deviation −5.2 million reads). Preceding the mapping step, reads were subjected to trimming and quality filtering processes to enhance data accuracy and reliability. RNA-seq results were uniformly of high quality (>90% mapping rate). The replicates displayed high QC quality with regards to the mapped count data and correlation, and thus each concentration for every chemical was analyzed using all four replicated samples. The 4 biological replicates displayed high QC quality with regards to the mapped count data and correlation, and thus all 4 replicates for every condition were analyzed using DESeq2 for differential gene expression and VAST-tools for alternative splicing analysis. FeatureCounts (Liao et al., 2014), a part of the Subread package, was employed to quantify the read counts for each gene from the aligned RNA-seq data. The analysis involved specifying the aligned BAM files as input, and the resulting count matrix was utilized for downstream expression analysis.

Differential gene expression analyses

The open-source statistical program R was used to run the DESeq2 library (Love et al., 2014), and normalized log2 (gene count+1) transformed data for BMD analyses were generated with DESeq2 for all analyses. Genes that harbored less than 10 gene counts were removed from the dataset prior to DeSeq2 analysis. To capture genes altered by exposure to compound with high confidence, we applied corrected p-value (false discovery rate) of less than 0.05 (FDR<0.05) and a fold change greater than 1.5-fold (|FC|>1.5). The Venn diagrams were generated by using the DeepVenn package (Hulsen, 2022). The gene ontology analysis was performed with the gProfiler software (Kolberg et al., 2023).

Analysis of Splicing events

Alternative splicing analysis of RNA-Seq data was performed with VAST-Tools (version 2) (Tapial et al., 2017), a comprehensive software suite designed for the identification and quantification of various alternative splicing events from RNA-Seq data. VAST-Tools employs a rigorous methodology that aligns the reads to the reference genome and transcriptome, accurately identifying splice junctions and quantifying alternative splicing events, including intron retention, exon skipping, alternative 5′and 3′splice site usage, and mutually exclusive exons.

To identify differentially regulated alternative splicing events across the different chemical treatments and doses, we utilized the “diff” module of VAST-Tools. This module compares the Percentage Spliced In (PSI) values, which represent the inclusion ratio of a particular alternative splicing event, between the treatment and control conditions. We then applied a threshold of |dPSI| > 5%, which means that only alternative splicing events with an absolute difference in PSI values greater than 5% between the treatment and control conditions were considered significant.

Additionally, we employed the MV_dPSI_at0.95CI > 0 criterion, which ensures that the 95% confidence interval for the mean difference in PSI values removes non-confident events. This step further enhances the statistical robustness of the analysis by accounting for biological variability and ensuring that the observed differences in alternative splicing are statistically significant.

By combining these two criteria, we identified alternative splicing events that were consistently and significantly altered across the different chemical treatments and doses. The raw codes and usage of VAST-tools can be obtained from https://github.com/vastgroup/vast-tools.

Benchmark dose modeling

The BMDExpress2 package was used to identify different BMD models for the ontology enrichment genes and pathways (Phillips et al., 2019). A complete description of the BMD modeling approach was also previously described (Black et al., 2022). All models were fitted assuming constant variance as the data were log2 transformed after the normalization in DESeq2. A final best fitting model was determined by first determining the best fitting polynomial model by nested Chi square test. The best fit polynomial model was then compared to the remaining models to select a final overall best fitting model by Akaike information criteria (AIC). A BMD factor of 1 standard deviation (SD) was used. The option for flagging the Hill model with “k” parameter less than ⅓ of lowest positive dose was utilized, and best model selection with flagged Hill model was done by selecting the next best model with p-value >0.05. First, any model with extrapolated BMD value beyond the highest or the lowest concentrations used in the dose-response experiment was rejected. Similarly, any best fitting model with a goodness of fit p-value less than 0.1 was rejected. Finally, any model which had a BMDU/BMDL ratio greater than 40 was rejected to avoid BMD estimates with excessively large 95% confidence intervals.

Ontology enrichment was performed using the publicly available GO ontology database. BMDExpress2 performs a conventional over-representation test of the query gene lists relative to the pathway annotated gene lists in Gene Ontology. For BMD analyses, the data were prefiltered using the ANOVA test (p ≤ 0.05) and a fold-change of ≥1.5 or ≤1.5. Data were then filtered with the settings in BMDExpress v2.3 of best BMDU/BMDL <40, and best fitPvalue >0.1. We set a significance threshold for the enrichment of a pathway having at least 5 of our query genes found amongst the pathway elements, and with a Fisher’s Exact test (two-tailed) p-value from the over-representation test that was less than 0.05. A pathway BMD value is computed as the median gene based BMDU/BMD/BMDL values for the genes found amongst the pathway category elements.

For intron retention BMD modeling, we used the number of transcript-mapped reads to calculate the expression levels of individual introns, and they were modeled at the “intron” level in BMDExpress2, similar to the approach for gene expression counts. Specifically, the intron retention mapped read counts were log2 transformed and median-normalized per sample prior to BMD modeling. This matches the preprocessing for gene counts. After the BMD calculation step, the corresponding genes were used to perform the BMD pathway analyses.

Significantly enriched pathways and computation of pathway-based POD values were obtained by selecting genes from the ontology enrichment analyses. For ontology enrichment results, we summarized the 20 most sensitive enriched Gene Ontology (GO) pathways (p < 0.05, genes that passed all filters >5). Comparisons were restricted to median pathway values or averages of the median BMD values of genes with a pathway as our experience with other data (not shown) has shown these are consistently more conservative than pathway mean values.

The GO analysis presented in Figure 2 was conducted based on data from single dose concentrations, involving the categorization of genes as either up- or downregulated. This initial analysis aimed to capture early molecular responses to the various chemicals studied. In contrast, the subsequent BMD pathway analysis, conducted using BMDExpress2, took a different approach. It centered on calculating the increase or decrease in gene expression that aligns with one of the mathematical models tested, thereby providing a more comprehensive representation of dose-response relationships across the entire concentration range. By utilizing BMDExpress2, this analysis aimed to identify critical benchmark concentration values (BMDL, BMD, and BMDU) that serve as thresholds indicative of the chemicals’ impact on molecular pathways. This transition from single dose GO (molecular, biological, and functional domains) analysis to BMD pathway analysis generates a deeper and more nuanced understanding of the molecular changes induced by these oxidative stressors. The network analysis was performed by using the CytoScape software (v3.9.1).

Results

Cell viability of A549 cells to different oxidative stress inducers

In this study, we assessed 5 oxidative stress effectors using toxicogenomics analyses (Figure 1A, Supplementary Table S1). To assess the sensitivity of the A549 cell line to each chemical, and to determine the concentration range for transcriptomic studies, we performed cell viability assays (Figures 1B–F). This experiment allowed us to identify the upper concentration for each of the five chemicals, where the concentrations for the transcriptomic experiments consisted of half-log dilution series from that maximum concentration. These values were 10-4 M for Paraquat and Prochloraz, 10-5 M for Sulforaphane and DNCB, and 10-3 M for Diazinon (Figures 1B–F). For RNA-seq experiments, the maximum concentration, where the cell viability decreased by 15%–25%, was selected to prevent confounding cell death and cytotoxicity and the loss of viable RNA. The 24-h exposure was chosen as the time point since this time point displayed intermediate effects when compared to shorter and longer time points. Altogether, these experiments determined the viability rates of A549 cells upon exposure to different oxidative stress inducers and allowed us to choose the optimum dose for dose-curve transcriptomics experiments.

Figure 1

Figure 1. Cell viability upon treatment of oxidative stress modulators. (A) Schematic overview of the experimental procedure in this study. (B–F) A549 cell viability assay for (A) DNCB, (B) Sulforaphane, (C) Prochloraz, (D) Paraquat and (E) Diazinon. Blue = 1 h, red = 6 h, green = 24 h and purple = 48 h.

RNA-seq and differential gene expression

RNA-Seq results were uniformly of high quality with an average of −40 million mapped reads per replicate. The replicates displayed high QC quality with regards to the mapped count data and correlation, and thus each concentration for every chemical was analyzed using all 4 replicated samples. Differential gene expression results (see Methods for details) for A549 cells show that at highest dose concentration, there were hundreds to thousands of differentially expressed genes (DEGs) for all five chemicals. We observed a clear dose-response of DEGs with diazinon and prochloraz with the A549 cells (Figure 2A; Table 1).

Figure 2

Figure 2. Transcriptomic analysis of oxidative stress modulators. (A) Dot plot showing the number of upregulated (red) and downregulated (blue) genes for the 5 chemicals across the doses tested. (B) Heatmap showing the -logFDR value of the gene ontology pathways for the lowest chemical concentrations in which there was observed differential gene expression. (C) Upset plot showing the number of upregulated genes at the highest dose concentration for each oxidative stress modulator. The bars on the left represent the total number of upregulated genes for each compound, while the bars on top indicate the number of genes shared among different combinations of compounds. The inset shows the enriched Gene Ontology (GO) terms for the overlapping upregulated genes across most compounds, highlighting their association with oxidative stress response, and signaling pathways. (D) Upset plot showing the number of downregulated genes at the highest dose concentration for each oxidative stress modulator. The inset shows the enriched Gene Ontology (GO) terms for the overlapping downregulated genes across most compounds, highlighting their association with DNA metabolism, cell cycle processes, and chromosome organization.

Table 1

Table 1. The number of genes dysregulated for each condition.

Next, in order to assess the pathways that were altered via the treatment of each chemical, we performed GO analysis of DEGs identified with the lowest dose concentration in which significantly differentially expressed genes and pathways were identified for each compound and cell line (Figure 2B).

The 31.6 µM diazinon treatment in A549 cells revealed alteration of genes involved in GO terms such as the “apolipoprotein binding” or “cholesterol homeostasis”. 10 µM DNCB treatment revealed pathways such as “respiratory pathway IV” and “mitochondrial envelope” in A549 cells, indicating the perturbation of respiratory pathways. Next, GO analysis of 10 µM sulforaphane treatment revealed the perturbation of “protein binding” or “cell surface” pathways. The lowest doses of which prochloraz treatment showed a DEG and perturbed pathways were 10 µM for A549 cells. GO terms altered in A549 cells in this condition included “circulatory system process” “transcriptional activity” and “cell-matrix adhesion”. Finally, 31.6 µM paraquat treatment in A549 cells revealed pathways such as “long chain fatty acid metabolism” or “signaling receptor binding” in A549 cells (Figure 2B).

In summary, these findings indicate that examining lower doses where DEGs are observed unveils diverse pathways affected. An exception to this would be DNCB, which led to aberrant expression of genes within the oxidative respiratory pathway. We identified upregulation of expected oxidative stress marker genes, such as SOD1, HMOX1 and NQO1, and relevant pathways related to oxidative stress responses. Therefore, the differential gene expression and GO pathway analysis using transcriptomics reveal overlapping as well as differing number of DEGs and enriched molecular, cellular and functional GO pathway categories across the 5 oxidative stress modulators, providing a highly detailed information about the molecular changes associated with each compound (Figure 2B).

To analyze the common cellular responses between the chemicals, we assessed the extent of commonly regulated up- and downregulated genes in the highest dose condition of each chemical that are unique or shared among the different compounds (Figures 2C,D). A substantial number of genes were commonly upregulated across multiple chemicals, with a core set of 26 genes shared by all compounds except DNCB (Figure 2C). GO analysis of the overlapping upregulated genes revealed a significant enrichment of terms associated with oxidative stress response, such as “cellular response to oxygen levels,” “cellular response to decreased oxygen levels,” and “cellular response to hypoxia” (Figure 2C, inset). In contrast, the analysis of overlapping downregulated genes (Figure 2D) revealed 36 genes that commonly overlap across all compounds except DNCB. GO analysis of the overlapping downregulated genes highlighted terms related to DNA metabolism, cell cycle processes, and chromosome organization (Figure 2D, inset), indicating potential disruptions in these cellular processes.

Benchmark concentration modeling and POD derivation

To gain insight into the pathways that are affected as a function of increasing dose concentrations, we used BMD analyses to derive POD values from the transcriptomic data of each chemical. Given the varying approaches found in the literature (Harrill et al., 2021; Matteo et al., 2023) and the evolving nature of transcriptomics data analysis, we chose to conduct a comprehensive and impartial analysis by employing two distinct methods. First, we applied a significance criterion for BMD pathways, where each pathway contained more than 5 genes and displayed an FDR value of less than 0.05 (see Methods) and selected the 20 lowest BMD pathways. For the second method, we analyzed the BMD values for individual intron retention events.

First, we applied a significance criterion for BMD pathways, where each pathway contained more than 5 genes and displayed an FDR value of less than 0.05 (see Methods) and selected the 20 lowest BMD pathways. This approach led to hundreds of pathways being identified (Figure 3B). Analysis of the average lower benchmark dose (BMDL), BMD and upper benchmark dose (BMDU) values of the lowest 20 pathways revealed values that were consistent with the cell viability assays (Figure 1). Since only sulforaphane and DNCB had two enriched pathways with A549 cells, a single summary POD was derived based on the median value for all pathways for these results. The average BMDL values for the 20 most significant pathways were compared to the average BMDL values for all significant pathways. In all cases, the average BMDL for the 20 most significant pathways was lower than the average BMDL for the full set of significant pathways. We observed BMDL values of 19.9, 3.3, 2.43, 0.5, and 0.29 μM for diazinon, paraquat, prochloraz, DNCB, and sulforaphane, respectively. As a result, this analysis highlighted several known features of these compounds and unearthed unexpected MOA’s (see Discussion).

Figure 3

Figure 3. BMD analysis using gene expression. (A) BMD Median Accumulation plot for the genes in each tested condition. The ANOVA test (p < 0.05) and a fold-change of >1.5 or < −1.5 were used to pre filter the data. In BMDExpress v2.3, data were also post-filtered for best BMDU/BMD 40 and best fitPvalue >.1. (B) Upset plot showing the number of BMD pathways that are either unique to each condition or overlapping with multiple conditions. The lines that connect the dots below indicate an overlap between the conditions, and the bar graph on top indicates the number of overlapping pathways. The bar graph on the left indicates the number of observed BMD pathways for each condition. (C) Network graph showing the connectivity of the BMD pathways that display the highest rate of overlap among the tested conditions.

Next, we visualized the BMD accumulation plots for genes that exhibited responses below the maximum tested dose for each condition (Figure 3A). A549 cells treated with DNCB displayed the gene distribution with the lowest BMD accumulation, whereas cells treated with diazinon had the highest (Figure 3A). We next analyzed the BMD pathways that were commonly shared among all the chemical conditions and identified that pathways such as “regulation of serine/threonine kinase activity” or “metabolic” and “apoptotic process” that were overlapping (Figures 3B,C, Supplementary Figure S2). A wide range of cell cycle and DNA damage endpoints were most predominant among the pathways observed at the low concentration range. Notable cellular responses to oxidative or chemical stress pathways were apparent in cells exposed to diazinon, paraquat, prochloraz, and DNCB at the BMD and BMDL levels observed. Only paraquat demonstrated major chemical oxidative stress response pathways observed.

Our analysis highlights shared toxicity pathways while also revealing some substantial differences between the chemicals (Supplementary Figure S2). Core cell cycle control and DNA damage pathways showed sensitivity to all compounds, implying general reactive toxicity. However, paraquat uniquely disrupted MAPK signaling (Liu et al., 2022), potentially reflecting oxidative inhibition of phosphatases. Diazinon’s impacts on immune pathways align with its acetylcholinesterase inhibition (Colović et al., 2013). Sulforaphane and DNCB showed earlier perturbation of translation and RNA processing. Prochloraz’s low-dose effects on reproductive pathways match its endocrine disruption (Vandenberg et al., 2012). Thus, despite common downstream outcomes like cell cycle dysfunction, signatures of distinct molecular initiating events emerged. Overall, combined assessment of shared and unique BMD-perturbed pathways facilitates MOA elucidation, aiding chemical characterization and risk analysis.

These analyses altogether highlight the specific and shared pathways that are perturbed with each chemical and underscore the importance of careful BMD analysis pipelines for the identification of MOAs.

Alternative splicing analysis

Alternative splicing generates multiple mRNA isoforms from a single gene through mechanisms like exon skipping, intron retention, and alternate splice site usage (Figure 4A), allowing individual genes to produce multiple transcript isoforms, expanding proteomic diversity (Blencowe, 2006; Lee and Rio, 2015; Tapial et al., 2017; Park et al., 2018; Ule and Blencowe, 2019). This process is regulated by splice sites, splicing factors, and RNA-binding proteins that control exon inclusion/exclusion and intron excision. Major splicing patterns include cassette exon skipping, intron retention, and use of alternate 5′/3′ splice sites (Figure 4A). Alternative splicing can influence the relative abundance of different transcript isoforms originating from the same gene, thereby impacting the composition of the transcript pool without necessarily altering the overall gene expression levels. Despite its prevalence in eukaryotes, and its frequent perturbation in several developmental diseases and cancer (Scotti and Swanson, 2015), the role of alternative splicing in toxicology remains under-explored (Zaharieva et al., 2012; Villaseñor-Altamirano et al., 2019; Banerjee et al., 2020).

Figure 4

Figure 4. Splicing changes occur prior to changes in gene expression. (A) Schematic representation of the major alternative splicing events observed. (B) Bar plot showing the number of exon skipping (EX), intron retention (IR), alternative 3′ splice site (ALT3), and alternative 5′ splice site (ALT5) events along with the number of differentially expressed genes (DEG) for the five chemicals across the doses. The splicing events are observed at concentration where there aren’t any DEGs.

Here, we utilized the RNA-seq data to systematically investigate splice variants associated with each chemical (see Methods) using the widely utilized VAST-Tools software (Tapial et al., 2017). VAST-Tools quantifies alternative splicing events by counting RNA-seq reads mapping uniquely to splice junctions and involved exons, then calculates an inclusion ratio for each event representing the relative abundance of isoforms with or without the event. Alternative splicing analysis revealed extensive changes in events such as exon skipping (EX) and intron retention (IR) across the tested chemicals. All chemicals elicited robust changes in the four splicing event types often at low doses (Figure 4B, Supplementary Figure S3A). For instance, diazinon triggered the inclusion of a high number of IR and EX events by 0.316 μM, preceding toxicity. DNCB similarly resulted in marked IR and EX changes at 0.001 μM, below concentrations impairing viability or gene expression. Most compounds displayed preferential modulation of certain event types, illustrating the complexity of splicing regulation. However, IR and EX were consistently the predominant response, emerging as generalizable markers of chemical perturbation. In contrast, ALT3 and ALT5 events were detected in lower numbers. Thus, it may be possible that all types of splicing events could be playing a combinatorial role in oxidative stress. While many chemicals elicited splicing perturbations without significant differential expression, sulforaphane uniquely induced gene expression alterations alongside splicing effects. Overall, alternative splicing analysis provides another window to query chemical effects overlooked by standard expression profiling.

Benchmark concentration modeling using intron retention events

Of the splicing alterations elicited by chemical toxicants, intron retention (IR) and exon skipping (EX) events were the most common splicing alterations across all compounds tested, highlighting it as a possible marker of cellular stress (Vu et al., 2022). Though under-appreciated relative to exon skipping, IR is now recognized to play key regulatory roles, influencing mRNA localization, stability, translation, and protein production (Wong et al., 2016; Jacob and Smith, 2017; Monteuuis et al., 2019). Retained introns introduce premature termination codons, enabling rapid message turnover. Moreover, the unique sequences and motifs within retained introns can affect splicing, RNA localization, transcription, and translation via RNA-binding proteins and chromatin modifiers.

We therefore leveraged IR profiles to model dose-response relationships and benchmark concentrations. The BMD pathway accumulation plot revealed numerous intron retention events surpassing benchmark response cutoffs, even at low concentrations (Figures 5A,B). For instance, critical mRNA processing pathways emerged at a BMD of −4 μM for several of our test compounds. Upset plot analysis highlighted extensive overlap in BMD pathways across chemicals, centered on key processes like splicing and translation (Figure 5C). Interestingly, except for diazinon, the mean BMDs derived from intron retention events were lower than those from differential expression for the same compounds (Figures 5D–H). Accordingly, the number of BMD pathways detected by gene expression profiling partially overlapped with BMD pathways detected by IR profiling (Supplementary Figures S3B–F). This supports intron retention’s sensitivity, capturing pathway perturbations undetected by expression profiling. Overall, applying benchmark concentration modeling to global intron retention patterns revealed unexpected pathways and POD. Splicing-derived BMD modeling expands the utility of toxicogenomics.

Figure 5

Figure 5. BMD Analysis using intron retention events. (A) BMD Median Accumulation plot, generated by using the intron retention events, for each tested condition. (B) Bubble plot showing the minus log10 p-value of lowest 50 BMD pathways as a function of Median BMD levels. (C) Upset plot showing the overlap of BMD pathways determined by intron retention events. (D–H) Scatter plot showing the Mean BMD values that have been generated either by gene expression profiling (y-axis) or splicing analyses (x-axis) for (D) Diazinon, (E) Prochloraz, (F) Sulforaphane, (G) DNCB, and (H) Paraquat.

The BMD analysis uncovered biological pathways and points of departure from IR that differed from those associated with gene expression. Our results reveal retained introns as promising targets meriting extensive analysis to elucidate chemical MOAs and enhance risk assessment.

To elucidate shared and unique biological changes across chemicals, we further constructed a pathway network to visualize the overlapping BMD pathways (Figure 6). Network visualization of pathways modulated at benchmark concentration levels revealed extensive connectivity and shared biology across the mechanistically distinct chemicals. Numerous pathways were jointly perturbed by multiple compounds, centered on key processes like mRNA processing, vesicle trafficking, and post-translational modifications. For instance, critical splicing and polyadenylation pathways were common to prochloraz, DNCB, and sulforaphane, three otherwise divergent oxidative toxicants. Similarly, regulation of ubiquitination and transport vesicles emerged as conserved hubs, suggesting generalized cellular stress effects. Furthermore, DNA repair pathways unique to sulforaphane emerged at higher BMD levels, providing chemical-specific insights. However, closer examination also revealed chemical-specific interactions reflecting unique alterations, such as DNA repair for sulforaphane or neuron development for DNCB. In addition, calcium signaling pathways showed lower median BMDs for DNCB versus prochloraz, implying greater potency despite common modulation. This network approach enhances compound-specific potency assessments, elucidating unique molecular initiation events and their dose-dependencies.

Figure 6

Figure 6. Network analysis of BMD pathways determined by intron retention events. Network graph showing the unique, as well as overlapping BMD pathways that have been determined by intron retention event modeling. The color of the edges indicates the Median BMD values.

Collectively, these results demonstrate that combining BMD modeling with network biology provides both a broad overview of common toxicological mechanisms, as well as nuanced insights into the precise pathway perturbations underlying each compound’s toxicity. This systems-level perspective enhances molecule-specific risk assessment.

Discussion

In this study, we sought to identify the pathways affected in a dose response manner by various oxidative stress modulators. The range of BMD values obtained for the chemicals showed a correlative trend with the doses where ∼15–25% cell death occurred in viability experiments (Figure 1). This trend is consistent with the fact that the DEGs also followed the same trend as the viability assays (Figure 2A). GO analysis using fixed, early dose concentrations revealed several mitochondrial and cellular respiration-associated pathways consistent with the known effects of these chemicals (Figure 2B).

This study demonstrates the potential of combining gene expression and splicing analysis to elucidate chemical mechanisms of action (MOAs) and PODs. While expression profiling revealed overt toxicity thresholds, systematic splicing characterization uncovered lower-dose RNA processing disruption, affirming RNA homeostasis as a toxicological MOA. Remarkably, intron retention emerged as the primary splicing alteration, providing a generalized indicator of cellular stress across structurally diverse toxicants, as also suggested by earlier reports (Hadar et al., 2022; Vu et al., 2022). By modeling dose-dependent intron retention profiles, we extracted benchmark concentration pathways and PODs at lower concentrations than those affecting gene expression. The minimal, but significant overlap between expression- and splicing-derived pathways affirms retained introns as informative, complementary biomarkers. Rather than superseding expression results, splicing analysis may add useful information.

The modulation of splicing events, encompassing intron retention, exon skipping, and 3′or 5′splice site selection, plays important roles in gene expression regulation (Jacob and Smith, 2017). Perturbations induced by environmental factors or toxicants can lead to aberrant splicing patterns, generating diverse mRNA isoforms, some of which may encode dysfunctional proteins that disrupt cellular processes. Altered splicing can also result in the accumulation of unstable RNA species, triggering stress responses. Moreover, changes in splicing may yield protein isoforms with varying activities and interactions, contributing to cellular dysfunction and toxicity. These alterations serve as valuable biomarkers, offering insights into the molecular mechanisms underlying toxicity and providing early indicators for comprehensive risk assessment in toxicology.

Our integrated transcriptomics approach provides broad, pathway-level insights through expression analysis, supplemented by event-level perturbations from splicing. By fusing both techniques, we obtained a more complete perspective into chemical MOAs, capturing low-abundance anomalies antecedent to downstream endpoints. Both techniques provide complementary insights. Gene expression defines chemical impacts on functional protein pathways that influence cellular outcomes. Meanwhile, splicing analyses reveals the fine-tuned regulatory disruptions underlying these transcriptomic shifts. Integrating both is crucial for a mechanistic understanding of how chemicals initiate molecular perturbations that propagate into overt toxicity. Future work should balance sensitive pathway examination with systems-level perspectives and investigate whether changes in alternative splicing is a hallmark within different toxicological contexts.

We identified intricate dose-dependent modulation of alternative splicing events (Supplementary Figure S3). With Diazinon and Sulforaphane, there is a balanced up- and downregulation of IR and EX events across all the doses. DNCB induces a robust downregulation of IR and EX events across all doses. Paraquat treatment leads to IR events that are consistently upregulated especially at lower doses, with balanced distribution at higher doses. Moreover, in the context of Prochloraz treatment, there is a dose-dependent decrease in the number of upregulated and an increase in the number of downregulated IR and EX events. These observations highlight the complex interplay between different splicing events and their dose-dependent modulation, which may be influenced by the specific mechanisms of action and cellular responses elicited by each oxidative stress modulator. The clear differences in the patterns of up- and downregulation across compounds underscore the need for a comprehensive analysis of alternative splicing events to fully understand the splicing dysregulation associated with oxidative stress and its implications for cellular function, toxicity and risk assessment. Comparison of altered pathways detected by differential gene expression and IR analysis reveals several common pathways that are co-regulated by these processes (Supplementaryu Table S4).

Differential basal expression of genes and pathways involved in mediating or counteracting chemical MOAs likely contributes to cell-type specific susceptibility (Black et al., 2023). Furthermore, compounds may preferentially disrupt pathways that are highly expressed in certain cell contexts due to tissue-specific roles, while not affecting cell types lacking associated gene programs. Elucidating how chemical perturbations interface with cell-intrinsic expression landscapes will shed light on selective sensitivities and improve cross-cell-type extrapolation.

Assessing significant BMD pathways which were commonly identified to be overlapping among the cell line-chemical combinations revealed that these pathways do not necessarily harbor the lowest BMD or BMDL (i.e., BMDL = POD) levels but are shared among the oxidative stressors and the 2 cell lines used (Figure 3). These analyses revealed several converging pathways that include metabolic processes, phosphorylation of proteins (i.e., serine/threonine kinases in cell cycle control) or regulation of gene expression by gene silencing. A detailed, gene-level BMD analysis of the top three overlapping pathways (Figure 3) indicates that there is a range of genes with differing BMD values that drive the alterations of these pathways.

Intron retention-based BMCs tended to be lower than viability or gene expression-based BMCs across most chemicals. However, there is a high correlation between the median BMC values of lowest BMC pathways and those that have been derived from cell viability assays (Krebs et al., 2020) (Supplementary Figure S4). This highlights the sensitivity of splicing alterations like intron retention in capturing early molecular changes at doses below those impairing viability or gene expression. However, further evaluation using in vitro to in vivo extrapolation (IVIVE) is needed to determine if intron retention provides PODs that better reflect in vivo dose-response compared to gene expression alone. Nonetheless, our results demonstrate the value of intron retention profiling for interpreting chemical MOA through the unique pathways and dose-dependencies revealed. The non-overlapping sets of BMC pathways from gene expression vs. splicing analyses substantiate intron retention as a complementary indicator of biological impacts at lower concentrations.

Future work focusing on the validation and functional characterization of individual alternative splicing events is required to identify individual markers for oxidative stress or other toxicology applications. Follow-up work involving ROS assays and DNA damage assays like COMET to directly quantify oxidative stress upon chemical treatment will be paramount. It will be important to correlate altered splice isoform levels with these functional readouts to determine which splicing events are most associated with oxidative stress and toxicity. Knockdown or overexpression studies focused on shifting isoform ratios will provide further evidence linking alternative splicing to oxidative stress. In addition, studies characterizing the involvement of other splicing events, such as exon skipping, will be important to identify novel biomarkers.

Further visualization of the pathways indicated several processes that were specific to individual chemical-cell line combinations or shared by many chemicals (Supplementary Figure S1). For instance, A549 cells treated with DNCB showed several pathways related to oxygen transport. On the other hand, “Transcriptional Regulation by TP53” and “D-loop structure” pathways were commonly shared with many compounds. Previous research has highlighted the p53 pathways as a marker for oxidative stress (Liu et al., 2020). Similarly, D-loops structures, due to their relaxed structures, have been shown to be more prone to oxidative damage than other mitochondrial DNA regions (Rothfuss et al., 2010).

Taken together, these results demonstrate that splicing analysis may provide another complementary biomarker to standard expression profiling. The non-overlapping pathway sets highlight retained introns as a complementary indicator of biological impacts induced by chemical exposure. Overall, coupling intron retention modeling with gene expression analysis provides a systems perspective.

Data availability statement

The raw sequencing and processed data in this study have been submitted to Gene Expression Omnibus (GEO) with the accession number GSE241672 and can be accessed via this link: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE241672.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

AB: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–original draft, Writing–review and editing. MB: Data curation, Formal Analysis, Investigation, Methodology, Software, Writing–original draft. RS: Formal Analysis, Investigation, Writing–original draft. SS: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Writing–original draft. PM: Conceptualization, Data curation, Formal Analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing–original draft. AN: Conceptualization, Formal Analysis, Funding acquisition, Investigation, Writing–original draft.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the American Chemistry Council’s Long-range Research Initiative.

Acknowledgments

We would like to thank Melvin Andersen for critical feedback.

Conflict of interest

Authors ARB, MB, RS, SS, PM, and AN were employed by the ScitoVation LLC.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2024.1389095/full#supplementary-material

References

Baltazar, M. T., Cable, S., Carmichael, P. L., Cubberley, R., Cull, T., Delagrange, M., et al. (2020). A next-generation risk assessment case study for coumarin in cosmetic products. Toxicol. Sci. 176, 236–252. doi:10.1093/toxsci/kfaa048

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, M., Ferragut Cardoso, A. P., Lykoudi, A., Wilkey, D. W., Pan, J., Watson, W. H., et al. (2020). Arsenite exposure displaces zinc from ZRANB2 leading to altered splicing. Chem. Res. Toxicol. 33, 1403–1417. doi:10.1021/acs.chemrestox.9b00515

PubMed Abstract | CrossRef Full Text | Google Scholar

Baralle, F. E., and Giudice, J. (2017). Alternative splicing as a regulator of development and tissue identity. Nat. Rev. Mol. Cell Biol. 18, 437–451. doi:10.1038/nrm.2017.27

PubMed Abstract | CrossRef Full Text | Google Scholar

Black, M. B., Efremenko, A. Y., McMullen, P. D., Rasim Barutcu, A., and Nong, A. (2023). Assessing the basal gene expression of cancer cell lines for in vitro transcriptomic toxicology screening. bioRxiv.

Google Scholar

Black, M. B., Stern, A., Efremenko, A., Mallick, P., Moreau, M., Hartman, J. K., et al. (2022). Biological system considerations for application of toxicogenomics in next-generation risk assessment and predictive toxicology. Toxicol. Vitro 80, 105311. doi:10.1016/j.tiv.2022.105311

CrossRef Full Text | Google Scholar

Blencowe, B. J. (2006). Alternative splicing: new insights from global analyses. Cell 126, 37–47. doi:10.1016/j.cell.2006.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Braunschweig, U., Barbosa-Morais, N. L., Pan, Q., Nachman, E. N., Alipanahi, B., Gonatopoulos-Pournatzis, T., et al. (2014). Widespread intron retention in mammals functionally tunes transcriptomes. Genome Res. 24, 1774–1786. doi:10.1101/gr.177790.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Chepelev, N. L., Moffat, I. D., Labib, S., Bourdon-Lacombe, J., Kuo, B., Buick, J. K., et al. (2015). Integrating toxicogenomics into human health risk assessment: lessons learned from the benzo[a]pyrene case study. Crit. Rev. Toxicol. 45, 44–52. doi:10.3109/10408444.2014.973935

PubMed Abstract | CrossRef Full Text | Google Scholar

Colović, M. B., Krstić, D. Z., Lazarević-Pašti, T. D., Bondžić, A. M., and Vasić, V. M. (2013). Acetylcholinesterase inhibitors: pharmacology and toxicology. Curr. Neuropharmacol. 11, 315–335. doi:10.2174/1570159X11311030006

PubMed Abstract | CrossRef Full Text | Google Scholar

EFSA Scientific Committee Hardy, A., Benford, D., Halldorsson, T., Jeger, M. J., Knutsen, K. H., et al. (2017). Update: use of the benchmark dose approach in risk assessment. EFSA J. 15, e04658. doi:10.2903/j.efsa.2017.4658

PubMed Abstract | CrossRef Full Text | Google Scholar

Farmahin, R., Williams, A., Kuo, B., Chepelev, N. L., Thomas, R. S., Barton-Maclaren, T. S., et al. (2017). Recommended approaches in the application of toxicogenomics to derive points of departure for chemical risk assessment. Arch. Toxicol. 91, 2045–2065. doi:10.1007/s00204-016-1886-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Franzosa, J. A., Bonzo, J. A., Jack, J., Baker, N. C., Kothiya, P., Witek, R. P., et al. (2021). High-throughput toxicogenomic screening of chemicals in the environment using metabolically competent hepatic cell cultures. NPJ Syst. Biol. Appl. 7, 7. doi:10.1038/s41540-020-00166-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadar, S., Meller, A., Saida, N., and Shalgi, R. (2022). Stress-induced transcriptional readthrough into neighboring genes is linked to intron retention. iScience 25, 105543. doi:10.1016/j.isci.2022.105543

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrill, J., Shah, I., Setzer, R. W., Haggard, D., Auerbach, S., Judson, R., et al. (2019). Considerations for strategic use of high-throughput transcriptomics chemical screening data in regulatory decisions. Curr. Opin. Toxicol. 15, 64–75. doi:10.1016/j.cotox.2019.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrill, J. A., Everett, L. J., Haggard, D. E., Sheffield, T., Bundy, J. L., Willis, C. M., et al. (2021). High-throughput transcriptomics platform for screening environmental chemicals. Toxicol. Sci. 181, 68–89. doi:10.1093/toxsci/kfab009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hukkanen, J., Lassila, A., Päivärinta, K., Valanne, S., Sarpo, S., Hakkola, J., et al. (2000). Induction and regulation of xenobiotic-metabolizing cytochrome P450s in the human A549 lung adenocarcinoma cell line. Am. J. Respir. Cell Mol. Biol. 22, 360–366. doi:10.1165/ajrcmb.22.3.3845

PubMed Abstract | CrossRef Full Text | Google Scholar

Hulsen, T. (2022). DeepVenn -- a web application for the creation of area-proportional Venn diagrams using the deep learning framework Tensorflow.js. Available at: http://arxiv.org/abs/2210.04597 (Accessed June 30, 2023).

Google Scholar

Jacob, A. G., and Smith, C. W. J. (2017). Intron retention as a component of regulated gene expression programs. Hum. Genet. 136, 1043–1057. doi:10.1007/s00439-017-1791-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jensen, S. M., Kluxen, F. M., and Ritz, C. (2019). A review of recent advances in benchmark dose methodology. Risk Anal. 39, 2295–2315. doi:10.1111/risa.13324

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, K. J., Auerbach, S. S., and Costa, E. (2020). A rat liver transcriptomic point of departure predicts a prospective liver or non-liver apical point of departure. Toxicol. Sci. 176, 86–102. doi:10.1093/toxsci/kfaa062

PubMed Abstract | CrossRef Full Text | Google Scholar

Kappus, H. (1987). Oxidative stress in chemical toxicity. Arch. Toxicol. 60, 144–149. doi:10.1007/BF00296968

PubMed Abstract | CrossRef Full Text | Google Scholar

Kolberg, L., Raudvere, U., Kuzmin, I., Adler, P., Vilo, J., and Peterson, H. (2023). g:Profiler-interoperable web service for functional enrichment analysis and gene identifier mapping (2023 update). Nucleic Acids Res. 51, W207–W212. doi:10.1093/nar/gkad347

PubMed Abstract | CrossRef Full Text | Google Scholar

Krebs, A., Nyffeler, J., Karreman, C., Schmidt, B. Z., Kappenberg, F., Mellert, J., et al. (2020). Determination of benchmark concentrations and their statistical uncertainty for cytotoxicity test data and functional in vitro assays. ALTEX 37, 155–163. doi:10.14573/altex.1912021

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y., and Rio, D. C. (2015). Mechanisms and regulation of alternative pre-mRNA splicing. Annu. Rev. Biochem. 84, 291–323. doi:10.1146/annurev-biochem-060614-034316

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi:10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Fan, L., Lu, C., Yin, S., and Hu, H. (2020). Functional role of p53 in the regulation of chemical-induced oxidative stress. Oxid. Med. Cell. Longev. 2020, 6039769. doi:10.1155/2020/6039769

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Yang, H., and Liu, Z. (2022). Signaling pathways involved in paraquat-induced pulmonary toxicity: molecular mechanisms and potential therapeutic drugs. Int. Immunopharmacol. 113, 109301. doi:10.1016/j.intimp.2022.109301

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550–621. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Lukaszewicz, A., Cwiklinska, M., Zarzecki, M., Szoka, P., Lachowicz, J., and Holownia, A. (2019). Cytotoxicity, oxidative stress, and autophagy in human alveolar epithelial cell line (A549 cells) exposed to standardized urban dust. Adv. Exp. Med. Biol. 1176, 101–108. doi:10.1007/5584_2019_387

PubMed Abstract | CrossRef Full Text | Google Scholar

Matteo, G., Leingartner, K., Rowan-Carroll, A., Meier, M., Williams, A., Beal, M. A., et al. (2023). In vitro transcriptomic analyses reveal pathway perturbations, estrogenic activities, and potencies of data-poor BPA alternative chemicals. Toxicol. Sci. 191, 266–275. doi:10.1093/toxsci/kfac127

PubMed Abstract | CrossRef Full Text | Google Scholar

Mittal, M., Siddiqui, M. R., Tran, K., Reddy, S. P., and Malik, A. B. (2014). Reactive oxygen species in inflammation and tissue injury. Antioxid. Redox Signal. 20, 1126–1167. doi:10.1089/ars.2012.5149

PubMed Abstract | CrossRef Full Text | Google Scholar

Monteuuis, G., Wong, J. J. L., Bailey, C. G., Schmitz, U., and Rasko, J. E. J. (2019). The changing paradigm of intron retention: regulation, ramifications and recipes. Nucleic Acids Res. 47, 11497–11513. doi:10.1093/nar/gkz1068

PubMed Abstract | CrossRef Full Text | Google Scholar

National Toxicology Program (2018) NTP research report on national toxicology program approach to genomic dose-response modeling: research report 5. Durham (NC): National Toxicology Program.

Google Scholar

Niechoda, A., Roslan, J., Maciorowska, K., Rosłan, M., Ejsmont, K., and Holownia, A. (2023). Oxidative stress and activation of H2A.X in lung alveolar epithelial cells (A549) by nanoparticulate carbon black. Respir. Physiol. Neurobiol. 316, 104140. doi:10.1016/j.resp.2023.104140

PubMed Abstract | CrossRef Full Text | Google Scholar

Pagé-Larivière, F., Crump, D., and O’Brien, J. M. (2019). Transcriptomic points-of-departure from short-term exposure studies are protective of chronic effects for fish exposed to estrogenic chemicals. Toxicol. Appl. Pharmacol. 378, 114634. doi:10.1016/j.taap.2019.114634

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, E., Pan, Z., Zhang, Z., Lin, L., and Xing, Y. (2018). The expanding landscape of alternative splicing variation in human populations. Am. J. Hum. Genet. 102, 11–26. doi:10.1016/j.ajhg.2017.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, J. R., Svoboda, D. L., Tandon, A., Patel, S., Sedykh, A., Mav, D., et al. (2019). BMDExpress 2: enhanced transcriptomic dose-response analysis workflow. Bioinformatics 35, 1780–1782. doi:10.1093/bioinformatics/bty878

PubMed Abstract | CrossRef Full Text | Google Scholar

Pouzou, J. G., Kissel, J., Yost, M. G., Fenske, R. A., and Cullen, A. C. (2020). Use of benchmark dose models in risk assessment for occupational handlers of eight pesticides used in pome fruit production. Regul. Toxicol. Pharmacol. 110, 104504. doi:10.1016/j.yrtph.2019.104504

PubMed Abstract | CrossRef Full Text | Google Scholar

Rothfuss, O., Gasser, T., and Patenge, N. (2010). Analysis of differential DNA damage in the mitochondrial genome employing a semi-long run real-time PCR approach. Nucleic Acids Res. 38, e24. doi:10.1093/nar/gkp1082

PubMed Abstract | CrossRef Full Text | Google Scholar

Scotti, M. M., and Swanson, M. S. (2015). RNA mis-splicing in disease. Nat. Rev. Genet. 17, 19–32. doi:10.1038/nrg.2015.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Tapial, J., Ha, K. C. H., Sterne-Weiler, T., Gohr, A., Braunschweig, U., Hermoso-Pulido, A., et al. (2017). An atlas of alternative splicing profiles and functional associations reveals new regulatory programs and genes that simultaneously express multiple major isoforms. Genome Res. 27, 1759–1768. doi:10.1101/gr.220962.117

PubMed Abstract | CrossRef Full Text | Google Scholar

Ule, J., and Blencowe, B. J. (2019). Alternative splicing regulatory networks: functions, mechanisms, and evolution. Mol. Cell 76, 329–345. doi:10.1016/j.molcel.2019.09.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Upadhyay, S., Vaish, S., and Dhiman, M. (2019). Hydrogen peroxide-induced oxidative stress and its impact on innate immune responses in lung carcinoma A549 cells. Mol. Cell. Biochem. 450, 135–147. doi:10.1007/s11010-018-3380-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Vandenberg, L. N., Colborn, T., Hayes, T. B., Heindel, J. J., Jacobs, D. R., Lee, D.-H., et al. (2012). Hormones and endocrine-disrupting chemicals: low-dose effects and nonmonotonic dose responses. Endocr. Rev. 33, 378–455. doi:10.1210/er.2011-1050

PubMed Abstract | CrossRef Full Text | Google Scholar

Villaseñor-Altamirano, A. B., Watson, J. D., Prokopec, S. D., Yao, C. Q., Boutros, P. C., Pohjanvirta, R., et al. (2019). 2,3,7,8-Tetrachlorodibenzo-p-dioxin modifies alternative splicing in mouse liver. PLoS One 14, e0219747. doi:10.1371/journal.pone.0219747

PubMed Abstract | CrossRef Full Text | Google Scholar

Vu, T. D., Ito, N., Oshima, K., Maruko, A., Nishi, A., Mizoguchi, K., et al. (2022). Intron retention is a stress response in sensor genes and is restored by Japanese herbal medicines: a basis for future clinical applications. Gene 830, 146496. doi:10.1016/j.gene.2022.146496

PubMed Abstract | CrossRef Full Text | Google Scholar

Wong, J. J.-L., Au, A. Y. M., Ritchie, W., and Rasko, J. E. J. (2016). Intron retention in mRNA: No longer nonsense: known and putative roles of intron retention in normal and disease biology. Bioessays 38, 41–49. doi:10.1002/bies.201500117

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaharieva, E., Chipman, J. K., and Soller, M. (2012). Alternative splicing interference by xenobiotics. Toxicology 296, 1–12. doi:10.1016/j.tox.2012.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Park, C., Bennett, C., Thornton, M., and Kim, D. (2021). Rapid and accurate alignment of nucleotide conversion sequencing reads with HISAT-3N. Genome Res. 31, 1290–1295. doi:10.1101/gr.275193.120

PubMed Abstract | CrossRef Full Text | Google Scholar