Abstract Background Hub1, a conserved ubiquitin-like protein, is essential for pre-mRNA splicing and transcriptional regulation in Saccharomyces cerevisiae. Despite its known functions, the genome-wide effects of Hub1 overexpression remain largely uncharacterized. This study investigates the transcriptomic and splicing landscape changes triggered by Hub1 overexpression using an integrative bioinformatic approach. Results We analyzed RNA-seq data from the [32]GSE84215 dataset, employing differential expression, alternative splicing, functional enrichment, and network-based methods. DESeq2 identified 3,915 differentially expressed genes (DEGs; 1,964 upregulated, 1,951 downregulated, padj < 0.05), demonstrating extensive transcriptional reprogramming. Principal component analysis revealed that Hub1 overexpression explained 98% of transcriptional variance, indicating its dominant regulatory influence. Using rMATS, we detected seven exon skipping events, with DYN2 showing significant differential splicing (FDR = 0.0481, ΔPSI = − 0.036). MaxEntScan analysis confirmed that DYN2’s 5′ splice site is significantly weaker than canonical yeast splice sites (score = − 18.32, p = 0.03), consistent with Hub1’s role in facilitating non-consensus splicing. Functional enrichment analyses revealed metabolic reprogramming, with upregulated pathways including biosynthesis of secondary metabolites and carbon metabolism, while growth-related processes like ribosome biogenesis and cell cycle were downregulated. Gene Set Enrichment Analysis (GSEA) further supported stress response activation (p53 signaling, NES = 1.255) and cell cycle suppression (NES = − 0.692). Weighted Gene Co-expression Network Analysis (WGCNA) identified 61 co-expression modules, with the brown module highly correlated with Hub1 overexpression (r = 0.99, p < 0.001) and enriched in biosynthetic and proteasome pathways. Protein-protein interaction network analysis revealed 35 Hub1 interactors, including spliceosomal components, reinforcing its central role in RNA processing. Conclusion Our findings reveal that Hub1 overexpression drives coordinated transcriptional and post-transcriptional changes, promoting metabolic reprogramming while specifically modulating splicing of genes with weak splice sites like DYN2. These results establish Hub1 as a dual regulator linking transcriptional control with splicing precision, suggesting a regulatory mechanism that enhances cellular adaptability under stress conditions. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-025-12006-w. Keywords: Hub1, Saccharomyces cerevisiae, RNA-seq, Alternative splicing, DYN2, Transcriptomics Introduction Post-transcriptional regulation plays a pivotal role in shaping cellular responses to environmental and physiological stimuli [[33]1]. A central mechanism within this regulatory layer is pre-mRNA splicing, which enables the generation of multiple transcript isoforms from a single gene, thereby expanding proteomic diversity and functional complexity [[34]2]. This process is tightly regulated by various proteins, including a distinct class known as ubiquitin-like proteins (UBLs) [[35]3]. Unlike classical ubiquitin, which typically mediates protein degradation via covalent modification, UBLs such as Hub1 (UBL5 in higher eukaryotes) exert their influence through non-covalent interactions, especially with components of the RNA processing machinery [[36]4, [37]5]. Hub1 is an evolutionarily conserved UBL that modulates splicing fidelity and transcript maturation across eukaryotes [[38]4]. In Saccharomyces cerevisiae, Hub1 enhances the splicing of pre-mRNAs with weak or non-consensus splice sites by interacting with spliceosomal proteins like Snu66 [[39]6]. Although S. cerevisiae contains relatively few intron-bearing genes (~ 5%), many of these genes feature non-canonical splice sites, making them highly reliant on accessory factors like Hub1 [[40]7]. Through this mechanism, Hub1 increases the efficiency and flexibility of the splicing machinery, facilitating the inclusion of otherwise poorly recognized exons [[41]8]. Beyond its role in splicing, Hub1 has been implicated in diverse biological processes, including mitochondrial maintenance, transcriptional control, and cellular stress adaptation [[42]9]. Studies have shown that Hub1 activity is modulated under stress conditions such as DNA damage and oxidative stress, suggesting that it may serve as a molecular bridge linking environmental sensing to RNA processing efficiency. This regulatory versatility makes Hub1 an important factor in maintaining cellular homeostasis [[43]5]. Despite extensive characterization of its splicing-related functions, the transcriptome-wide effects of Hub1 overexpression remain incompletely understood. While previous studies have identified individual Hub1-regulated genes, a systems-level understanding of how elevated Hub1 levels reshape the expression and splicing landscape is still lacking [[44]10]. Given the conservation of Hub1 across species and its emerging relevance in human disease contexts, such as splicing disorders and stress-related pathologies, elucidating its genome-wide regulatory influence in a model organism is both timely and important [[45]4]. The budding yeast S. cerevisiae provides an ideal system for investigating these mechanisms. Its compact and well-annotated genome, coupled with low intron density and genetic tractability, allows for precise dissection of RNA processing events with minimal background complexity [[46]11]. Additionally, the conservation of fundamental splicing components between yeast and humans ensures that insights gained from yeast studies can inform broader biological understanding [[47]12]. In this study, we perform a comprehensive integrative analysis of publicly available RNA-seq data from S. cerevisiae strains overexpressing Hub1. Our pipeline incorporates multiple layers of transcriptomic analysis, including differential gene expression profiling, detection of alternative splicing events, splice site strength evaluation, functional enrichment, weighted gene co-expression network analysis (WGCNA), and protein-protein interaction (PPI) network construction [[48]13–[49]15]. This multifaceted approach enables a detailed assessment of both transcriptional and post-transcriptional changes induced by Hub1 overexpression. The analysis reveals that Hub1 drives widespread transcriptional reprogramming, significantly affecting genes involved in metabolic processes, ribosome biogenesis, and stress response [[50]16, [51]17]. Simultaneously, Hub1 selectively modulates splicing events in a subset of genes characterized by weak splice sites such as DYN2, a gene involved in cytoskeletal organization and membrane dynamics [[52]18]. Moreover, co-expression and PPI network analyses position Hub1 as a central regulatory node with connections to key pathways involved in RNA processing and cellular adaptation [[53]19, [54]20]. Overall, this study positions Hub1 as a dual-function regulator that integrates transcriptional control with splicing precision. By coordinating these two regulatory layers, Hub1 enhances the cell’s ability to fine-tune gene expression programs in response to internal and external cues, underscoring its role in maintaining cellular homeostasis and adaptive capacity. Methods Data acquisition and quality assessment RNA-seq data for Saccharomyces cerevisiae overexpressing Hub1 and wild-type controls were retrieved from the NCBI Gene Expression Omnibus (GEO) under accession number [55]GSE84215 (accessed on July 5, 2025) [[56]21]. The dataset comprised six paired-end RNA-seq samples: three biological replicates for the xrn1 rat1-1 control condition (SRR3822520, SRR3822521, SRR3822522) and three for the Hub1-overexpression condition (SRR3822526, SRR3822527, SRR3822528). Initial quality assessment was performed using FastQC v0.11.9 to evaluate read quality, adapter contamination, and base composition [[57]22]. All samples passed quality thresholds with Phred scores > 30 for > 95% of bases and no significant adapter contamination, eliminating the need for additional preprocessing. Read alignment and processing High-quality reads were aligned to the Saccharomyces cerevisiae reference genome (Ensembl R64-1-1, release 114) using HISAT2 v2.2.1with default parameters optimized for eukaryotic transcriptomes [[58]23]. The HISAT2 index was constructed from the corresponding genome FASTA and GTF annotation files (Saccharomyces_cerevisiae.R64-1-1.114). Alignments were performed using the --dta (downstream transcriptome analysis) option to preserve alignment information necessary for transcript assembly and splicing analysis. The resulting SAM files were converted to BAM format, coordinate-sorted, and indexed using SAMtools v1.10 [[59]24]. Alignment statistics were assessed to ensure > 85% mapping rates across all samples. Gene expression quantification and normalization Gene-level read quantification was conducted using featureCounts from the Subread package v2.0.3 [[60]25]. Parameters were set for paired-end mode (-p), unique mapping only (-Q 10), and strand-specific counting when applicable. The Ensembl GTF annotation file served as the gene model reference for read assignment. Resulting raw count matrices were filtered to remove genes with fewer than 10 total reads across all samples to reduce noise from lowly expressed transcripts. For transcript-level analysis, Salmon v1.10.1 was employed using the Ensembl R64-1-1 transcriptome FASTA to build the index [[61]26]. Quantification was performed in paired-end quasi-mapping mode with the --validateMappings flag to improve alignment accuracy and --gcBias to correct for GC content bias. Output files containing transcript-level counts and TPM (Transcripts Per Million) values were imported into R using tximport v1.18.0 for downstream analysis [[62]27]. Differential gene expression analysis Differential expression analysis was performed using DESeq2 v1.30.1 [[63]13] in R v4.1.1. Raw count matrices were normalized using the median-of-ratios method to account for library size differences and compositional bias. Dispersion parameters were estimated using empirical Bayes shrinkage with the fitType="parametric” option. A negative binomial generalized linear model was fitted to model count data, and statistical significance was assessed using the Wald test. Multiple testing correction was applied using the Benjamini-Hochberg method, with genes considered significantly differentially expressed at adjusted p-value < 0.05 [[64]28]. Log2 fold change shrinkage was applied using the apeglm method to reduce variance for lowly expressed genes. Principal component analysis (PCA) was performed on variance-stabilized transformed (VST) data using the top 500 most variable genes to assess sample clustering and batch effects. Visualization included volcano plots (log2 fold change vs. -log10 adjusted p-value) and MA plots (mean expression vs. log2 fold change) to evaluate the distribution and magnitude of expression changes [[65]29]. Alternative splicing event detection Alternative splicing analysis was conducted using rMATS-turbo v4.1.2 [[66]14], specifically designed for detecting differential splicing events from replicated RNA-seq data. Analysis parameters included paired-end reads (read length = 150 bp), minimum read coverage of 10 for inclusion/exclusion counts, and false discovery rate (FDR) threshold of 0.05. The analysis focused on five canonical splicing event types: skipped exons (SE), alternative 5’ splice sites (A5SS), alternative 3’ splice sites (A3SS), mutually exclusive exons (MXE), and retained introns (RI). Percent spliced-in (Ψ) values were calculated for each condition, and differential Ψ values (ΔΨ) quantified the magnitude of splicing regulation. Significant events were identified based on FDR < 0.05 and |ΔΨ| > 0.03. Sashimi plots were generated using rmats2sashimiplot to visualize junction read coverage and validate splicing events [[67]30]. Complementary exon usage analysis To provide an alternative perspective on splicing regulation, differential exon usage analysis was performed using DEXSeq v1.48.0 [[68]31]. This method models exon usage independently of predefined splicing events and can detect subtle changes in exon inclusion. The analysis used the same BAM files and GTF annotation, with significance determined at adjusted p-value < 0.05 after Benjamini-Hochberg correction. Splice site strength assessment Splice site motif strength was evaluated using MaxEntScan to assess whether Hub1-regulated exons contained suboptimal splice sites [[69]32]. Genomic coordinates for alternatively spliced exons were extracted, and flanking sequences were retrieved using BEDTools v2.30.0 [[70]33]. Specifically, 9 bp sequences surrounding 5’ donor sites and 23 bp sequences around 3’ acceptor sites were extracted from the reference genome. These sequences were scored using MaxEntScan’s maximum entropy models (score5.pl and score3.pl), with lower scores indicating weaker splice site strength. Statistical comparison against canonical splice site distributions was performed using two-tailed t-tests [[71]34]. Functional enrichment analysis Functional annotation of differentially expressed genes was performed using clusterProfiler v4.8.3 in R [[72]35]. Gene identifiers were converted from Ensembl to Entrez format using the org.Sc.sgd.db annotation package. Over-representation analysis (ORA) was conducted for Gene Ontology (GO) terms across three domains: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF), using the enrichGO() function. KEGG pathway enrichment was performed using enrichKEGG() with organism code ‘sce’ for Saccharomyces cerevisiae. Significance was determined at FDR-adjusted p-value < 0.05. Gene Set Enrichment Analysis (GSEA) was conducted using gseGO() and gseKEGG() functions on a pre-ranked gene list sorted by log2 fold change values from DESeq2 [[73]36]. GSEA parameters included 1000 permutations and significance threshold of FDR < 0.2 to capture subtle but coordinated pathway-level changes [[74]37]. Weighted gene Co-expression network analysis Co-expression network analysis was performed using WGCNA v1.71 to identify modules of co-regulated genes associated with Hub1 overexpression [[75]15]. Low-variance genes below the 25th percentile was filtered out, and variance-stabilized counts from DESeq2 were used as input. A signed network was constructed with soft-thresholding power set to 12, determined by scale-free topology model fit criteria (R² > 0.8). Network construction and module detection were performed using the blockwiseModules() function with parameters: minimum module size = 30, deepSplit = 2, and mergeCutHeight = 0.25. Module eigengenes were calculated and correlated with the Hub1-overexpression trait using Pearson correlation. Modules with |correlation| > 0.8 and p-value < 0.05 were considered significantly associated with the experimental condition. Protein-Protein interaction network analysis Protein-protein interaction networks were constructed using STRINGdb v2.10.0 to examine interaction patterns among Hub1 and differentially expressed genes [[76]19]. Hub1 and its first-degree interactors were retrieved with a minimum combined confidence score of 400, encompassing experimental, database, and computational evidence. The resulting network was visualized using ggraph v2.0.5, with nodes colored by log2 fold change values to highlight expression patterns within the interaction network. Results Differential gene expression analysis To explore the global transcriptional response to Hub1 overexpression, RNA-seq data from three biological replicates per condition (Hub1-overexpressing vs. control) were analyzed using DESeq2. After normalization and statistical testing, a total of 3,915 genes were found to be differentially expressed with an adjusted p-value < 0.05. Of these, 1,964 genes were significantly upregulated, and 1,951 were downregulated, indicating a balanced but extensive transcriptional shift in response to Hub1. The volcano plot (Fig. [77]1A) visualizes the distribution of all genes based on statistical significance (–log10 adjusted p-value) and magnitude of expression change (log2 fold change). A substantial number of genes clustered at both ends of the x-axis demonstrate the strong upregulation (red dots) and downregulation (blue dots) of transcripts in the Hub1-overexpressing condition. Several of the most strongly regulated genes, such as YJR005C-A (log2FC = + 4.97), YOR382W (+ 4.11), and YOR383C (+ 3.55), appear in the upper-right quadrant of the plot, highlighting their elevated expression and statistical significance. On the opposite end, genes such as YHR124W and YGR250C displayed notable downregulation. To further examine gene expression trends relative to average abundance, an MA plot was constructed (Fig. [78]1B). This plot shows the relationship between mean expression (log counts per million) and log2 fold change for each gene. Significantly upregulated genes (in red) and downregulated genes (in blue) are symmetrically distributed around the zero-fold change line, confirming a global and balanced regulatory effect. Genes with high expression levels and extreme fold changes, such as YEL065W, YEL064C, and YMR058W, stand out as potential targets for further functional investigation. Together, these results suggest that Hub1 overexpression induces a strong and widespread effect on yeast gene expression, impacting both low- and high-abundance transcripts and implicating genes involved in RNA metabolism, stress response, and transcriptional regulation. Fig. 1. [79]Fig. 1 [80]Open in a new tab Differential gene expression analysis A Volcano plot showing differentially expressed genes between Hub1-overexpressing and control yeast strains. B MA plot of mean expression versus log2 fold change in Hub1-overexpressing yeast To highlight the most significant transcriptional changes, we listed the top 10 differentially expressed genes ranked by adjusted p-value and log2 fold change (Table [81]1). These genes showed substantial upregulation in Hub1-overexpressing samples compared to control, with many involved in RNA processing, stress response, and cell cycle regulation. The inclusion of genes such as YOR382W, YJR005C-A, and YMR059W underscores the strong transcriptional reprogramming driven by Hub1 overexpression. To further examine the expression dynamics of the most significantly regulated genes, a heatmap was constructed using the top 30 differentially expressed genes (DEGs), ranked by adjusted p-value and absolute log₂ fold change. Normalized expression values were transformed to Z-scores to highlight differential trends across samples. As shown in Fig. [82]2, hierarchical clustering revealed clear and consistent separation between Hub1-overexpressing and control groups. All three biological replicates per condition clustered tightly, confirming sample reproducibility. Genes upregulated in Hub1-OE samples exhibited coordinated increased expression, while downregulated genes were consistently suppressed in control samples. Notably, several genes from the top 10 DEG list (Table [83]1), including YJR005C-A, YOR382W, and YMR058W, showed strong upregulation. This reinforces the transcriptomic signature induced by Hub1 overexpression and supports the robustness of the differential expression analysis. Table 1. Top 10 most significantly differentially expressed genes (DEGs) in Hub1-OE vs. control samples Gene ID log₂ Fold Change (log2FC) Base Mean Adjusted p-value (padj)* YDR533C + 2.26 7,674 0 YOR382W + 4.11 10,351 0 YOR383C + 3.55 17,547 0 YMR059W + 2.70 4,820 0 YMR058W + 2.55 7,802 0 YJR005C-A + 4.97 8,186 0 YER145C + 2.26 11,851 0 YEL065W + 3.77 5,780 0 YEL064C + 3.16 5,241 0 YHL047C + 2.38 5,876 0 [84]Open in a new tab *Genes are ranked based on statistical significance (adjusted p-value) and magnitude of expression change Fig. 2. [85]Fig. 2 [86]Open in a new tab Heatmap of the top 30 differentially expressed genes in Hub1-OE and control samples. Each row represents a gene and each column a sample. Expression values were normalized and scaled (Z-scores). Genes were selected based on adjusted p-value and log₂ fold change from DESeq2 analysis Principle component analysis (PCA) To evaluate the global transcriptional variance and assess sample clustering, principal component analysis (PCA) was conducted on variance-stabilized transformed (VST) gene expression data derived from the DESeq2 pipeline. The PCA plot (Fig. [87]S1) demonstrates a clear separation between the Hub1-overexpressing (Hub1-OE) and control (xrn1 rat1-1) groups along the first principal component (PC1), which explains 98% of the total variance. The second principal component (PC2) accounts for only 1% of the variance, indicating that Hub1 overexpression is the predominant source of transcriptional variation. Biological replicates within each group clustered tightly, suggesting high experimental reproducibility and minimal batch or technical variability. The distinct separation observed along PC1 reinforces the conclusion that Hub1 overexpression induces a substantial and consistent shift in the yeast transcriptome. Functional enrichment and pathway analysis of DEGs To gain functional insights into the biological roles affected by Hub1 overexpression, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using differentially expressed genes (DEGs; adjusted p-value < 0.05). GO enrichment was conducted across three major domains: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). As shown in Fig. [88]3, the top enriched BP terms included proteolysis, phosphorus metabolic process, and small molecule biosynthetic process, indicating increased activity in protein degradation and metabolic regulation. In the CC category, enriched terms such as vacuolar membrane and endoplasmic reticulum subcompartment suggest that Hub1 influences membrane-associated organelle function. Within the MF domain, terms like oxidoreductase activity and transmembrane transporter activity were prominent, reflecting changes in redox biology and molecular transport mechanisms. Fig. 3. [89]Fig. 3 [90]Open in a new tab GO enrichment analysis of differentially expressed genes (DEGs) upon Hub1 overexpression in Saccharomyces cerevisiae.Dot plots display the top 10 significantly enriched GO terms for Biological Process (left), Cellular Component (middle), and Molecular Function (right) To further explore pathway-level effects, KEGG pathway enrichment was performed separately for upregulated and downregulated genes using the clusterProfiler package [[91]35]. As shown in Fig. [92]4, activated pathways in Hub1-OE samples included biosynthesis of secondary metabolites, carbon metabolism, and purine metabolism, indicating upregulation of biosynthetic and metabolic functions. Conversely, suppressed pathways such as ribosome, cell cycle - yeast, and steroid biosynthesis suggest downregulation of growth-related and translational machinery. These findings demonstrate that Hub1 overexpression induces a major transcriptomic shift favoring metabolic reprogramming while attenuating cell cycle and protein synthesis processes. Fig. 4. [93]Fig. 4 [94]Open in a new tab KEGG pathway enrichment analysis of DEGs in Hub1-overexpressing yeast. Dot plot showing the top 15 enriched KEGG pathways among upregulated (left) and downregulated (right) gene sets Gene set enrichment analysis To gain insights into the systemic effects of Hub1 overexpression, Gene Set Enrichment Analysis (GSEA) was performed using a ranked gene list based on DESeq2-derived log2 fold changes. This analysis aimed to identify whether predefined gene sets (KEGG pathways) were enriched among the genes most strongly affected by Hub1 overexpression. As shown in Fig. [95]5, several biologically relevant pathways exhibited notable enrichment patterns. The p53 signaling pathway and Endocytosis were modestly enriched in upregulated genes, with normalized enrichment scores (NES) of 1.255 and 1.342, respectively, suggesting an association with stress response and membrane trafficking. Ubiquitin-mediated proteolysis showed a nearly neutral enrichment (NES = 0.001), while the Cell Cycle pathway displayed a modest negative enrichment (NES = − 0.692), indicating potential transcriptional repression of cell division–related genes under Hub1-OE. These enrichment dynamics support the hypothesis that Hub1 modulates global transcriptional programs, including proteostasis, stress signaling, and cell cycle control, thereby reshaping the physiological state of Saccharomyces cerevisiae. Fig. 5. [96]Fig. 5 [97]Open in a new tab Gene Set Enrichment Analysis (GSEA) of four representative KEGG pathways under Hub1 overexpression.Each panel shows an enrichment plot for one KEGG pathway: A p53 signaling, B Endocytosis, C Ubiquitin-mediated proteolysis, and D Cell Cycle. Green lines indicate enrichment score (ES) curves calculated over the ranked log2 fold change–based gene list Alternative splicing analysis To assess the impact of Hub1 overexpression on pre-mRNA splicing, we performed a comprehensive analysis using rMATS v4.1.2, comparing three biological replicates of Hub1-overexpressing (Hub1-OE) yeast to three controls. Among the five canonical types of alternative splicing (AS) events, only exon skipping (SE) events were detected in the Saccharomyces cerevisiae transcriptome (Fig. [98]6A). A total of seven SE events were identified, consistent with the known low splicing complexity in yeast, where ~ 95% of genes lack introns. No events were detected for alternative 3′ or 5′ splice sites (A3SS, A5SS), mutually exclusive exons (MXE), or retained introns (RI), indicating that Hub1 may selectively modulate exon skipping under the tested conditions. To identify significantly regulated AS events, we constructed a volcano plot based on inclusion level difference (ΔΨ) and false discovery rate (FDR). Among the seven SE events, only one gene DYN2 exhibited statistically significant differential splicing with FDR = 0.0481 and ΔΨ = − 0.036 (Fig. [99]6B). This suggests a modest reduction in exon inclusion in the Hub1-OE condition and a potential target-specific role of Hub1 in modulating splicing precision. Fig. 6. [100]Fig. 6 [101]Open in a new tab Alternative splicing analysis using rMATS identifies a single significant SE event. A Bar plot showing the number of detected events across five alternative splicing types. B Volcano plot of SE events comparing Hub1-OE to control conditions. The DYN2 gene was the only event surpassing the significance threshold (FDR < 0.05, |ΔΨ| > 0.03), indicating reduced exon inclusion. C Heatmap of Z-score normalized exon inclusion levels for all SE events. Samples cluster by condition, with DYN2, RPL7A, and SUS1 showing notable splicing differences. We next visualized exon inclusion patterns across all SE events using a Z-score normalized heatmap (Fig. [102]6C). The heatmap revealed clear sample clustering by condition, with replicates from Hub1-OE and control groups forming distinct groups. While DYN2 was the only statistically significant event, other genes such as RPL7A and SUS1 also showed consistent exon usage differences, hinting at subtle splicing modulation. To validate the splicing alteration in DYN2, we generated a sashimi plot depicting exon junction read coverage and inclusion levels across samples (Fig. [103]7). The exon displayed high inclusion in both groups (93–96% in controls vs. 95–100% in Hub1-OE), supporting the rMATS prediction despite the small ΔΨ value. Although the change is modest, the consistent pattern and statistical support suggest a reproducible splicing effect associated with Hub1 activity. Fig. 7. [104]Fig. 7 [105]Open in a new tab Sashimi plot validating DYN2 exon skipping event in Hub1-overexpressing yeast.Exon junction usage in DYN2 was visualized using sashimi plots generated from BAM files. Read density and junction arcs indicate inclusion of the alternatively spliced exon. Hub1-OE samples show slightly increased exon inclusion (~95–100%) compared to controls (~93–96%), consistent with rMATS output (ΔΨ = –0.036; FDR = 0.0481) Exon usage analysis To evaluate alternative splicing events modulated by Hub1 overexpression, we performed a comprehensive exon usage analysis using DEXSeq. This revealed a distinct pattern of differential exon regulation between Hub1-OE and control yeast strains. The MA plot (Fig. S2A) illustrates the global distribution of exon usage changes. A total of 27 exons were significantly downregulated and 21 exons upregulated (adjusted p-value < 0.05), while 970 exons showed no significant change. Noteworthy among the top differentially spliced regions were mitochondrial genome loci (Q0050 + Q0060 + Q0065 + Q0070 + Q0045 + Q0055), ribosomal subunits (YNL112W, YBR116C + YBR117C), and overlapping non-coding ORFs (YPR160W-A + YPR160W), implicating Hub1 in the regulation of both energy metabolism and ribosome assembly at the splicing level. The hierarchical clustering heatmap (Fig. S2B) of normalized exon counts confirmed consistent and condition-specific patterns across biological replicates. Exons such as E002 and E004 showed upregulation in Hub1-OE samples, while E003 and E001 were downregulated, suggesting a coordinated shift in exon inclusion or exclusion. Together, these findings support the hypothesis that Hub1 overexpression influences splicing fidelity and targets specific exon regions, especially in mitochondrial and ribosomal gene contexts. Weighted gene co-expression network analysis To explore transcriptional co-regulation in response to Hub1 overexpression, we performed Weighted Gene Co-expression Network Analysis (WGCNA) on the filtered expression matrix of 2,178 genes across six RNA-seq samples. After removing low-variance genes, a soft-thresholding power of 12 was selected to construct a signed scale-free network. Hierarchical clustering based on the topological overlap matrix (TOM) identified 61 co-expression modules, each assigned a unique color label (Fig. [106]8A). To assess the relationship between gene modules and experimental conditions, we correlated module eigengenes with a binary trait variable distinguishing Hub1-overexpressing and control samples. The module–trait relationship heatmap revealed several modules with significant correlations (Fig. [107]8B). Notably, the brown module showed the strongest positive association with Hub1 overexpression (r = 0.99, p < 0.001), followed by the pink (r = 0.98, p < 0.001) and turquoise (r = 0.83, p < 0.01) modules. Conversely, the darkolivegreen (r = − 0.91, p < 0.001) and yellow (r = − 0.81, p < 0.001) modules exhibited strong negative correlations, indicating potential downregulation of associated gene networks. To uncover functional insights, we performed KEGG pathway enrichment analysis on genes from the brown module, which was most strongly associated with Hub1 overexpression. The enriched pathways included biosynthesis of secondary metabolites, biosynthesis of amino acids, proteasome, and carbon metabolism (adjusted p < 0.01), reflecting a transcriptional shift toward enhanced metabolic and proteolytic activity (Fig. [108]8C). These results suggest that Hub1 overexpression may promote biosynthetic and growth-related programs while repressing mitochondrial energy production, consistent with the observed suppression of oxidative phosphorylation pathways in negatively correlated modules. Fig. 8. [109]Fig. 8 [110]Open in a new tab Weighted gene co-expression network analysis (WGCNA) of Hub1-OE transcriptome. A Cluster dendrogram of genes based on topological overlap matrix (TOM). Each leaf represents a gene, and branches are grouped into co-expression modules labeled by unique colors. B Heatmap of Pearson correlation coefficients between module eigengenes and experimental trait (Hub1 overexpression vs. control). Red and blue indicate positive and negative correlations, respectively. C KEGG pathway enrichment dot plot for the brown module. Dot size reflects the number of genes; color intensity indicates adjustedp-value. Top enriched pathways include biosynthetic and proteasome-related functions. Isoform-level quantification analysis To explore isoform-specific expression changes associated with Hub1 overexpression, transcript-level quantification was performed using Salmon and normalized to TPM values. Among the analyzed transcripts, DYN2 (YDL007W) a gene encoding a dynamin-like GTPase involved in microtubule organization and endocytosis was found to be significantly upregulated. As shown in Fig.[111]9, the mean TPM of DYN2 increased from 100.0 ± 0.8 in control samples to 117.1 ± 0.5 in Hub1-overexpressing (Hub1OE) samples, corresponding to a 17.1% increase (log₂FC = 0.228). This difference was statistically significant (p = 5.74 × 10⁻⁶, two-tailed unpaired t-test), indicating robust transcriptional activation of DYN2 under Hub1-OE conditions. Both the bar plot (Fig.[112]9A) and box plot (Fig.[113]9B) illustrate the consistent upregulation with low intra-group variability (CV < 1%). These results suggest a functional link between Hub1 expression and DYN2 transcriptional regulation, potentially implicating Hub1 in pathways related to cytoskeletal dynamics and membrane trafficking. Fig. 9. [114]Fig. 9 [115]Open in a new tab DYN2 is significantly upregulated in response to Hub1 overexpression. A Bar plot of DYN2 transcript abundance (TPM) in control vs. Hub1OE samples. Error bars indicate standard deviation (n = 3). B Box plot of DYN2 TPM values across replicates. Tight clustering of values indicates high technical reproducibility. Asterisks (***) indicate statistical significance at p < 0.001 Splice site motif analysis To investigate whether Hub1 preferentially regulates exons with weak splice site signals, we performed motif strength analysis using MaxEntScan on differentially used exons, including the Hub1-responsive DYN2 exon. Splice site sequences (9 bp for 5′ donor, 23 bp for 3′ acceptor) were extracted from the Saccharomyces cerevisiae genome and evaluated using MaxEntScan scoring models. The 5′ donor site of the DYN2 exon scored − 18.32, which was significantly lower than the average canonical 5′ donor site score (mean = − 13.88, SD = 3.20), with a p-value of 0.03 (Fig. [116]10). This suggests that the DYN2 exon has a notably weak donor site. The 3′ acceptor site of DYN2 scored − 5.52, also weaker than canonical scores (mean = − 7.88, SD = 3.15), though this difference did not reach statistical significance (p = 0.12). In addition to DYN2, we evaluated five other Hub1-regulated exons (YDL246C_E1, YDL243C_E1, YDR387C_E1, YDL094C_E1, and YDR438W_E1). All exhibited low MaxEnt scores at either the donor or acceptor sites (Table [117]S1). These findings support a model in which Hub1 facilitates the recognition of non-consensus or suboptimal splice sites, consistent with its known role in modulating spliceosome flexibility. Fig. 10. [118]Fig. 10 [119]Open in a new tab MaxEntScan splice site strength scores for Hub1-regulated exons. Boxplots show distributions of canonical yeast 5′ donor and 3′ acceptor splice site scores (gray). Red dots indicate the DYN2 exon’s donor (−18.32) and acceptor (−5.52) scores Protein- protein interaction analysis To explore the broader regulatory context of Hub1 (YPL254W), we constructed a protein–protein interaction (PPI) network using high-confidence STRING interactions (combined score > 400). This network focused on first-degree interactors of Hub1 and incorporated expression data from our RNA-seq analysis. The resulting subnetwork revealed 35 direct interactors (Fig. [120]11), many of which are differentially expressed under Hub1 overexpression. Notably, YDR224C and YLL039C were among the top upregulated interactors, with log₂ fold changes of + 0.69 and + 0.65, respectively, while YPL235W and YJL076W showed modest downregulation (log₂FC = − 0.33 and − 0.31). Several interactors, including YBR111W-A (a homolog of SNU66) and YHR099W (a homolog of SPP381), are known components of the spliceosomal machinery, reinforcing Hub1’s proposed role in splicing fidelity and exon definition. In addition, interactors like YDR145W, a chromatin remodeling factor, suggest a multi-layered role of Hub1 in transcriptional and post-transcriptional gene regulation. These findings highlight that Hub1 operates within a functionally coherent network, potentially linking RNA processing with chromatin state, and support a mechanistic model where Hub1 modulates non-consensus splicing through interaction with key spliceosomal and regulatory proteins. Fig. 11. [121]Fig. 11 [122]Open in a new tab Protein-protein interaction (PPI) network of Hub1 (YPL254W) and its first-degree interactors. The network displays Hub1 (YPL254W) at the center along with 35 first-degree interacting proteins identified from the STRING database (combined score > 400). Nodes represent proteins, and edges indicate predicted interactions based on experimental and computational evidence. Node color reflects log₂ fold change in expression between Hub1-overexpressing and control cells, with red indicating higher expression and blue indicating lower expression Discussion The integrative transcriptomic and splicing analysis presented in this study elucidates the multifaceted role of Hub1 overexpression in Saccharomyces cerevisiae, revealing its profound impact on both gene expression and alternative splicing. By leveraging RNA-seq data from the [123]GSE84215 dataset and employing a comprehensive bioinformatics pipeline, we identified 3,915 differentially expressed genes (DEGs) and seven exon skipping events, with a significant alteration in the DYN2 gene (FDR = 0.0481, ΔPSI = − 0.036) (Figs. [124]6B and [125]7). These findings underscore Hub1’s dual role as a regulator of transcriptional reprogramming and splicing fidelity, particularly in modulating non-consensus splice sites, and provide novel insights into its contributions to cellular homeostasis and stress adaptation. The extensive transcriptional changes observed, with 1,964 upregulated and 1,951 downregulated genes (Fig. [126]1A; Table [127]1), highlight Hub1’s capacity to orchestrate a balanced yet widespread shift in the yeast transcriptome. The top differentially expressed genes, such as YJR005C-A (log2FC = + 4.97) and YOR382W (log2FC = + 4.11) (Table [128]1), are implicated in RNA processing and stress response, aligning with Hub1’s known functions in RNA metabolism and cellular adaptation [[129]16]. The principal component analysis further confirmed the robustness of these changes, with PC1 explaining 98% of the transcriptional variance and clearly separating Hub1-overexpressing samples from controls (Fig. [130]S1). This suggests that Hub1 overexpression is a dominant driver of global gene expression dynamics, consistent with its role as a ubiquitin-like modifier that interacts non-covalently with regulatory complexes [[131]38]. The heatmap analysis of the top 30 differentially expressed genes (Fig. [132]2) demonstrated clear hierarchical clustering between experimental conditions, with all biological replicates showing consistent expression patterns. This reproducibility across replicates strengthens confidence in the observed transcriptional signature and supports the conclusion that Hub1 overexpression induces systematic changes in gene expression rather than stochastic effects. Functional enrichment analyses provided deeper insights into the biological processes modulated by Hub1. Gene Ontology analysis revealed enriched terms such as proteolysis, phosphorus metabolic process, and vacuolar membrane function (Fig. [133]3), indicating that Hub1 influences protein degradation and organelle-associated processes. Similarly, KEGG pathway enrichment highlighted the upregulation of biosynthetic pathways (e.g., secondary metabolites, carbon metabolism) and downregulation of growth-related pathways (e.g., ribosome, cell cycle) (Fig. [134]4). These findings suggest that Hub1 overexpression induces metabolic reprogramming, potentially redirecting cellular resources toward stress response and biosynthetic activity while attenuating proliferation. Gene Set Enrichment Analysis (GSEA) further supported this metabolic shift, with modest enrichment of stress-related pathways like p53 signaling (NES = 1.255) and endocytosis (NES = 1.342), alongside negative enrichment of cell cycle genes (NES = − 0.692) (Fig. [135]5). These patterns align with prior reports linking Hub1 to stress adaptation and proteostasis [[136]9, [137]38], and they suggest a coordinated regulatory program that balances cellular growth and survival under altered Hub1 levels. The downregulation of ribosome biogenesis and cell cycle pathways, combined with upregulation of stress response mechanisms, indicates that Hub1 overexpression may trigger a cellular state reminiscent of environmental stress adaptation [[138]39]. The alternative splicing analysis revealed a targeted effect of Hub1 on exon skipping, with DYN2 emerging as the only statistically significant event (FDR = 0.0481, ΔPSI = − 0.036) among seven detected exon skipping events (Fig. [139]6A and B). The modest reduction in DYN2 exon inclusion, validated by sashimi plots (Fig. [140]7), underscores Hub1’s role in modulating splicing precision, particularly for genes with weak splice sites [[141]40]. MaxEntScan analysis confirmed that the 5′ splice site of the DYN2 exon is significantly weaker than canonical yeast splice sites (score = − 18.32, p = 0.03) (Fig. [142]10), consistent with Hub1’s established function in enhancing non-consensus splice site recognition [[143]6]. The complementary DEXSeq analysis identified 48 differentially used exons (Figure S2A), with hierarchical clustering revealing consistent condition-specific patterns (Fig. S2B). Notably, several mitochondrial genome loci and ribosomal subunits showed differential exon usage, suggesting that Hub1’s splicing effects extend beyond nuclear genes to organellar transcripts. This broader splicing impact aligns with Hub1’s known role in facilitating splicing of genes with suboptimal splice sites [[144]10]. The splice site strength analysis of Hub1-regulated exons (Table [145]S1) revealed that most targeted exons contain weak donor or acceptor sites, supporting a model where Hub1 facilitates splicing of suboptimal sites through interactions with spliceosomal components like Snu66 [[146]41]. Given that only ~ 5% of S. cerevisiae genes contain introns [[147]42], the detection of even a small number of splicing events is notable and suggests that Hub1’s splicing role is highly selective yet functionally important. The significant splicing alteration in DYN2, a gene encoding a dynamin-like GTPase involved in microtubule organization and endocytosis [[148]43], points to a functional link between Hub1 activity and cytoskeletal dynamics. The observed 17.1% increase in DYN2 transcript abundance (log2FC = 0.228, p = 5.74 × 10⁻⁶) (Fig. [149]9) suggests that Hub1 not only modulates DYN2 splicing but also enhances its transcriptional output. This dual regulation may amplify DYN2’s role in membrane trafficking and cytoskeletal remodeling, processes critical for stress adaptation and cellular homeostasis [[150]44]. The enrichment of endocytosis pathways in GSEA (Fig. [151]5) further supports this connection, suggesting that Hub1 may coordinate splicing and expression changes to enhance endocytic pathways under stress conditions. Weighted Gene Co-expression Network Analysis (WGCNA) provided a systems-level perspective on Hub1’s regulatory influence, identifying 61 co-expression modules, with the brown module showing the strongest correlation with Hub1 overexpression (r = 0.99, p < 0.001) (Fig. [152]8B). This module’s enrichment in biosynthetic and proteasome-related pathways (Fig. [153]8C) reinforces the notion that Hub1 drives metabolic and proteolytic reprogramming. The negative correlation of modules like darkolivegreen and yellow with Hub1 overexpression suggests a coordinated suppression of specific cellular programs, potentially including mitochondrial energy production as part of a broader stress response strategy [[154]6]. The tight clustering of genes within co-expression modules (Fig. [155]8A) and their strong correlation with experimental conditions demonstrate that Hub1 overexpression affects not just individual genes but entire regulatory networks [[156]45]. This systems-level impact underscores Hub1’s role as a master regulator capable of orchestrating complex transcriptional programs [[157]46]. The protein-protein interaction (PPI) network analysis illuminated Hub1’s regulatory context, identifying 35 direct interactors, including spliceosomal components like YBR111W-A (SNU66 homolog) and YHR099W (SPP381 homolog) (Fig. [158]11). The differential expression of these interactors, such as the upregulation of YDR224C (log2FC = + 0.69) and YLL039C (log2FC = + 0.65), suggests that Hub1 operates within a functionally coherent network that integrates splicing and transcriptional regulation [[159]6]. The presence of chromatin remodeling factors among the interactors further indicates that Hub1’s effects may extend to epigenetic regulation, consistent with its reported roles in nuclear envelope maintenance and chromatin remodeling [[160]47]. This multi-layered regulatory network underscores Hub1’s versatility as a ubiquitin-like modifier and its capacity to bridge post-transcriptional and transcriptional control [[161]5]. The coordinated expression changes among Hub1 interactors suggest that Hub1 overexpression may amplify its regulatory effects through positive feedback loops within its interaction network. The findings of this study have significant implications for understanding Hub1’s role in eukaryotic gene regulation. The conservation of Hub1 across species suggests that its regulatory mechanisms in S. cerevisiae may inform studies in higher eukaryotes, where Hub1 (UBL5) is implicated in splicing and stress response [[162]48]. The specific modulation of DYN2 splicing and expression highlights a potential regulatory axis linking Hub1 to cytoskeletal dynamics and stress adaptation, which warrants further investigation in other model systems. The metabolic reprogramming observed under Hub1 overexpression, characterized by enhanced biosynthetic activity and reduced growth-related processes, suggests that Hub1 may function as a molecular switch that adjusts cellular priorities in response to stress or environmental changes [[163]49, [164]50]. This regulatory mechanism could be particularly important in organisms facing fluctuating environmental conditions, where rapid adjustment of gene expression programs is crucial for survival [[165]51]. While this study provides comprehensive insights into Hub1’s regulatory role, several limitations must be acknowledged. The analysis relies on a single experimental dataset, potentially limiting the generalizability of findings across different growth conditions or stress scenarios [[166]52]. The lack of experimental validation through RT-PCR or functional assays represents a significant limitation [[167]53]. Future studies should prioritize experimental confirmation of the identified splicing changes and functional characterization of their cellular consequences [[168]54]. Additionally, investigation of Hub1’s regulatory mechanisms under diverse stress conditions would provide deeper insights into its adaptive functions [[169]55]. The modest effect size of the DYN2 splicing change, while statistically significant, raises questions about biological relevance that can only be addressed through functional validation [[170]56]. Complementation studies using splicing-defective DYN2 variants could determine whether subtle splicing changes translate to measurable phenotypic effects. Conclusion This study establishes Hub1 as a dual regulator of transcriptional reprogramming and alternative splicing in Saccharomyces cerevisiae, coordinating metabolic adaptation and stress response through integrated regulatory mechanisms. The extensive transcriptional changes (3,915 DEGs) coupled with selective splicing modulation demonstrate Hub1’s capacity to orchestrate complex cellular reprogramming while maintaining homeostatic balance. The identification of DYN2 as a primary splicing target, characterized by weak splice sites and dual transcriptional-splicing regulation, exemplifies Hub1’s role in expanding regulatory precision through multiple control layers. The systems-level analysis reveals Hub1’s integration within highly correlated gene networks and protein interaction complexes, positioning it as a central node in cellular regulatory hierarchies. These findings advance our understanding of ubiquitin-like protein function beyond traditional protein modification roles, revealing Hub1’s capacity to integrate transcriptional and post-transcriptional regulation. The conservation of Hub1 across eukaryotes suggests that these regulatory mechanisms may be fundamental to cellular adaptation strategies, with potential implications for understanding disease processes involving splicing dysregulation. Supplementary Information [171]12864_2025_12006_MOESM1_ESM.docx^ (230.6KB, docx) Additional file 1: Supplementary Figures and Tables. Contains supplementary figures S1-S2 and supplementary table S1 with detailed legends and methodological descriptions. Acknowledgements