Long-Read Isoform Sequencing Reveals Aroclor1260-Induced Isoform Usage in Mouse Livers

TL;DR

Exposure to Aroclor1260 changes how certain genes are expressed in mouse livers, contributing to liver disease.

Contribution

This study is the first to use long-read sequencing to identify differential transcript usage in PCB-induced liver disease.

Findings

Isoform changes in Adpgk, Blvra, Mup2, and Ndufaf6 are linked to metabolic pathways in MASLD.

PCB exposure alters transcript isoforms of genes involved in lipid metabolism and oxidative stress.

Network analysis connects these isoform changes to selenoprotein modifications observed in proteomics.

Abstract

Background/Objectives: Long-term exposure to polychlorinated biphenyls (PCBs), including the mixture of PCBs in Aroclor1260 (Ar1260), results in metabolic dysfunction-associated steatotic liver disease (MASLD) in mice and humans. While the effects of PCBs on gene expression are well-documented using short-read RNA sequencing, the regulatory roles of alternative splicing (AS) and differential transcript usage (DTU) are uncharacterized. AS has been implicated in MASLD. Previously, we reported that chronic (34 wks.) exposure of normal, low-fat-diet (LFD)-fed male mice to Ar1260 resulted in 12 hepatic RNA modifications. Proteomic analysis of these same liver samples identified Ar1260 exposure-associated changes in selenoproteins: GPX4 and SELENBP2 were increased and SELENOS and SELENOF were reduced. Methods: Here we used long-read isoform sequencing (IsoSeq) to identify DTU in four genes in the Ar1260-exposed livers: Adpgk, Blvra, Mup2, and Ndufaf6. Results: Network analysis of the corresponding proteins revealed a strong association with pathways relevant to MASLD including lipid metabolism, glycolysis, and oxidative stress. Conclusions: These findings suggest that PCB exposure alters the transcript isoform landscape of key metabolic genes involved in MASLD.

Figures

Full text

1. Introduction

Chronic caloric overload results in impaired glucose and lipid metabolism in the liver, resulting in lipid accumulation in hepatocytes (steatosis) as the starting point for the spectrum of metabolic dysfunction-associated steatotic liver disease (MASLD) [1]. The prevalence of MASLD is increasing, with 36–48% of the U.S. population affected [2]. MASLD ranges from steatosis to steatohepatitis, characterized by hepatocellular cell death, inflammation, and fibrosis with accumulation of proinflammatory monocyte-derived macrophages [3] to cirrhosis and hepatocellular cancer (HCC) [4]. There is only one therapeutic agent, Resmetirom, that is FDA-approved for treating MASLD [5], and semiglutide was recently approved for treating metabolic-associated steatohepatitis (MASH) [6]. In addition to a high-fat/Western diet (HFD), metabolic syndrome, Type II diabetes mellitus (T2DM) [7], gene variants (SNPs) [8], and gut microbiome dysbiosis [4,9] act as risk factors in MASLD development. Furthermore, exposure to metabolism-disrupting environmental chemicals, e.g., polychlorinated biphenyls (PCBs), is associated with MASLD [10]. Our research group coined the term ‘TASH’ (toxicant-associated steatohepatitis) [11] in the context of environmental pollutant exposures, including PCBs. Subsequently, TASH has been identified in the context of additional environmental pollutant exposures, including PCBs, and it has been established that environmental chemicals interact with HFD to promote steatohepatitis in a ‘two hit’ model [12].

We recently reported that chronic (34 wks.) exposure in male mice fed a normal/ low-fat diet (LFD) to Aroclor 1260 (Ar1260), a non-dioxin-like mixture of PCBs, resulted in steatohepatitis and MASLD [13]. Network analysis integrating 12 identified hepatic RNA modifications altered by Ar1260 exposure revealed specific epitranscriptomic changes linked to NRF2 (Nfe2l2). Changes in NRF2 protein abundance were validated. The results of this study demonstrated that Ar1260 exposure altered the liver epitranscriptome in pathways associated with MASLD.

Post-transcriptional alternative splicing (AS) allows a single gene to produce multiple mRNA variants, thus expanding proteomic diversity from a limited number of genes [14]. Differential transcript usage (DTU) refers to changes in the relative abundance of different transcript isoforms derived from the same gene between varying conditions, independent of changes in the overall gene expression level. Isoform switching is a specific type of DTU where the dominantly expressed isoform shifts from one experimental condition to another [15]. These changes can carry significant functional consequences, such as a shift from a protein-coding transcript to a non-coding one, or a transition between isoforms with different functional properties [16].

Environmental stressors, including various toxicants, have been shown to rapidly induce transcript-specific changes in splicing patterns [17]. Recent investigations have demonstrated that AS is dysregulated in MALSD [18]. Specifically, we reported that co-exposure of male C57Bl6/J mice fed an HFD to Ar1260 and PCB126 altered the abundance of splicing factors (SFs) and induced differential AS events (ASEs) in mouse livers, impacting genes associated with MASLD pathways [19]. This aberrant splicing can lead to the production of protein variants that drive key pathological features, including steatosis, inflammation, and fibrosis [19].

Isoform switching has functional implications in toxicology and the progression of liver disease. For example, 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD), a compound structurally and mechanistically related to dioxin-like PCBs, induces hepatic metabolic reprogramming by increasing the expression of pyruvate kinase isoform M2 (PKM2) [20]. PKM2 exhibits lower catalytic activity compared to other isoforms and redirects glycolytic flux towards pathways that support antioxidant capacity, a mechanism that contributes to the progression of steatosis to steatohepatitis and HCC [20]. A comprehensive understanding of these isoform-level changes is crucial because they represent a layer of gene regulation not captured by traditional short-read RNA-seq gene expression analysis, yet they directly impact protein function and disease pathogenesis. Given that AS is a critical post-transcriptional regulator of gene expression and protein function, and environmental stressors, including PCBs, are known to induce AS/isoform switching and alter splicing factors in the liver, a critical question arises regarding whether Ar1260-induced isoform switches directly contribute to pathology and disease progression (e.g., promoting fibrosis or carcinogenesis).

Pacific Biosciences (PacBio) Isoform Sequencing (IsoSeq) has revolutionized our understanding of isoform diversity using single-molecule real-time (SMRT) technology to generate full-length reads of mRNA [21]. This technology directly sequences complete cDNA molecules, often achieving lengths of 10–15 kilobases (kb), thus offering a significant advantage over traditional short-read sequencing technologies [22]. PacBio’s HiFi sequencing can achieve an accuracy exceeding 99.999% by generating circular consensus sequences from multiple passes around a circular DNA template [23]. Traditional short-read RNA sequencing inherently struggles to resolve full-length transcript structures, especially for long or complex transcripts necessitating computational assembly of fragmented sequences, a process prone to errors, particularly when dealing with complex isoform mixes [24]. Consequently, short-read methods often fail to accurately infer isoform-level resolution and quantify isoform switching, whereas IsoSeq’s long reads span entire transcripts, enabling resolution of isoform structures, precise identification of various alternative splicing events (e.g., intron retention, exon skipping, and alternative 5′/3′ splice sites), and accurate detection of alternative transcription start and end sites [24]. This capability extends to the discovery of novel isoforms and genes that are entirely missed by short-read approaches and the ability to identify DTU, which is not observable using short-read data, providing a more comprehensive view of gene regulation [25].

Here we used PacBio IsoSeq to sequence full-length transcripts without assembly to identify novel isoforms, characterize complex alternative splicing events, and accurately quantify DTU and isoform switching in mouse liver after chronic (34 wk.) Ar1260 exposure in mice fed a normal chow LFD. Livers were in three treatment groups: (1) vehicle control; (2) Ar1260-exposed livers with low NAFLD Activity Score (NAS) [26,27]; or (3) high NAS (>3.0) [14]. This project tested the hypothesis that Ar1260 exposure generates novel isoforms in mouse livers in pathways associated with MASLD, including apolipoproteins. The objective of the proposed research was to identify the isoform diversity in the livers of mice with MASLD after chronic Ar1260 exposure.

2. Materials and Methods

2.1. Experimental Design, Liver Tissue Collection, and RNA Isolation

The experimental design is described in [14]. The animal protocol was ratified by the University of Louisville Institutional Animal Care and Use Committee (Approval code 21871 on 15 July 2021) [28]. Briefly, eight-week-old male C57Bl/6J mice were fed ad libitum a normal LFD (20%, 69.8%, and 10.2% of total calories were from protein, carbohydrate, and fat, respectively; TekLad 06416 (Envigo, Indianapolis, IN, USA)). After two weeks, Ar1260 was administered via a single oral gavage at a dose of 20 mg/kg. The duration of exposure, due to entero-hepatic PCB recycling [29], extended to 34 weeks and resembled a chronic exposure setting. Prior to euthanasia, mice were fasted for 6 h. Liver tissue was rapidly excised from each mouse and weighed to record any changes in organ weight. Immediately following excision, tissue samples were flash-frozen in liquid nitrogen and preserved in RNAlater Stabilization Solution (ThermoFisher, Waltham, MA, USA). Total RNA was isolated using TRIzol (ThermoFisher).

2.2. PacBio IsoSeq Library Preparation and Sequencing

RNA samples were assessed for purity on the NanoVue Plus Spectrophotometer (Biochrom, Holliston, MA, USA). Quantification was performed by Qubit Fluorometer using the RNA High Sensitivity assay (ThermoFisher, Q32855). RNA integrity was assessed on the Agilent 2100 Bioanalyzer using the RNA 6000 Nano kit (Agilent, 5067-1511, Santa Clara, CA, USA) (Table S1). Barcoded cDNA was prepared from 300 ng total RNA utilizing the IsoSeq Express 2.0 Kit (PacBio, 103-071-500, Menlo Park, CA, USA) according to the manufacturer. The barcoded amplified cDNA was quantified by Qubit with the 1X dsDNA HS Assay Kit (ThermoFisher, Q33231) and quality was assessed on the Agilent Tapestation 4150 using the D5000 DNA kit (D5000 DNA reagents 5067-5589, D5000 DNA screentape 5067-5588). cDNA Barcode information is recorded in Table S2. Barcoded Kinnex Array libraries were prepared from 56 ng equimolar pooled barcoded cDNA according to “Preparing Kinnex Libraries using the Kinnex full-length RNA kit, 103-238-700, Rev03 (MAR2024)” utilizing the Kinnex PCR 8-Fold kit (PacBio, 103-071-600) and Kinnex concatenation kit (PacBio, 103-071-800). The final libraries were quantified by Qubit and quality assessed on the Tapestation using the Genomic DNA kit (Genomic DNA reagents 5067-5366, Genomic DNA screentape 5067-5365).

A total of 150 pM of the library was prepared for sequencing with the following calculations generated with SMRT Link v13.1 software (PacBio) for Kinnex full-length RNA library type and applications. Briefly, the library was annealed to the Kinnex sequencing primer and bound to polymerase prior to purification with SMRTbell cleanup beads (PacBio, 100-265-900) utilizing the Revio Polymerase Kit (102-739-100, PacBio). Sequencing was performed on a 25 M SMRT cell on a Revio platform (PacBio) and data collected through a 24 h video with adaptive loading enabled. Primary analysis of raw data to HiFi reads was performed on-instrument, after which data were exported to a local server. The raw data of the IsoSeq are available in the Gene Expression Omnibus (GEO) database: RNA accession number GSE308889.

2.3. Bioinformatics Workflow for DTU Analysis

Following sequencing, the resulting flnc.bam files were aligned to the mouse genome version GRCm39.fa using the PBMM2 aligner (v1.10.0). IsoSeq3 (v4.0.0) was used to collapse redundant transcripts into unique isoform models and measure abundance for each unique isoform. The option “do-not-collapse-extra-5exons” was used as recommended for bulk IsoSeq samples [30]. The maximum allowed 3′ difference was 100 bases on the same exon. The maximum allowed 5′ difference was 1000 bases on the same exon. Pigeon (v1.0.0) was used to classify the full-length transcripts against a reference genome to identify known and novel transcripts. Pigeon was also used to filter isoforms for primers annealing to downstream A-rich regions (IntraPriming), isoforms containing a single exon (Mono-Exonic), artifacts of reverse transcriptase (RT) template switching, and low sample coverage for splice sites that are not in the known canonical set of splice sites (LowCoverage/non-canonical). PacBio transcript IDs differ across samples for the same transcript. Therefore, a fasta file was generated for filtered isoforms using gffread (v0.12.6), converted to sam format using minimap2 (v2.24), then converted to md tagged sam files using samtools (v1.3), and added to a Talon SQLite database. Talon was used to match full-length transcripts across individual samples and identify them with a common transcript ID. The unified raw count matrix was filtered to include only transcripts with a sum of counts greater than or equal to 10 across the 12 samples. After normalization, SQANTI3 (v3.3) classification was performed on filtered annotations [31] and DTU analysis was performed using DRIMseq (v1.32.0) [32]. This method employs a Dirichlet-multinomial model to analyze transcript-level read counts, effectively accounting for both overdispersion and the correlated nature of transcript proportions within a gene. Dispersion/precision is estimated using a Maximum Likelihood (ML) approach, followed by an Empirical Bayes (EB) moderation/shrinkage. Controlling the false discovery rate (FDR). This is achieved by adjusting the raw p-values using the Benjamini–Hochberg (BH) procedure. The count matrix was filtered to include transcripts with a sum of counts greater than or equal to 10 across the 12 samples (Table S3) and normalized using trimmed mean of M-values (TMM) (Table S4). To overcome a known issue with DRIMseq when one experimental group has all zero counts for a transcript, but sufficient counts from the other experimental group and sufficient counts for the gene, a pseudocount of 0.1 was added to counts for genes with either experimental group in a comparison having zero counts. DRIMseq analysis was followed by an analysis in StageR (v1.26.0) to control overall error rate when dealing with both genes and transcripts [33].

2.4. STRING Database Protein–Protein Interaction Analysis

To gain preliminary insights into the functional network and potential protein–protein interactions of the genes exhibiting significant differential transcript usage (DTU), we performed an analysis using the Search Tool for the Retrieval of Interacting Genes/Proteins database (STRING v12.0). The genes analyzed were found to have a statistically significant DTU (q < 0.05) between the control and Ar1260-exposed groups in the Iso-Seq analysis. The analysis was restricted to the Mus musculus (Mouse) organism to ensure accurate mapping and species-specific interaction data. The STRING network was constructed using the following specific parameters to ensure focus and reliability: full string network, medium interaction confidence threshold (0.400), all prediction methods, and no more than 10 interactions per primary DTU gene. STRING automatically performed a Gene Ontology and KEGG pathway enrichment analysis on the resulting interaction network. The results were filtered using a false discovery rate (FDR) to identify enrichment within the network of DTU-affected genes.

3. Results

3.1. PacBio IsoSeq Data Quality

To investigate hepatic DTU resulting from PCB exposure, we used PacBio to sequence the livers of male C57Bl/6J mice fed an LFD with a single oral exposure the non-dioxin-like Ar1260. To account for inter-individual variability in chemical response, livers were sorted based on the severity of injury using the NAS. Livers with a histological NAS > 3.0 were classified as ‘high NAS’ as described in [13]. Four controls, four low NASs, and four high NASs were selected for Isoform sequencing. The sample size is sufficient because the focus is on the magnitude of the observed biological effect and the statistical methods used, which are designed to compensate for lower N. Following sequencing, 12 files holding full-length non-chimeric (FLNC) reads were obtained, 4 replicates each for each of the experimental groups. Sample descriptions can be found in Table 1. The bioinformatics workflow for DTU analysis is a multi-step process designed to extract comprehensive and accurate information from full-length transcript reads (Figure 1A). The average number of reads per sample was ~5.5 million and the average number of transcripts was ~500,000 (Table 2).

SQANTI3 systematically classifies gene isoforms into known and novel categories. Known isoforms are further delineated as full splice matches (FSMs), which completely align with existing splice site annotations, or incomplete splice matches (ISMs), representing partial splice site matches or alternative 3′ and 5′ ends (Figure 1B). Beyond known forms, SQANTI3 identifies novel isoforms previously unannotated in the genome. These novel classifications encompass transcripts with novel splice junctions, intron retention, mono-exon retention, and novel splice sites, providing a comprehensive characterization of the isoform landscape (Figure 1B). The analysis of the transcriptome landscape identified a substantial number of unique isoforms: 12,501 annotated genes and 433 novel genes (Figure 1C). Further classification identified splice junctions: known canonical (95,471), known non-canonical (35), novel canonical (9787), and novel non-canonical (5) (Figure 1D).

SQANTI3 performs quality control and curation of IsoSeq data, which often contains technical artifacts that can complicate downstream analysis. It classifies transcripts based on a range of parameters to filter out low-quality isoforms and ensure a reliable transcriptome. Here, SQANTI3 addressed reference transcript redundancy, a common issue in gene annotation that occurs when multiple reference transcripts in a database represent the exact same gene product or splicing variant by creating a non-redundant set of known isoforms. By filtering out technical artifacts from long-read sequencing data, SQANTI3 provides a more accurate baseline for comparison (Figure S1A–C). Additionally, SQANTI3 identified and removed intra-priming artifacts, which occur when the oligo-dT primer used during reverse transcription binds to an internal adenine-rich region of a transcript instead of the genuine poly-A tail, resulting in truncated transcripts [32] (Figure S1D–G). SQANTI flagged and filtered out RT switching that occurred when the RT enzyme switched templates during cDNA synthesis (Figure S1E), non-canonical splice junctions (Figure S1F), and nonsense-mediated decay (Figure S1G). By applying these specific filters, SQANTI3 cleaned the dataset, providing a more accurate and high-quality set of isoforms for further analysis [31].

3.2. Transcriptome Classification and Isoform Novelty

The analysis of the PacBio IsoSeq data using SQANTI3 revealed a complex and diverse transcriptome, with isoforms classified into several structural and genomic categories. The majority of genes expressed a low number of isoforms (Figure 2A), although a subset of genes demonstrated a high degree of AS. When distinguishing between known and novel transcripts, we found that a substantial number of isoforms per gene were classified as novel (Figure 2B), indicating significant new discoveries within the transcriptome. Further analysis of transcript length across these structural categories revealed that novel isoforms for known genes (NNC and NIC) generally had a different length distribution compared to the annotated (FSM and ISM) transcripts (Figure 2C).

The overall isoform population was categorized into four major structural groups (Figure 2D). A large proportion of the isoforms were classified as full splice matches (FSMs), i.e., they aligned perfectly with a known reference transcript. A smaller but significant fraction were categorized as incomplete splice matches (ISMs), representing truncated versions of annotated isoforms, possibly due to technical artifacts or genuine AS events. The remaining isoforms were classified as novel, with the novel not in catalog (NNC) category comprising transcripts with one or more new splice junctions, and the novel in catalog (NIC) category representing new combinations of existing splice sites. Detailed breakdowns of the structural categories provided a more granular view of the transcriptome. The FSM category (Figure S2A) and ISM category (Figure S2B) showed a diversity of known transcripts. The NNC category (Figure S2C) highlighted transcripts with previously uncharacterized splice sites, while the NIC category (Figure S2D) pointed to novel AS patterns.

In terms of genomic location, most isoforms were found within the boundaries of known genes (genic genomic, Figure S2E). However, a small percentage of transcripts were mapped to other genomic regions. These included isoforms transcribed from the opposite strand of a gene (antisense, Figure S2F), transcripts spanning across two or more adjacent genes (fusion, Figure S2G), and a minor but important proportion of isoforms originating from intergenic regions (intergenic, Figure S2H), suggesting the presence of novel genes or long noncoding RNAs (lncRNAs) not yet annotated in the reference genome. These results provide a comprehensive catalog of all detected isoforms, including the annotation and structural classification of novel isoforms across all samples. The presence of these previously unannotated transcripts provides foundational evidence that Ar1260 exposure significantly affects the transcriptional output by potentially activating new splicing events. While this initial analysis does not distinguish between experimental groups, the identification of these novel isoforms serves as a strong indication of transcriptional plasticity with potentially significant functional implications in disease progression.

3.3. Ar1260 Exposure Results in DTU

We performed three-dimensional gene-level principal component analysis between three sample groups, corn oil (control), Ar1260 with low NAS, and Ar1260 with high NAS. The primary source of variance in the data was the biological differences between the control and Ar1260 high NAS (Figure 3A). To identify DTU, the DRIMseq package was utilized. The analysis was conducted in a two-stage hierarchical testing procedure to ascertain statistical significance. An initial likelihood ratio test was performed for each gene in two pairwise comparisons (Ar1260 versus control, and Ar1260 high NAS versus control) to evaluate the null hypothesis of no DTU between conditions (Tables S5 and S6). For all comparisons, genes with a p < 0.05 proceeded to the next stage of analysis. A second likelihood ratio test was performed at the individual transcript level to identify which specific isoforms exhibited significant proportional changes (Tables S7 and S8). To maintain control over the FDR across both stages of this analysis, the StageR method was applied. StageR integrates the results of the two-stage testing procedure to provide a single, adjusted p-value for each gene and transcript, ensuring robust FDR control. A significance threshold of an adjusted p-value < 0.05, calculated using the Benjamini–Hochberg procedure to control the FDR, was applied to determine genes and transcripts with significant DTU.

The DRIMseq analysis, followed by StageR, identified several genes with significant DTU across comparisons. A summary of these findings is provided in Table 3 and Table 4, which detail the number of significant genes and transcripts for each pairwise comparison. In the Ar1260 (low + high NAS, combined) versus control comparison, we found significant DTU in six genes, affecting a total of 20 transcripts. These genes included Adpgk (two transcripts), Blvra (two transcripts, one of which was a novel isoform), Cops7a (seven transcripts, five of which were novel), Hsd17b6 (three transcripts, two novel), Mup2 (four transcripts, three novel), and Ndufaf6 (two transcripts, one novel) (Figure 3B). In Ar1260 high NAS versus control, DTU was identified in nine genes, affecting a total of 28 transcripts. These genes included Clk1 (five transcripts, four novel), Cops7a (three transcripts, two novel), Eif4a2 (four transcripts, three novel), Ewsr1 (two transcripts, one novel, Mup2 (four transcripts, three novel), Mup3 (two transcripts), Nat9 (three transcripts, one novel), Tef (two transcripts, one novel), and Vapa (three transcripts, two novel) (Figure 3C). Custom DTU tracks were generated from BED files where read counts were normalized by total gene expression to reflect isoform proportions. (Figures S3–S15).

3.4. Isoform Switching with Ar1260 Exposure Is Independent of Gene Expression or Global Alternative Splicing Changes

We previously performed short-read Illumina RNA-seq on the livers of the chronic Ar1260-exposed mice and analyzed the differential gene expression (DEG) only in relation to global RNA modification changes with Ar1260 exposure [13]. Here we further analyzed our RNA-seq data deposited in the GEO database (GSE17327) [13] to examine DEG in all livers exposed to Ar1260 (combined low + high NAS) (Figure 4A) and Ar1260-exposed livers with high NASs (Figure 4B) compared with control. To investigate the full transcriptional impact of Ar1260, we integrated RNA-Seq and IsoSeq data. This combined approach allowed for a comprehensive analysis of both total gene expression and isoform-specific regulation. Our analysis revealed a notable finding: a subset of genes demonstrated DTU in response to Ar1260 exposure, yet these same genes showed no significant changes in overall gene expression. This suggests that Ar1260 can modulate gene function by altering transcript splicing without affecting the total number of transcripts produced.

To evaluate the effect of Ar1260 on AS, we performed differential splicing analysis using rMATS (replicate multivariate analysis of transcript splicing, v4.1.2) on the short-read RNA-seq data [34]. Our short-read dataset contained an average of 3 to 5 million mapped reads per sample, sufficient for the effective detection and quantification of splice junctions and individual splicing events by rMATS. The rMATS tool identifies and quantifies different types of alternative splicing events from RNA-seq data, such as exon skipping, intron retention, and alternative splice sites, providing a statistical framework to determine significant changes between sample groups [34]. Ar1260 treatment minimally affected splicing in genes with DTU, with only two genes, Mt1 and Naa10, exhibiting significant changes in retained intron events (Tables S9 and S10).

Further investigation into our RNA-seq data revealed that Ar1260 treatment had a limited impact on the expression of known alternative splicing factors. In Ar1260 versus the control, only three alternative splicing factor genes were significantly upregulated: Srsf5 (serine and arginine rich splicing factor 5), Htatsf1 (HIV-1 Tat specific factor 1), and Polr2a (RNA polymerase II subunit A), which plays a role in transcription and splicing (Figure 4C). In the high NAS subgroup, two alternative splicing factor genes were upregulated (Srsf5 and Polr2a), while four were downregulated: Snrpb (small nuclear ribonucleoprotein polypeptides B and B1), Yap1 (Yes-associated protein 1), Puf60 (Poly binding-splicing factor PUF60), and Zfp207 (zinc finger protein 207) (Figure 4D). These results suggest that, while Ar1260 has a minimal effect on the overall landscape of hepatic AS events, the expression of a small, specific set of splicing factors is significantly affected by Ar1260, particularly in the high NAS subgroup. This selective modulation may contribute to the subtle isoform differences observed.

3.5. Correlation of DTU with Liver Pathology and Metabolic Alterations

To understand the potential impact of Ar1260-induced DTU on liver pathology, we utilized STRINGdb (v12.0) [35] to identify correlations between the altered genes and known hepatic injury pathways. Our analysis focused on investigating how changes in gene isoforms might affect protein–protein interactions.

This analysis revealed that, compared to the control, four of the genes with DTU in mice exposed to Ar1260, Adpgk, Blvra, Mup2, and Ndufaf6, showed significant enrichment in several KEGG pathways [36] associated with hepatic injury. Specifically, interactions involving Adpgk were enriched in pathways related to glycolysis and gluconeogenesis, as well as the insulin signaling pathway and Type II diabetes (Figure 5A). The altered isoforms of Blvra were associated with pathways for drug metabolism (cytochrome P450), bile secretion, and other general metabolic pathways (Figure 5B). The protein interactions of Mup2 were linked to the metabolism of xenobiotics and the inflammatory mediator regulation of Transient Receptor Potential (TRP) channels (Figure 5C). Finally, Ndufaf6 and its interactions with Ndufa10 were significantly enriched in pathways associated with MASLD (Figure 5D). While the other two genes identified to exhibit DTU in Ar1260 compared to the control group (Cops7a and Hsd17b6) were not enriched in pathways directly involved in MASLD, they were enriched in pathways indirectly associated with metabolism-associated disorders (Figure S16). Genes with DTU in the Ar1260-exposed high NAS versus control subgroup demonstrated interactions with Ugt2b1 enriched in pathways associated with liver dysfunction, including metabolism of xenobiotics, inflammatory mediator regulation of TRP channels, and metabolic pathways (Figure S17).

We analyzed the amino acid sequences encoded by identified transcripts to determine how DTU would alter protein products. In Ar1260-exposed mouse livers, the Adpgk gene primarily utilized the ENSMUST00000026266.9 transcript, which resulted in a protein isoform lacking a glutamine residue at position 313 (Figure S18). This may affect protein structure. For the Blvra gene, the canonical transcript (ENSMUST00000002064.15) was more abundant with Ar1260 exposure. However, IsoSeq analysis identified a novel transcript in control samples featuring a truncated 5′ end of exon 1 (Figure S19). The result of this truncation was an unmodeled region of the protein (Figure S19C); thus, the potential impact is unknown. For Mup2, the established transcript ENSMUST00000074700.9 was dominant in unexposed mice, but three novel, truncated transcripts (TALNOT000300526, TALNOT000300560, and TALNOT000300598) were utilized in the Ar1260-exposed group. These novel transcripts possess fewer exons, which may affect chemical and insulin receptor interaction (Figure S20), although no one has examined MUP2 direct interactions in mouse livers. Finally, Ndufaf6 showed a shift from the known, control-dominant transcript (ENSMUST00000058183.9) to a novel transcript (TALONT000158320) in the Ar1260-exposed group, characterized by the omission of exon 6 (Figure S21). While the amino acid sequence structure is unavailable, the skipped exon is expected to alter protein structure.

These results demonstrate that Ar1260-induced changes in gene isoforms correlate with key metabolic and liver disease pathways, suggesting that Ar1260 effects on DTU may contribute to hepatic pathology by disrupting protein–protein interactions linked to glycolysis, drug metabolism, and MASLD.

4. Discussion

In the development and progression of human metabolic diseases, a complex interplay of genetic predisposition, environmental factors, and lifestyle choices leads to disruptions in key metabolic pathways [37,38]. The molecular mechanisms by which PCB exposure affects the liver, leading to MASLD and its progression to MASH, remain to be resolved. We report that a single oral exposure of male C57Bl/6J mice to the PCB mixture Ar1260 alters liver pathology (MASLD and MASH) and the hepatic proteome [28], transcriptome [39], and epitranscriptome [13], depending on the diet of the mice and time of overall PCB exposure.

AS is implicated in MASLD [40]. Here, we examined the impact of chronic Ar1260 exposure in mice fed a normal (LFD) diet on AS and isoform switching in the liver. The observation that DTU is independent of DGE underscores a specific and important regulatory mechanism. While DGE measures changes in the total abundance of transcripts from a given gene, DTU quantifies changes in the relative proportions of gene-specific different transcript isoforms. In isoform switching, a gene continues to use the same canonical set of exons/introns, but the ratio between the resulting full-length transcripts changes [41]. Here, we observed a change in the relative usage of existing transcripts (an isoform switch) where DTU is significant, but the underlying set of splice sites used did not cross the statistical threshold for significance. Our findings indicate that Ar1260 exposure promotes the production of a different set of isoforms from the same gene locus, independent of global alternative splicing patterns, potentially leading to the synthesis of distinct protein products or altering the regulatory capacity of noncoding transcripts. This effect on isoform balance is overlooked by standard, short-read RNA-Seq analyses, which focus exclusively on total gene expression. Therefore, the integration of IsoSeq data was critical to uncover this nuanced level of transcriptional regulation induced by Ar1260.

Ar1260 exposure resulted in significant DTU in four genes: Adpgk, Blvra, Mup2, and Ndufaf6. Exposure to Ar1260 resulted in the predominant expression of novel or distinct transcript isoforms—such as an exon-truncated Mup2 and an exon-skipped Ndufaf6—which are predicted to fundamentally alter the amino acid sequence and potentially compromise protein function, particularly domain integrity. The protein interaction network of the four proteins encoded by these genes suggests a multi-faceted dysregulation of key metabolic and cellular processes that may contribute to MASLD (Figure 5). The altered expression of Ndufaf6 transcripts may play a role in previously observed Ar1260-induced hepatic mitochondrial dysfunction [42]. NDUFAF6 is a critical assembly factor for mitochondrial Complex I [43]. Impaired Complex I function is a hallmark of mitochondrial stress, resulting in reduced ATP production, increased reactive oxygen species (ROS) generation, and the accumulation of unoxidized fatty acids in MASLD [44,45]. The Ar1260-mediated perturbation of Ndufaf6 transcript usage provides a mechanistic link between the Ar1260 MASLD phenotype.

DTU of the Adpgk and Blvra transcripts may play a role in liver metabolic disruption in response to Ar1260. ADPGK is a key enzyme in glycolysis, catalyzing the phosphorylation of glucose to glucose-6-P [46]. Altered APDGK perturbs hepatic glucose homeostasis, resulting in increased de novo lipogenesis and impaired insulin responsiveness, thereby exacerbating the hepatic and extra-hepatic metabolic syndrome features in MASLD [47]. Proteins interacting with BLVRA (biliverdin reductase A) were enriched in several key pathways associated with hepatic injury, including xenobiotic metabolism by cytochrome P450s, bile secretion, and steroid hormone biosynthesis. Deregulation of steroid hormone biosynthesis can impact lipid and glucose metabolism, key drivers of MASLD [48]. For example, an altered estrogen-to-androgen ratio leads increased fat accumulation in the liver [48]. Steroid hormone biosynthesis is linked to cholesterol metabolism that contributes to dyslipidemia in MASLD [49]. BLVRA is a stress-response protein and the product of its enzymatic activity, bilirubin, is an antioxidant [50]. BLVRA is a cell surface receptor that activates PI3K, acts as a nuclear transcription factor, has important roles in mitochondrial metabolism and in inflammation, and is a potential therapeutic target (reviewed in [51,52]). The dysregulation of Blvra isoform usage in response to Ar1260 exposure may compromise the liver’s ability to manage oxidative stress, creating a cellular environment conducive to MASLD progression.

The Ar1260-mediated DTU of hepatic Mup2 also points to a disruption in lipid and energy metabolism. Proteins interacting with MUP2 are associated with inflammatory mediator regulation of TRP channels that mediate immune and inflammatory responses [53]. TRP channels are expressed in hepatocytes and other liver cell types, where they play a role in cellular processes and the pathogenesis of HCC [54]. The observation that Ar1260-mediated DTU of Mup2 links to disrupted lipid/energy metabolism and inflammatory signaling via TRP channels provides a clear molecular mechanism for the severe tissue damage seen in the ‘high NAS’ phenotype. Major urinary proteins (MUPs) secreted by the liver act as transporters of small lipophilic molecules linked to regulating systemic energy expenditure and lipid homeostasis [55]. Given the central role of lipid accumulation in MASLD [56], an alteration in the expression or function of Mup2 due to DTUs could directly impact lipid transport and storage within the liver, contributing to the development of hepatic steatosis.

Limitations

Although our IsoSeq analysis of Ar1260-exposed mouse livers provided novel insights into the DTU of key hepatic genes with known functional roles in MASLD, there are limitations to this analysis. While IsoSeq offers high-quality, full-length transcript data, it may not capture all low-abundance isoforms or subtle splicing events across the entire transcriptome. The sequencing depth, while robust, represents a single snapshot in time and may not fully reflect the dynamic, temporal changes in transcript usage mediated by the circadian clock [57,58] that occur during the progression of MASLD. Our starting material was polyA-selected RNA from livers, which represents a bulk sample from the entire organ. Consequently, any observed changes in gene expression or RNA processing represent a cumulative, average effect across the whole liver, not cell-specific changes. Single-cell or single-nucleus RNA-seq will be informative to identify liver and immune cells responsive to Ar1260. Our analysis focused on DTU in four specific genes. While their roles in metabolism and mitochondrial function are well-established, it is likely that additional genes and their isoforms also contribute to the observed MASLD phenotype. Future research should aim to integrate IsoSeq data with proteomics and metabolomics to provide a more comprehensive view of the liver response to PCB exposure. Conclusively determining how the observed isoform switches directly translate to changes in the final protein structure or function for the four genes highlighted will require in-depth domain analysis (loss/gain), prediction of protein–protein interaction site alteration, and integrating the available proteomics data from similar treatments. Validating the functional significance of the identified isoforms and exploring the cell-specific effects of these changes will be critical steps toward a more complete understanding of how these legacy contaminants drive liver disease. Finally, examining the expression of these four DTUs in human MASLD and MASH samples is required for further translational relevance.

5. Conclusions

The combined protein interaction and KEGG pathway analysis of Adpgk, Blvra, Mup2, and Ndufaf6 provides a systems-level analysis of how Ar1260 exposure disrupts core hepatic metabolic processes in this mouse model of Ar1260-exposure-induced MASLD. This analysis was supported by our IsoSeq data and demonstrated that DTU in these genes likely results in altered protein isoforms. These isoforms are predicted to have distinct structural and functional properties, thereby re-wiring protein–protein interactions and causing changes in cellular metabolism and homeostasis. The functional associations of NDUFAF6 with oxidative phosphorylation, ADPGK with glycolysis, and both BLVRA and MUP2 with oxidative stress and lipid metabolism, respectively, indicate that Ar1260 exposure disrupts a complex network of interacting processes. The direct molecular link between the IsoSeq-identified isoforms and the network-level changes provides a compelling hypothesis for the role of Ar1260 in MASLD pathogenesis. Future studies will focus on validating the existence of these specific protein isoforms and assessing their impact on protein–protein binding affinities and enzymatic activity.