Validation of alternative splicing in the xenobiotic-induced PRH model

Before initiating a global splicing analysis, we aimed to validate our in vitro model system by confirming the expression and splicing of known in vivo Cyp2b2 variants induced by xenobiotics. This step was essential for establishing physiological relevance and ensuring our model’s reliability for downstream transcriptome-wide splicing studies. In this study, we focused on recapitulating the xenobiotic-induced expression and splicing patterns of the Cyp2b gene family in vivo [18], specifically the Cyp2b2v splice variant. This focus serves as a quality control measure for our primary rat hepatocyte (PRH) model, providing physiological relevance to our global splicing analysis. The Cyp2b2v variant includes a small addition of eight amino acids (VSPAWMRE), resulting from a cryptic 5’-splice site in intron 5. Detailed structural and sequence information regarding the Cyp2b2 gene and CYP2B2v protein is available in the Additional Materials (Supplemental Fig. 1). To achieve our goals, we compared the bioactivity of three discrete xenobiotics (Fig. 1), focusing on their targeted effects on Cyp2b2v splicing and broader impacts on the global transcriptome. This analysis included comparisons of two MDB compounds (C1-MDB and C6-MDB) and carbamazepine (CBZ), each known to be mixed inducers of multiple CYP enzymes from CYP Families 1–3 [18, 19, 59].

Fig. 1

Xenobiotics used to Probe Structure-Activity Relationships and Global Splicing Modulation in Primary Rat Hepatocytes. The figure illustrates the three xenobiotics investigated for their effects on alternative splicing across the primary rat hepatocyte (PRH) transcriptome. Two methylene-dioxybenzene (MDB) compounds, C1-MDB (4-n-methyl-methylene-dioxybenzene) with a short aliphatic chain and C6-MDB (4-n-hexyl-methylene-dioxybenzene) with an extended chain, were included to explore structure-activity relationships influencing alternative gene splicing and expression. These MDBs are known to induce differential splicing of Cyp2b genes, providing a basis to examine their broader splicing effects across other gene families. Additionally, carbamazepine (CBZ), an FDA-approved drug with mixed CYP induction properties, acts as a comparative standard by modulating the expression of several CYP isoforms overlapping with those influenced by the MDB compounds. This selection allows us to probe whether structural variations drive unique splicing patterns, despite similar changes in gene expression

Here, we developed a semiquantitative endpoint PCR assay for rat Cyp2b1, Cyp2b2, and Cyp2b2v, multiplexed with the SDHA housekeeping gene, to validate the basal expression level and inducibility of the Cyp2b2v splice variant transcript in PRHs grown in sandwich culture. Because some in vitro hepatocyte cultures fail to replicate in vivo expression or splicing [60], confirming Cyp2b2v splice variant expression served as a quality control for our global splicing analysis, as its absence could indicate broader issues with splicing machinery in our system. To detect this variant, primers were designed to target non-homologous regions of the Cyp2b2v mRNA, which differs from the wild-type by only 24 bp, resulting in distinct amplicon sizes of each gene (Supplemental Table 1; Additional Materials).

Test concentrations for the three xenobiotics (Fig. 1) and the vehicle (DMSO; 0.5-1%) were based on published studies [18] and PRH toxicity dose-response profiling (100–500 µM), as highlighted in Supplemental Figs. 2 and 3 (Additional Materials). Optimal treatment (100–250 µM) and vehicle control (0.05%) levels were selected based on pilot studies of fresh PRHs in a 24-well, sandwich culture format. A modest, but insignificant difference in Cyp2b2v mRNA was detected at 1% DMSO, leading us to select 0.5% as the optimal vehicle level for these studies (Supplemental Fig. 4). Multiplexed endpoint PCR analysis of Cyp2b1, Cyp2b2, and Cyp2b2v gene expression showed clear separation of amplified PCR bands: Cyp2b1 (110 bp), Cyp2b2 (394 bp), and Cyp2b2v (843 bp), from the SDHA housekeeping gene (218 bp). These distinct and diagnostic bands enable semi-quantitative comparisons of CYP expression across treatment groups using CYP/SDHA ratios. (Fig. 2a). Treatment of PRHs (24 h) with 250 µM C1-MDB C1 did not significantly alter basal Cyp2b1 mRNA expression levels, similar to vehicle (0.5% DMSO). However, treatment with 250 µM C6-MDB resulted in an approximately two-fold increase in Cyp2b1 transcript (*p < 0.05). Treatment with 100 µM CBZ resulted in an increase in Cyp2b1 message that was significantly increased from both vehicle (**p < 0.005) and C6 MDB (*p < 0.05) as shown in Fig. 2b.

Fig. 2

Endpoint PCR Analysis of Cyp2b1, Cyp2b2 and Cyp2b2v mRNA Expression in Fresh Sprague-Dawley PRHs. (a) Results of semiquantitative, endpoint PCR using fresh PRHs grown in sandwich culture, for rat CYPs 2b1, 2b2, and 2b2v. The succinate dehydrogenase (SDHA) gene was selected as the optimal housekeeping gene in this study, as compared to ACTB and GAPDH (not shown). (b) CYP2B1 mRNA levels were normalized to SDHA and quantified using Image J software analysis using 3 replicate images. Statistically significant changes are denoted (** indicates p < 0.01; Student T-test) and standard deviation is shown. (c) Cyp2b2 WT and (d) Cyp2b2v mRNA levels were also normalized and quantified in a similar way using Image J analysis. Statistically significant changes (* indicates p > 0.05, ** indicates p < 0.01, *** indicates p < 0.005) are denoted and standard deviation is shown for each treatment group

Expression levels of mRNA encoding rat Cyp2b2 were unchanged comparing untreated and vehicle treated (0.5% DMSO) PRHs. Treatment with 250 µM C1-MDB resulted in a significant, nearly 5-fold, increase in wild-type Cyp2b2 message (***p < 0.001) compared to 0.5% DMSO vehicle control. Treatment with 250 µM C6-MDB significantly increased Cyp2b2 transcript over vehicle (****p < 0.0001) nearly 10-fold, and over C1-MDB (*p < 0.05) over 2-fold. Similarly, 100 µM CBZ treatment induced a significant, near 15-fold, increase in wild-type Cyp2b2 mRNA over vehicle (****p < 0.0001), and an approximate 3-fold increase in C1 MDB treatments (*p < 0.05) as shown in Fig. 2c.

With respect to the Cyp2b2v splice variant, vehicle treated (0.5% DMSO) and untreated PRHs did not differ significantly in Cyp2b2v transcript level. Treatment with 250 µM C1-MDB resulted in a significant ~ 3-fold increase in Cyp2b2v message (**p < 0.005) compared to vehicle control. Both 250 µM C6-MDB and 100 µM-CBZ treatment induced Cyp2b2v mRNA by over 10-fold compared to the vehicle control (****p < 0.0001) and over 3-fold compared to 250 µM C1 MDB (**p < 0.005) (Fig. 2d).

The mRNA expression ratios of wild-type Cyp2b1 to Cyp2b2v and wild-type Cyp2b2 to Cyp2b2v were quantified using endpoint PCR. As detailed in Supplemental Fig. 5 (Additional Materials), Cyp2b2v expression was consistently high, reaching approximately 90 to 100% of wild-type Cyp2b2 levels across all controls (Untreated, 0.88 ± 0.46; Vehicle, 0.94 ± 0.46) and most treatment groups (C1-MDB, 0.92 ± 0.45; CBZ, 0.98 ± 0.50). However, C6-MDB exposure resulted in a two-fold reduction in Cyp2b2v levels (0.55 ± 0.44). The Cyp2b1 to Cyp2b2v ratio, initially 12–15 times higher in controls (Untreated, 12.71 ± 2.83; Vehicle, 15.55 ± 3.91), decreased substantially upon treatment: C1-MDB (3.76 ± 1.30), CBZ (3.49 ± 1.08), and most notably C6-MDB (1.62 ± 0.17). This reduction indicates that although both Cyp2b1 and Cyp2b2 may be induced by xenobiotics, Cyp2b2 induction outpaces Cyp2b1. The pronounced decrease in the Cyp2b2v variant ratio to both Cyp2b1 and Cyp2b2 with C6-MDB suggests possible the involvement of a specific xenosensor (e.g., AHR, ) whose activation may function to alter Cyp2b gene transcript stability in a unique way. These findings underscore the complex regulatory interactions that modulate global splice variant production during xenobiotic exposures, and these results provide a reference point for comparing splicing ratios in RNA-seq analyses.

Global differential gene expression profiling of xenobiotic-induced PRHs using RNAseq

After validating the expression of CYPs 2B1, 2B2 and 2B2v in our PRH sandwich culture model, we conducted an Illumina-based RNAseq Analysis of global gene expression and splicing changes in PRHs treated with 250 µM C1-MDB, C6-MDB and 100 µM CBZ. Each treatment group (untreated, vehicle (0.05% and 1%), C1-MDB, C6-MDB, and CBZ) was sequenced in quadruplicate. Principal component analysis (PCA) of each biological replicate and treatment group showed strong correlation with limited outliers (Fig. 3a). A Pearson correlation plot further confirmed this, with high correlation (R² > 0.9) across all samples, except for one outlier (C6-MDB-2) shown in Supplemental Fig. 6a (Additional Materials). A Venn Diagram generated using the Comparative Toxicogenomics Database (CTD) [61] illustrates the full spectrum of differentially expressed genes across all treatments in our RNAseq analysis, including genes unique to each compound and those shared among each pair (Fig. 3b). We identified over 10,000 genes active in our four primary treatment groups (Vehicle (0.5% DMSO), C1-MDB, C6-MDB, and CBZ), as highlighted at the center of the Venn diagram, with each treatment modulating ~ 50 to 150 unique genes.

Fig. 3

Global RNAseq Analysis of Differential Gene Expression in C1-MDB, C6-MDB, and CBZ-treated PRHs. (a) Principal Component Analysis (PCA) of RNA-seq data from untreated, control (VEH_05 and VEH_1), and treated (C1-MDB, C6-MDB, and CBZ) primary rat hepatocyte groups, showing strong clustering across groups. Notable outliers include C6-MDB replicate 2 along PC1 and variance in C1-MDB along PC2. To maintain a comprehensive dataset, all replicates were included in downstream analyses to capture the full scope of gene expression and splicing changes in the PRH model. (b) Venn diagram depicting the overlap of differentially expressed genes (DEGs), relative to control (VEH_05), across treatments (C1-MDB, C6-MDB, and CBZ) based on DESeq2 analysis. The central overlap represents over 10,000 commonly expressed genes, though not all are differentially regulated. Each xenobiotic modulated approximately 50 to 150 unique genes, indicating treatment-specific regulation, while the central set reflects genes expressed under all conditions. (c) Venn diagram from DEG analysis using rSeqDiff, identifying 58 genes commonly altered across all treatments (C1-MDB, C6-MDB, CBZ) compared to control. Unlike the DESeq2 analysis, the control (VEH_05) is not shown, focusing instead on shared gene changes between treatments. C6-MDB and C1-MDB showed more substantial gene expression changes (488 and 368 genes, respectively) compared to CBZ (182). CBZ and C1-MDB shared the highest number of commonly modulated genes (136), compared to CBZ and C6-MDB (88) or C1-MDB and C6-MDB (92). Detailed control comparisons are provided in Supplemental Tables 3–6 (Additional Materials)

Differentially expressed gene (DEG) analysis using rSeqDiff [62] and DESeq2 [63] produced nearly identical results. In the DESeq2 dataset, 58 genes were commonly differentially expressed (log2 fold change > ± 1) across all three treatments (C1-MDB, C6-MDB, and CBZ) compared to control (0.5% DMSO). C6-MDB induced the largest expression changes (488 genes), followed by C1-MDB (368 genes) and CBZ (182 genes). Notably, CBZ and C1-MDB shared the most commonly modulated genes (134), compared to CBZ and C6-MDB (88) or C1- and C6-MDB (91) (Fig. 3c). GO Biological Process enrichment analysis of the 58 common DEGs highlighted roles in oxidative stress regulation, cytokine production, and inflammation, with key genes like FMO1, NOS2, HSPA1B, NOX4, and DUOX1. These findings suggest coordinated regulation of NADPH oxidase activity and cytokine signaling, which are crucial defenses against xenobiotics and infections (Supplemental Fig. 6b; Additional Materials). Using the Transcription Factor (PPI) pathway analysis module on the Enrichr bioinformatics platform [64,65,66], which leverages a literature-based protein interaction map, we identified a network of transcription factors (IRF1, IRF6, NFKB1, MYC, HSF1, NR2F2 (COUP-TF2)) driving gene expression changes in response to xenobiotic exposure. Each transcription factor in this module is associated with a set of interacting proteins, illustrating potential regulatory connections within the network (Supplemental Fig. 6c; Additional Materials).

DESeq2 analysis identified the most differentially expressed genes across the three xenobiotic treatments, detailed in Supplemental Table 2 (Additional Materials). C1-MDB strongly induced genes (up to 220-fold; Adgre1) linked to detoxification pathways (Ugt2b7), immune responses (Adgre1, Ly6al), and tissue repair mechanisms (Lrg1, Nos2), while suppressing genes (up to 943-fold; Rps9) linked to protein synthesis (Rps9), growth regulation (Igfbp6), neuronal development (Nxpe3), immune signaling (Cd209), and potassium ion transport (Kcnk1). C6-MDB significantly induced genes (up to 136-fold; Cyp2b1) involved in xenobiotic metabolism (Cyp2b1, Cyp2b2, Cyp1a1, Cyp3a23-3a1) and RNA processing (Rpph1) and suppressed those related to DNA binding and transcription (Zkscan8, up to 22-fold), cell adhesion (Lamb2), non-coding RNA regulation (Thumpd3-as1), and receptor tyrosine kinase signaling (Ephb6). CBZ, more similar to C6-MDB than C1-MDB overall, significantly induced genes (up to 71-fold; CYP2B1) involved in xenobiotic metabolism (Cyp2b1, Cyp2b2, Cyp2c6v1), rRNA processing (Rn28s), and immune modulation (Lexm) while suppressing genes (up to ~ 1600-fold; Rps9) linked to steroid hormone metabolism (Hsd17b1), protein synthesis (Rps9), cytoskeleton regulation (Rac2), leucine-rich repeat containing (Lrrc36), and transcriptional repression (Hopx).

Exposure to all three xenobiotics significantly induced genes involved in xenobiotic metabolism, immune responses, and tissue repair, while suppressing genes related to protein synthesis, most notably Rps9. This suppression of Rps9 may indicate a strategic cellular shift under xenobiotic stress, prioritizing detoxification, damage control, and survival pathways as part of the defensome response [67]. KEGG Pathway analysis of all differentially expressed genes induced (Supplemental Table 1; Additional Materials) or suppressed (Supplemental Table 4; Additional Materials) by all 3 treatments (C1-MDB, C6-MDB and CBZ) confirmed these trends, with the pathway for “steroid hormone biosynthesis” occurring in both lists, highlighting its complex association with xenobiotic exposures like endocrine disrupting chemicals (EDCs). GO Biological Process enrichment (Supplemental Table 5; Additional Materials) and KEGG Pathway analysis (Supplemental Table 6; Additional Materials) were also conducted on the 58 common genes modulated by all 3 xenobiotics, revealing significant involvement in crucial cellular response pathways, including the superoxide metabolic process (DUOX1, NOS2, NOX4) and the inflammatory response (CXCL9, CEBPB, NOS2, NOX4), which were highly enriched to mitigate oxidative stress and inflammation. Furthermore, the regulation of cytokine production (HIC2, NOS2, IL15, IRF8) and the positive regulation of wound healing (DUOX1, SMOC2) highlight the complex interplay between immune modulation and tissue repair. KEGG pathway analysis also pointed to significant modulation within the Pertussis (NOS2, IRF8), IL-17 signaling pathway (CEBPB, FOSB), and TNF signaling pathway (CEBPB, IL15), illustrating the xenobiotics’ impact on immune signaling, with potential implications for host defense mechanisms.

Detailed KEGG pathway analysis of differentially expressed genes (DEGs) induced or suppressed by C1-MDB, C6-MDB, or CBZ, individually, reveals distinct biological impacts of these xenobiotics on pathways linked to metabolic detoxification, cytokine signaling and the immune response. (Supplemental Tables 7–12; Additional Materials). These findings are expanded upon in the Supplemental Results (Additional Materials), providing a deeper discussion on the specific genes and pathways modulated by each xenobiotic. In summary, C1-MDB’s modulation of gene expression centered on inflammatory responses and chemotaxis, potentially enhancing immune defense mechanisms. C6-MDB targeted cholesterol metabolism and steroid biosynthesis, suggesting impacts on lipid homeostasis and hormonal processes. CBZ modulated genes linked to cholesterol catabolism and wound response.

A comprehensive GO Biological Process Analysis of all significant (p < 0.05) DEGs modulated by C1-MDB, C6-MDB, and CBZ revealed distinct biological effects of each xenobiotic (Supplemental Tables 13–15; Additional Materials). The GO Biological Process findings analysis are further detailed in the Supplemental Results (Additional Materials), but overall, enriched GO terms aligned with the KEGG analysis, highlighting specific effects of each xenobiotic on immune response, metabolic processes, and tissue homeostasis. The distinct cellular effects of C1-MDB, C6-MDB, and CBZ stem from their roles as mixed inducers of metabolic genes, each activating a distinct array of xenosensors and signaling networks that drive unique transcriptional responses in inflammation, metabolism and tissue repair.

Volcano Plot Analysis of metabolic gene modulation in primary rat hepatocytes treated with vehicle control (0.5% DMSO), C1-MDB (250 \(\:\mu\:\)M), C6-MDB (250 \(\:\mu\:\)M), and CBZ (100 \(\:\mu\:\)M) was performed using VolcaNoseR [68]. To provide a clear definition of “metabolic genes,” the primary Phase 0-III genes involved in xenobiotic disposition, as identified in our analysis, are listed in Supplemental Table 16 (Additional Materials). Differential expression profiles of significantly modulated (log2fold change > ± 1.5; p-value < 0.05) Phase 0 (NR genes, AHR related genes), Phase I (CYPs, FMO, etc.), Phase II (UGTs, GSTs, SULTs, etc.) and Phase III (SLC and ABC family transporters) genes were analyzed. Minimal transcriptional changes were observed for the vehicle control, with notable suppression of phase II (Gstm5) and phase III (Slc25a18) genes and induction of one Phase 0 (vitamin D receptor; Vdr) and one Phase III (Slc26a10) gene (Fig. 4a). In C1-MBD treated cells, differential expression of several Phase I-III metabolic enzymes was observed (Fig. 4b). This included significant induction of five Phase I enzymes, including CYPs 1a2, 2b1, 2b2, 2c6v1 and 4f37, along with one Phase II enzyme (Gsto2) and one phase III transporter (Slc26a24). Significant suppression was noted for the phase III transporter Slc22a13, as well as for Gstm5 and Slc25a18, which were also suppressed by the vehicle control (Fig. 4b). In contrast, C6-MDB treated primary rat hepatocytes resulted in stronger modulation of metabolic enzymes compared to C1-MDB. Differential expression of 13 CYP genes were detected in C6-MDB-treated cells, with notable induction of CYPs 1a1, 1a2, 2b1, 2b2, 2c6v1, 2c13, 3a23/3a1, 3a2, and 7a1, and suppression of CYPs 7b1, 8b1 and 46a1 (Fig. 4c). C6-MDB also induced significant differential expression in five phase II genes, with induction of Gsta2, Gsto2, Gsta5 and Gstt3, and suppression of Gstm5, which was suppressed in the control and all treatments except CBZ (Fig. 4c). C6-MBD treatments also strongly modulated phase III transporters and efflux pumps, inducing both ABC (Abcd2) and SLC family transporters (Slc4a1, Slc5a10, Slc16a5, Slc25a45, and Slco1a4) and suppressing one ABC transporter (Abc8a) and ten SLC transporters (Slc1a2, Slc1a5, Slc2a4, Slc2a9, Slc11a1, Slc13a5, Slc25a18, Slc34a2, Slc51a) (Fig. 4c). CBZ induced a similar set of phase I genes as the C6-MDB, with CYPs 2b1, 2b2, 2c6v1, 2c11, 2c13, 3a2, 3a23/3a1, and 7a1 (Fig. 4d). However, CBZ did not induce CYPs from the 1 A family (Cyps 1a1, 1a2). Additionally, CBZ suppressed several CYPs, including Cyp3a9 and the cholesterol metabolizing enzymes Cyp7b1 and Cyp8b1. Although Cyp8b1 fell slightly below our significance threshold (log2fc = -1.45) its concurrent suppression with Cyp7b1 mimics the response seen with C6-MDB, but not C1-MBD. In contrast, CBZ modulation of phase II enzymes differed from C6-MDB, with significant induction of Gsta5, Gsto2 and Ugt2b17, and no suppression of Gstm5. Unique to CBZ was the suppression of the phase III transporter Abcg1. CBZ also modulated 10 SLC transporters, inducing Slc5a10, Slc16a5, Slc25a45, and Slc26a5, while suppressing Slc2a6, Slc12a5, Slc13a2, Slc25a18, Slc26a9, and Slco1a2, showing only partial overlap with the phase III modulation by C6-MDB (Fig. 4d).

Fig. 4

Volcano Plot Analysis of Metabolic Gene Modulation in Vehicle Control, C1-MDB, C6-MDB and CBZ-treated PRHs. The differential expression of important Phase 0 (NR genes), Phase I (CYP genes), Phase II (UGTs, GSTs, SULTs, etc.) and Phase III (SLC and ABC transporters) genes were analyzed across each treatment group and the vehicle control. (a) Compared to untreated PRHs, cells treated with 0.5% (v/v) DMSO (vehicle) showed minimal induction of metabolic enzymes, but we detected a modest, but significant upregulation (log2fold change > ± 1.5; p-value < 0.05) of the vitamin D receptor (Vdr; NR1i1). The vehicle control also significantly modulated the expression of one phase II enzyme (Gstm5) and two, phase III transporters (Slc26a10 and Slc25a18) (b) For C1-MBD treatments, significant differential expression of phase I CYPs (1a2, 2b1, 2b2, 2c6v1, and 4f37 (all upregulated)), Phase II conjugating enzymes (Gsta5 (downregulated), Gsto2 (upregulated), and Phase III transporters (Slc22a13 and Slc25a18; downregulated; Slc26a10 upregulated) were detected. (c) In contrast, C6-MDB induced stronger modulation of multiple metabolic enzymes. Differential expression of 13 CYP genes were detected for C6-MDB, including strong induction of CYPs 1a1, 1a2, 2b1, 2b2, 2c6v1, 2c11, 2c13, 3a23/3a1, and 7a1, along with suppression of CYPs 7b1, 8b1 and 46A1. C6-MDB also induced significant differential expression in several phase II genes, including Gsta2, Gsta5, Gsto2, Gstt3 (upregulated) and Gstm5 (downregulate), and multiple phase III transporter genes, some of which were strongly induced (e.g. Slc16a5), while others were significantly suppressed (e.g. Slc1a2 and Slc1a5). (d) CBZ induced a similar set of phase I genes as C6-MDB, with CYPs 2b1, 2b2, 2c6v1 and 3a2 being most strongly induced. Several CYPs were also suppressed by CBZ, including CYPs 7a1, 8b1, and Cyp3a9. CBZ also strongly induced several phase II enzymes, including Gsta5, Gsto2, and Ugt2b17, without suppressing any. CBZ also had major effects on phase III enzyme expression, inducing or suppressing 10 distinct SLC and ABC transporter genes

RNAseq analysis of alternative gene splicing in xenobiotic-induced PRHs

To quantify alternative splicing across all genes in our PRH sandwich culture model, we developed a bioinformatics pipeline incorporating rSeqDiff [62], rMATs-turbo (rMATs 4.1.2; [69, 70]) and the Integrated Genomic Viewer (IGV; [71,72,73,74]). Although alternative splicing changes were detected using rSeqDiff (model 2; T-value), this tool does not generate a false discovery rate (FDR), essential for high-confidence quantification of splicing events. The absence of FDR values can lead to false positives, where statistical anomalies might be misinterpreted as genuine biological changes. To address similar challenges, Luo et al. (2020) propose a two-step mixed model approach to improve the reliability of differential alternative RNA splicing detection by incorporating a FDR estimation method [75]. Despite these limitations, we compared rSeqDiff T-values (model 2; splicing changes) with differential gene expression values (model 1; log2fold-Δ data) to visualize global trends between gene expression and alternative splicing among all 3 treatment groups and vehicle (0.5% DMSO) (see Supplemental Fig. 7a-7d; Additional Materials). Based on rSeqDiff analysis only, C1-MDB exposure was associated with novel splicing events in genes with both increased and decreased expression, while C6-MDB and CBZ treatments showed more pronounced alternative splicing linked to gene induction than suppression.

To address the limitations of rSeqDiff, we used rMATs-turbo for our global gene splicing analysis, applying a 5% false discovery rate (FDR) as a threshold for significance. rMATs-turbo detects five types of alternative splicing events in RNAseq data: alternate 3’ spice site (A3SS), alternative 5’ splice sites (A5SS), skipped exon (SE), mutually exclusive exon (MXE), and retained intron (RI). The difference in PSI (percent splice in) values between the control (vehicle) and each xenobiotic treatment reflects changes in the frequency of individual splicing events. A negative PSI difference indicates increased alternative splicing events in the treatment group relative to the control.

We analyzed global RNA splicing changes using the same RNAseq dataset as the DEG analysis, with cumulative rMATs-turbo results shown in Fig. 5a and b. In PRHs treated with 250 µM C1-MDB, 68,554 total alternative splicing events were detected, of which 6,525 were significant (p < 0.05). Approximately 40% were skipped exons (SE) (27,056 total; 3,112 significant), 19% were A5SS (13,071 total; 930 significant), 29% were A3SS (19,832 total; 1,557 significant), 5% were MXE (3,681 total; 398 significant), and 7% were RI (4,914 total; 28 significant). Similar trends appeared for PRHs treated with 250 µM C6-MDB, with 66,056 total alternative splicing events being detected, with 4,945 significant events. Of these 37% were SE (24,616 total; 2,274 significant), 21% were A5SS (13,758 total; 731 significant), 29% were A3SS (19,028 total; 1,133 significant), 6% were MXE (3,792 total; 352 significant), and 7% RI (4,862 total; 455 significant). The 100 µM CBZ treatment produced the highest number of splicing events (77,035), although fewer were significant (5,906) compared to C1-MDB (6,525). Among these, 39% were SE (30,167 total; 2,773 significant), 20% were A5SS (15,010 total; 838 significant), 29% were A3SS (22,058 total; 1,265 significant), 6% were MXE (4,316 total; 352 significant), and 7% were RI (5,484 total; 455 significant).

Fig. 5

rMATs analysis of Global Differential Splicing Events in Cryopreserved PRHs. A multivariate analysis of alternative splicing in within our RNA-seq data for MDB and CBZ treated PRHs was completed with the program rMATS 4.1.2. A) We detected over 66,000 individual splicing events in all 3 of our experimental treatments, with C1-MDB, C6-MDB, and CBZ. CBZ-treated cells induced over 77,000 splicing events, as compared to the vehicle control. B) There were fewer significant (p < 0.05) splicing events in total (Average = 5792 unique events across the 3 treatments), with C1-MDB (6525) treated PRHs, expressing slightly more variants than CBZ (5906) or C6-MDB (4945) treated cells. The most common events across all treatments were skipped exons (SE; Average = 2720 events per treatment), followed by alternative 3’ splice site selection (A3SS; Average = 1318 events per treatment), alternative 5’ splice site selection (A5SS; Average = 833 events per treatment), retained introns (RI; 508 events per treatment) and finally, mutually exclusive exon usage (MXE; Average = 413 events per treatment)

An overview of the most commonly spliced genes (in terms of total events) and the most highly spliced genes (in terms of the largest ΔPSI) are summarized for each of the 3 treatments (C1-MDB, C6-MDB and CBZ) in Supplemental Tables 17–22 (Additional Materials). We observed significant overlaps, along with some unique responses, in alternative splicing across the 3 treatments, with albumin (alb) being the most commonly spliced gene among each of the 3 treatments (31 individual splicing events; C1-MDB) and cumulatively (53 total splicing events; across all 3 treatment groups). The splicing of alb and Orm1 are notably modulated across all three treatments, indicating their critical role in response to xenobiotic stress, with alb undergoing predominantly MXE and SE events, and Orm1 experiencing a greater variety of splicing changes. Orm1 is part of the lipocalin family and acts as a transport protein in the blood, similar to albumin. NXPE4, a gene involved in neural development, potentially as a receptor or adhesion molecule, was found to be the most commonly spliced gene after C6-MDB treatment (15 events; SE), and was not detected in the other treatment groups. Furthermore, genes like Gstm3 and Fxyd1 were affected by multiple treatments but displayed unique splicing events for each of the 3 exposures, highlighting the complexity and dynamic nature of the splicing process.

GO Biological Process enrichment analysis of 833 genes modulated (> 10% change in a splicing event) by the C1-MDB treatment (Supplemental Table 23 Additional Materials) reveals significant modulation in biological processes related to RNA dynamics and chromatin remodeling, including RNA- and miRNA-processing, mRNA stability, and chromatin remodeling, highlighting a widespread impact on the regulation of gene expression. Overlap genes in these categories, such as AGO4, TUT4, and PUM1 are involved in RNA silencing, mRNA decay, and post-transcriptional modulations, underscoring the complex regulation of RNA-binding proteins and their modifiers in response to xenobiotic stress.

GO Biological Process enrichment analysis of the 862 genes influenced (x > 10% in PSI) by the C6-MDB treatment (Supplemental Table 24; Additional Information) identified significant alterations in processes related to RNA splicing, processing, and metabolic regulation, showcasing the treatment’s profound impact on the RNA lifecycle. The inclusion of processes related to monocarboxylic acid biosynthesis and vesicle-mediated transport between the ER and Golgi apparatus indicates a broader cellular response, affecting a more complex matrix of metabolic pathways and intracellular trafficking events than the C1-MDB treatment. Overlap genes like SRPK2, MBNL1, and HNRNPH1, which regulate discrete aspects of splicing, underscore the intricate regulation of RNA-binding proteins and spliceosomal components in response to xenobiotic stress. It is notable that C1-MDB had a significant effect on miRNA processing and chromatin remodeling, altering global gene expression, while C6-MDB had a more discrete effect on spliceosome-mediated RNA splicing events, highlighting a more focused impact on splicing and RNA metabolism.

GO Biological Process enrichment analysis of 702 genes affected (x > 10% in PSI) by CBZ treatment (Supplemental Table 25; Additional Information) reveal critical impacts on DNA modification, RNA processing, and lipid biosynthesis, including significant involvement in DNA alkylation and methylation, pointing to CBZ’s unique influence on genetic stability and epigenetic regulation. CBZ altered splicing for genes involved in RNA processing/splicing, emphasizing a broader impact on post-transcriptional regulation of gene expression compared to the MDB compounds. Notably, RNA-binding and RNA-modification genes like RBM39, TUT4, and SRPK1 showed increased splicing only after CBZ treatment, illustrating the complex selection pressure induced by CBZ. Unlike C6-MDB, which primarily targets splicing regulation, CBZ’s effects span both genetic and metabolic landscapes, indicating a multifaceted impact on cellular physiology. Furthermore, while CBZ and C1-MDB both impact the splicing of RNA processing genes, they diverge significantly in regulating genes involved in lipid metabolism and DNA modification.

Next, a paired KEGG Pathway Analysis of all significantly spliced genes modulated by C1-MDB, C6-MDB, and CBZ was completed, confirming distinct biological impacts of each xenobiotic on both gene splicing and expression (Supplemental Tables 26–34; Additional Materials). The details of the KEGG pathway analysis are discussed in greater detail in the Supplemental Results (Additional Materials), but in general, KEGG analysis reinforced the GO biological process analysis, revealing distinct biological impacts of each xenobiotic on immune response, metabolic processes, and tissue homeostasis. The KEGG analyses highlight a shared emphasis on splicing associated with metabolic regulation, particularly for genes linked to glycerophospholipid metabolism (CHKB, DGKA, LCAT). However, unique splicing responses were also observed, such as C1-MDB’s impact on genes linked to ubiquitin-mediated proteolysis (WWP1, CBLB) and AMPK signaling genes (STK11, TSC2). C6-MDB uniquely modulated the splicing of genes linked to mitophagy and mTOR signaling (PRR5, CAB39L), while CBZ’s influenced the splicing of genes in cell growth/cancer-related (EGF, RPS6KB2, PDGFA) and ErbB signaling (EGF, RPS6KB2) pathways, reflecting a more global impact on cell proliferation and communication, than either MDB compound. These findings highlight the distinct biological processes each xenobiotic affects and the nuanced cellular adaptations mediated by alternative splicing in specific signaling pathways.

While each xenobiotic modulates a discrete set of gene expression and splicing events, we did identify 66 common genes that were comparably spliced across all 3 treatments. GO Biological Process and KEGG pathway analysis of the 66 commonly spliced genes identified mutual alternative splicing patterns among cell growth and signal transduction genes (NCOA1, EGF, RPS6KB2), including genes linked to Rap1 signaling (EGF, CTNND1, PDGFA), and the activation of stress responses (NCOA1, FAM120B; EPB41L5, EPB41, SORBS1). While each xenobiotic uniquely modulates specific gene sets, all universally alter metabolism, stress response, and cell repair processes to some extent, highlighting shared cellular adaptations to xenobiotic exposures (as highlighted in Supplemental Tables 35 and 36; Additional Information).

To gain deeper insights and visualize the relationship between xenobiotic-induced changes in gene expression and alternative splicing, we performed a 3D-clustergram analysis that paired a KEGG pathway analysis for significantly spliced genes (PSI > 10%, FDR < 5%, p < 0.05; from rMATs-turbo; Enrichr), with differential gene expression data obtained from RNAseq. Log2 Fold change data from rSeqDiff was used in this analysis, as it quantified differential gene expression levels for a higher percentage (65%) of significantly spliced genes than DESeq2 (35%; see Supplemental Table 37; Additional Information). Figure 6a presents a 3D-clustergram for C6-MDB treated PRHs, highlighting 20 of the 862 total significantly spliced genes, and their associated biological pathways. Within this 3D landscape it is notable that the majority of genes undergoing splicing modulation after xenobiotic exposure are associated with a decrease in gene expression, as highlighted by the RAB7A gene, which displayed a ~ 12% reduction in splicing (SE event) after treatment, associated with an ~ 53% decrease in total mRNA (Log2 fold change = -1.11). This was in contrast to genes like HRAS, which demonstrated significant splicing increases (RI; 11% increase) linked to an ~ 46% increase in mRNA expression (Log2 fold change = 0.55). The 3D landscape plot for the C6-MDB treatment correlates with the KEGG pathway analysis shown in Supplementary Table 30 (Additional Information), and a 2D-contour plot analysis of the same dataset is provided in Supplemental Fig. 8a (Additional Information) to facilitate data interpretation. Comparable 2D-contour and 3D-clustergram plots for C1-MDB and CBZ treatment groups are also provided in Supplemental Figs. 9 and 10 (Additional Information), and show similar trends, with a preponderance of decreased gene expression and splicing in many commonly modulated pathways.

Fig. 6

3D-Clustergram Analysis of Global Alternative Splicing and Differential Gene Expression in Xenobiotic-induced PRHs. T