OBSCN undergoes extensive alternative splicing during human cardiac and skeletal muscle development

Abstract

Background

Highly expressed in skeletal muscles, the gene Obscurin (i.e. OBSCN) has 121 non-overlapping exons and codes for some of the largest known mRNAs in the human genome. Furthermore, it plays an essential role in muscle development and function. Mutations in OBSCN are associated with several hypertrophic cardiomyopathies and muscular disorders. OBSCN undergoes extensive and complex alternative splicing, which is the main reason that its splicing regulation associated with skeletal and cardiac muscle development has not previously been thoroughly studied.

Methods

We analyzed RNA-Seq data from skeletal and cardiac muscles extracted from 44 postnatal individuals and six fetuses. We applied the intron/exon level splicing analysis software IntEREst to study the splicing of OBSCN in the studied samples. The differential splicing analysis was adjusted for batch effects. Our comparisons revealed the splicing variations in OBSCN between the human skeletal and cardiac muscle, as well as between post-natal muscle (skeletal and cardiac) and the pre-natal equivalent muscle.

Results

We detected several splicing regulations located in the 5’end, 3’ end, and the middle of OBSCN that are associated with human cardiac or skeletal muscle development. Many of these alternative splicing events have not previously been reported. Our results also suggest that many of these muscle-development associated splicing events may be regulated by BUB3.

Conclusions

We conclude that the splicing of OBSCN is extensively regulated during the human skeletal/cardiac muscle development. We developed an interactive visualization tool that can be used by clinicians and researchers to study the inclusion of specific OBSCN exons in pre- and postnatal cardiac and skeletal muscles and access the statistics for the differential inclusion of the exons across the studied sample groups. The OBSCN exon inclusion map related to the human cardiac and skeletal muscle development is available at http://psivis.it.helsinki.fi:3838/OBSCN_PSIVIS/. These findings are essential for an accurate pre- and postnatal clinical interpretation of the OBSCN exonic variants.

Background

Discovered about two decades ago, the name of the gene Obscurin (OBSCN) refers to the challenges that were endured by the researchers in its initial detection and characterization [1]. These challenges were mainly caused by the large size and the relatively low abundance of the transcripts in most tissues. With 121 non-overlapping exons, OBSCN is indeed one of the genes that code for the largest mRNAs in the human genome. In human, it is expressed at highest levels in skeletal muscles, however it is also highly expressed in cardiac muscles (Fig. 1). With the relatively high number of exons, we speculate that OBSCN undergoes extensive alternative splicing, especially exon skipping. These events are known to result in different mRNAs that code for a family of proteins characterized as “obscurins” [1, 2]. The two largest OBSCN isoforms are obscurin A and obscurin B [3]. They are both characterized by a long stretch of more than 60 immunoglobulin (i.e. Ig) repeats that are interrupted by a calmodulin-binding IQ motif and several domains such as fibronectin type-III, SRC homology 3, Rho guanine nucleotide exchange factor and pleckstrin homology. The former (i.e., obscurin A), which is also the smaller of the two (i.e. ~ 720 kDa), features a non-modular C-terminal that includes several phosphorylation sites. In contrast, the C-terminus of obscurin B (i.e. ~ 870 kDa) features a fibronectin type-III domain, two additional Ig sites and two serine/threonine type kinase sites (for a comprehensive review see Grogan A., et al.) [3]. Furthermore, other smaller isoforms are reported to be highly abundant in cardiac muscles. Some of these isoforms solely feature an Ig site, a fibronectin type-III site, and one or two kinase domains [3,4,5].

Fig. 1

Expression of OBSCN in various tissues. Log₁₀ scaled expression of OBSCN across different tissues available in GTEX Analysis Release V8. The curves show the density of the expression values. The horizontal line within the black box shows the median. The black box plots extend vertically from the 25th to the 75th percentile. The outliers are expression levels higher or lower than 1.5 time the interquartile range. GTEx Portal (Analysis Release V8) on 09/13/23

Obscurin is associated with neuromuscular functions such as myofibrillogenesis, hypertrophic response and cytoskeletal arrangements [3]. Furthermore, obscurin interacts with the giant sarcomere protein titin [6]. In fact, obscurin/titin disrupting mutations are suggested to be the cause of several hypertrophic cardiomyopathies and muscular disorders [7,8,9]. Despite the important role of OBSCN in the skeletal muscle development and functioning, to our knowledge, the splicing regulation of OBSCN has not yet been studied in detail. In this study, we analyzed 75 pre- and post-natal skeletal muscle and heart (i.e. 45 postnatal muscles, 7 postnatal hearts, 20 fetal muscles and 3 fetal hearts) RNAseq data. The fetal samples have been denoted with “F”, and the muscle and heart samples are referred to as “M” and “H” respectively. However, to prevent confusion with “prenatal” we have labelled the postnatal samples with “A” as most of these samples are from adults. In particular, we studied the well-known isoforms of OBSCN to quantify and compare the inclusion levels of the OBSCN exons across the studied samples. This allowed us to characterize the splicing regulation of OBSCN in cardiac and skeletal muscle development by comparing pre- and post-natal samples. As the role of obscurin in muscle diseases is an emerging topic of study [8, 10], a detailed map of the splicing and the expression of various isoforms of this gene is of great importance for the correct interpretation of the effects of novel variants.

Results

The prenatal analysis consisted of analyzing inhouse RNA-Seq data from fetal skeletal muscles (n = 20) and fetal cardiac muscles (n = 2) from 2 different fetuses (Table 1). Additionally, it included analyzing publicly available skeletal muscle data from a 19 weeks female, skeletal muscle data from a 22 weeks male, cardiac muscle data from a 28 weeks female and cardiac muscle data from a 19 weeks female (for more information see ‘Methods’). For postnatal analysis, RNA-Seq data from an internal cohort of 44 individuals were analyzed (Table 1). We studied six OBSCN isoforms, four of which were curated mRNAs and featured NM REFSEQ IDs, i.e. ENST00000284548 (NM_052843), ENST00000422127 (NM_001098623), ENST00000570156 (NM_001271223) and ENST00000680850 (NM_001386125). Additionally, we analyzed isoforms ENST00000660857 and ENST00000493977 since they collectively featured four additional unique exons. Overall, we extracted 126 unique or 121 non-overlapping exons from these isoforms, out of which three exons were annotated as alternative first and four exons as alternative last. Additionally, the introns upstream of two exons were annotated with alternative 3’ splicing (Table 2).

Table 1 Internal cohort of individuals whose biopsies were collected for RNA-seq

Exon inclusion/skipping

We measured $\Psi$ inclusion levels of the 126 OBSCN exons for the 45 postnatal skeletal muscle, seven postnatal heart, 20 fetal skeletal muscle and three fetal heart samples. Heatmap (coupled with hierarchical clustering using the “Euclidean” distance) and PCA analysis of OBSCN exon inclusion PSI values showed that the samples did not group based on the sex and clinical diagnosis of the studied individuals (Fig. S1A,B, Additional File 1). However, overall, a clear distinction of pre- and postnatal skeletal muscles with a less distinction of pre- and postnatal cardiac muscles was observed (Fig. 2B, and Fig. S2, Additional File 1). We compared the samples in six different ways: the muscle samples to heart samples (i.e. denoted with M), postnatal (i.e. mostly from adult individual) samples to fetal samples (i.e. A), postnatal muscle samples to fetal muscle samples (i.e. AM/FM), postnatal heart samples to fetal heart samples (i.e. AH/FH), postnatal muscle samples to postnatal heart samples (i.e. AM/AH), and fetal samples to fetal heart samples (i.e. FM/FH) (Fig. 2A). We plotted the average $\Psi$ levels to detect loci within OBSCN where differential splicing was detected when fetal samples were compared to postnatal samples or heart samples were compared to skeletal muscles (Fig. 3). Furthermore, we visualized the distribution of the inclusion levels of the OBSCN exons in the studied samples (using box plots), for those exons whose at least two out of the six comparisons (Fig. 2A) produced significant results (FDR < 0.05) with $\Delta \psi>10$ (%)(Fig. 4). It is worth noting that for three of the studied exons (i.e. exons 48, 53 and 56) significant results were achieved in all comparisons except when postnatal hearts were compared to fetal hearts (Fig. 4, Table S1, and Fig. S3, Additional File 1).

Fig. 2

Sample comparisons in our study. A The RNA-Seq data in our study were analyzed for splicing by comparing the muscle samples to heart samples (denoted with M), mostly adult postnatal samples to fetal samples (denoted with A), mostly adult postnatal muscle samples to fetal muscle samples (denoted also with AM/FM), mostly adult postnatal heart samples to fetal heart samples (AH/FH), mostly adult postnatal muscle samples to mostly adult postnatal heart samples (AM/AH), and fetal samples to fetal heart samples (FM/FH). B Scatterplot shows the separation of the studied samples based on OBSCN exon inclusion PSI values, by illustrating PC1 vs PC2 (achieved from PCA analysis). The sample types have been labelled with different shapes and colours

Fig. 3

Inclusion levels of OBSCN exons: The plot illustrates the average $\Psi$ inclusion levels of the unique exons of OBSCN (in the 4 studied sample classes, i.e. postnatal muscles, postnatal hearts, fetal muscles and fetal hearts) and the exons are ordered by their start position (i.e. X- axis). Each dot shows the average of the $\Psi$ values for an OBSCN exon within a sample class. The dots are marked by the exon numbers. For those dots that are too close together to distinguish, the ranges of the exon numbers are stated. The average $\Psi$ measurements related to a sample class are connected via a line. As their $\Psi$ measurements are not accurate (due to the lack of exon-skipping sequence reads), the first and last exons are shown with red triangles and horizontal grey dashed-lines. The variance of the average $\Psi$ values (across the different sample classes) are shown (with a purple line) below, in the figure. Exon 126 (i.e. an alternative last exon) is omitted as its start coordinate is identical to that of exon 125

Fig. 4

Significantly differentially included OBSCN exons: A-O) Ordered by the exon number, the boxplots illustrate the distribution of the $\Psi$ levels of the exons that were detected as significant (FDR < 0.05) in at least three of the comparisons performed in the study. The box plots extend from the 25th to the 75th percentile, and the thick horizontal line represents the median. The whiskers of the box plots show 1.5 times the interquartile range. The outliers are values higher and lower than the interquartile range. P-R Sashimi plots, illustrate the exon-exon junctions observed in the RNAseq data, in regions flanking exons 17, 18, 48, 98 and 126. The samples with the nearest PSI values to the median PSI of the exons are chosen for the sashimi plots. S-T Relative expression levels of the mRNAs that include exons 17 and 18 (S) (based on primers matching junctions Ex16-Ex17 and Ex18-Ex19), and relative expression of the long and short OBSCN isoforms (T) (based on analyzing exon-exon junctions specific to these isoforms) are shown with bar plots. These values were measured by real-time polymerase chain reaction (i.e. RT-qPCR) in the postnatal and fetal skeletal muscles. The exon-exon junction levels have been normalized to the total OBSCN mRNA levels (i.e. inferred by Ex67-Ex68 junction in S and Ex5-Ex6 junction in T). The sample classes include: mostly adult postnatal muscles (AM), mostly adult postnatal hearts (AH), fetal muscles (FM), and fetal hearts (FH). The significant levels in the plots are shown using asterisks: P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***)

Extensive exon inclusion regulation was detected at several loci at the 5’ end (exons 17 and 18), the middle (exons 48–57) and the 3’ end of the gene (exons 97 ad 107) that were associated with skeletal muscle development (Figs. 3, 4). The inclusion (or usage) of exon 17 (FDR(AM/FM) = 0.00104, $\Delta \Psi$(AM/FM) = -19.2) and exon 18 (FDR(AM/FM) = 0.00236, $\Delta \Psi$(AM/FM) = -18.2) were noticeably lower in postnatal muscles compared to fetal muscles (Fig. 4A-B, 4S and Table S2). Interestingly, however, a similar effect was not seen in the postnatal heart compared to the fetal heart tissues (FDR(AH/FH) > 0.9, -5% < $\Delta \Psi (AH/FH)$ < 0%). This leads us to believe that the inclusion of these exons is regulated specifically during skeletal muscle development (and not during heart development). As a result of this exon inclusion decrease, upregulation of the canonical exon junctions chr1:228243458–228246528 (connecting exons 15 and 18) and chr1:228244571–228256651 (connecting exons 16 and 20) were observed in postnatal muscles compared to fetal muscles (P(AM > FM) = 1e − 04, 0.0243% ≤ $\Delta \Psi$ EJ(AM/FM) ≤ 0.0485%) (Fig. 5 A-B, Table S3).

Fig. 5

Normalized exon-exon junction levels of non-consecutive exons. Boxplots, illustrating the distribution of the normalized canonical (A-G) and non-canonical (H-M) junction levels of non-consecutive exons. If the exon-flanking 5’ or 3’ splice site is included in the reference (i.e. GENCODE) the name of the corresponding exon begins with “EX”; otherwise it starts with “ex”. The sample classes include: mostly adult postnatal muscles (AM), mostly adult postnatal hearts (AH), fetal muscles (FM), and fetal hearts (FH). The p-value and $\Delta EJ$ values for two sets of comparisons are listed below the box plots: postnatal muscle vs fetal muscle (AM > FM), and postnatal heart vs fetal heart (AH > FH). The Jonckheere Terpstra method was used to test the order and extract the significant results. The significant levels are shown using asterisks: P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***). The box plots extend from the 25th to the 75th percentile, and the thick horizontal line represents the median. The whiskers of the boxplots show 1.5 times the interquartile range. The outliers are values higher and lower than the interquartile range

Located in the central region of the OBSCN gene, exons 48–56 were included significantly less in postnatal muscles than in fetal muscles (FDR(AM/FM) < 0.05, $\Delta \psi$(AM/FM) < -10%) (Fig. 4). Although the inclusion levels of these exons were mostly lower in postnatal heart samples than in fetal heart samples (except for exon 52), the effects were milder and the false discovery rates were not significant (FDR(AH/FH) > 0.1, $\Delta \psi$, (AH/FH) < -10%) (Fig. 4, Table S1). We believe that the reason for observing a milder effect in heart is the small size of the fetal heart samples, as the P‐value for some of these effects are less than 0.05 even though their FDR values are not (Fig. 4 K-M). Furthermore, the inclusion levels for most of these exons (i.e. exons 48, 49, 52–56, as well as exons 57 and 58) were significantly higher in human muscle samples compared to human heart samples (FDR(M) < 0.05, $\Delta \Psi$(M) > 30%), suggesting that the detection of exon inclusion variations in the cardiac muscles are technically more challenging and require more sequence reads and biological replicates (Fig. 4). Concurrent to these findings, we also noticed significant increase of several canonical as well as a few noncanonical exon junctions in postnatal muscle samples compared to fetal muscle samples (Fig. 5). The upregulated canonical exon junctions were chr1:228288888 − 228292523 (connecting exons 47 and 49), chr1:228288888 − 228293353 (connecting exons 47 and 50), and chr1:228292161 − 228294152 (connecting exons 48 and 51) (Table S3). The upregulated non-canonical exon junctions were chr1:228293518 − 228294318 (overlapping exons 50 and 51), chr1:228295020 − 228300123 (overlapping exons 52 and 55), and chr1:228300020 − 228303814 (overlapping exon 56) (Table S4).

We developed an interactive visualization tool for the exon inclusion $\Psi$ values using the R Shiny package [11] which is available at http://psivis.it.helsinki.fi:3838/OBSCN_PSIVIS/. The software allows the users to zoom into more precise regions within the OBSCN gene to view the distribution of the inclusion levels (i.e. $\Psi )$ of the exons of interest and the measured statistics.

Alternative first/last exons and alternative 3’ splicing

We measured $\Psi$ values of the four alternative final exons (Table 2) and two alternative first exons. It is worth noting that the results for exons 126 and 125 (i.e. alternative last exons) were reported together as their differences are minor (Table S5). Our results showed that exon 98 (i.e. exon 95 of the meta-transcript) was included significantly less in the mRNAs of postnatal muscles compared to fetal muscles (FDR(AM/FM) = 3.574e − 19, $\Delta \Psi$(AM/FM) = − 36.74%) (Fig. 6). In contrast, exon 125 or 126 were included significantly more (FDR(AM/FM) = 3.419e − 19, $\Delta \Psi$(AM/FM) = 36.69%) in the mRNAs of postnatal muscles compared to fetal muscles (Fig. 6). As a consequence, the skipping of exon 97 (together with exon 98) was significantly upregulated in the postnatal muscles compared to fetal muscles (P(AM > FM) = 1e − 04, $\Delta$ EJ(AM/FM) = 0.34%) (Figs. 5J, 4T). These findings, together with the real-time polymerase chain reaction (RT-qPCR) results, indicate a higher abundance of the longer isoform obscurin-B in the postnatal skeletal muscles despite the higher abundance of the shorter isoform obscurin-A in the fetal skeletal muscles (Fig. 4T, Table S6).

Table 2 The studied OBSCN isoforms and exons

Fig. 6

The inclusion of the alternative last exons in OBSCN mRNAs. The boxplots, illustrate the distribution of the $\Psi$ levels of the alternative last exons. The sample classes for which the $\Psi$ values are shown are: muscle vs heart (M), mostly adult postnatal vs fetal (A), postnatal muscle vs fetal muscle samples (AM/FM), postnatal heart vs fetal heart samples (AH/FH), postnatal muscle vs postnatal heart samples (AM/AH), and fetal muscle vs fetal heart samples (FM/FH). The significant levels are shown using asterisks: P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***). They are coloured in red if $\left|\Delta \psi \right|\ge 10$. The box plots extend from the 25th to the 75th percentile, and the thick horizontal line represents the median. The whiskers of the boxplots show 1.5 times the interquartile range. The outliers are values higher and lower than the interquartile range. Since their start coordinates are identical, the $\Psi$ values for the exons 125 (i.e. 121a of meta-transcript, chr1: 228378618–228378874) and 126 (i.e. 121 of meta-transcript, chr1: 228378618–228378876) are indistinguishable, therefore reported together

We also studied the alternative 3’ splicing related to the exons 122 and 123 (i.e. 119 and 119a of meta-transcript, respectively) (Table 2). The inclusion levels of exon 123 (i.e. 119a of meta-transcript) were very low and the upstream intron was rarely spliced across our studied samples, suggesting that almost all mRNAs in our samples included the alternative exon 122 (i.e. 119 of meta-transcript) (Fig. S4, Additional File 1).

The affected protein domains

The exons 17 and 18 that were frequently skipped in the adult skeletal muscle samples are known to code for Ig domains (Table 2). Furthermore, the exons 48–56 that were less included in the postnatal skeletal and cardiac muscles compared to the equivalent prenatal samples, also code for Ig domains (Table 2). As mentioned earlier, our results showed higher abundance of the longer OBSCN isoform (e.g. obscurin-B) compared to the shorter isoform (e.g. obscurin-A) in postnatal skeletal muscles, even though the shorter isoform was more abundant in fetal skeletal muscles. Compared to the short isoform (i.e. obscurin-A), the long isoform (i.e. obscurin-B) features an addition fibronectin type-III domain, two additional Ig sites and two serine/threonine type kinase sites. These variations in the domains can change the chemical/physical properties of a protein and ultimately affect its function.

Regulation of OBSCN exon inclusion by the splicing factors

We examined the Spearman (rank) correlation of the expression of the significantly differentially expressed splicing factors (when comparing postnatal to fetal muscle samples) with the inclusion $\Psi$ values of the OBSCN exons that were significantly differentially included (in postnatal muscle vs fetal muscle) across the studied skeletal muscle samples. Several significant correlations (|rho|> 0.4, P < 0.05) were detected, e.g. expression of DHX15, THOC1, PRPF1 with inclusion levels of exon 17 and 49 (Fig. 7A-J, Table S7-S10). However, remarkably the expression of BUB3 was significantly correlated with the inclusion levels of most of the significantly differentially included exons (Fig. 7A-I). The BUB3 gene belongs to the budding uninhibited by benomyl (BUB) protein family and is involved in mitosis, aging, carcinogenesis, as well as splicing [12, 13]. Our results suggest the possibility of regulation of OBSCN splicing by BUB3 especially during muscle development.

Fig. 7

Correlation of the OBSCN exon inclusion levels with the expression of several splicing factors. A The rank correlation of the expression of the significantly differentially expressed splicing factors (comparing postnatal muscles vs fetal muscles) with the inclusion levels of the exons that were significantly differentially included (also when comparing postnatal muscles vs fetal muscles) have been shown as a matrix of circles (i.e. correlation plots). The significant (i.e. P < 0.05) correlations have only been shown. The size and the colour of the circles represent the correlation (i.e. rho value) of the corresponding splice factor (labelled on the row) with the corresponding OBSCN exon (labelled on the column). The rho values higher than 0.4 or lower than -0.4 have also been written. B-J For a highly correlated pair (i.e. |rho|> 0.5, P < 0.05), the two lines in the plot show the PSI values of the OBSCN exon, as well as the VST normalized expression levels of the splicing factor (scaled to 100) across the studied samples

Discussion

Alternative splicing plays an essential role in the regulation of gene expression during organ development in mammalians. It is known that throughout the different stages of human life, a great number of genes are differentially spliced, especially in tissues such as brain and heart [14]. Here we studied OBSCN, a gene associated with neuromuscular function that has 121 non-overlapping (or 126 unique) exons and codes for some of the largest mRNAs in the human genome. OBSCN is upregulated in aged skeletal muscle myofiber fragments (e.g. MF-IIsc) and RASA4 + myocytes (Table S10) [15]. However, we studied the splicing regulation of OBSCN during human skeletal and cardiac muscle development. Given the large number of exons in the gene, we hypothesized that it undergoes extensive alternative splicing regulation during muscle and heart development. As a result, we discovered several alternative splicing events in OBSCN associated with skeletal and cardiac muscle development. These mainly included cassette exons and alternative last exon usage events that were significantly differential in the postnatal human skeletal and cardiac muscles, compared to the equivalent prenatal tissues.

The splicing event that was most frequently differential across our pre- and postnatal cardiac, and skeletal muscle samples was exon inclusion (Figs. 3 and 4). The predominance of the exon inclusion (and exon skipping) was not surprising as in previous studies this splicing event has been the most frequent out of all the significant alternative splicing events found in mammals and vertebrates [16]. Also in-line with these findings, exon skipping has been reported as the most frequently regulated event during the development of seven organs (including heart) in six mammals (including human) and a bird [14]. In mice muscles, extensive differential gene expression and alternative splicing has been discovered to occur during the first two weeks after birth, with the vast majority of these alternative splicing events (i.e. 77%) being exon skipping [17].

In this study, we discovered extensive exon inclusion regulation at several loci, at the 5’ end, the middle and the 3’ end of OBSCN gene that are associated with cardiac or skeletal muscle development. Exons 48–56 of OBSCN were significantly less included in RNAs in the postnatal muscles compared to the fetal muscles (FDR(AM/FM) < 0.05, $\Delta \Psi$(AM/FM) < -10%) (Fig. 4C-M). A similar, albeit milder, effect was also seen in cardiac muscles (FDR(AH/FH) > 0.05, $\Delta \Psi$(AH/FH) < 0%) (Fig. 4C-M). It is worth mentioning that from this region of OBSCN (i.e. exons 48–56), exons 48–54 have previously been reported to undergo developmentally dynamic alternative splicing, especially during human heart development [14]. A dynamic alternative splicing event is a splicing event whose PSI changes in the studied tissue (e.g. heart) during human life is greater than 20% [14]. Another exon inclusion regulation, exhibiting strong effects in human muscle development and mild effects in human heart development, was seen in exons 17 and 18 (Fig. 4A,B). Additionally we noticed significant increase in several canonical as well as a few non-canonical exon junction levels in postnatal muscles compared to fetal muscles. To our knowledge the association of these alternative splicing events with human skeletal muscle development has not previously been reported. We also discovered several splicing factors (e.g. DHX15, THOC1, PRPF1, BUB3) whose expression levels were significantly correlated (P < 0.05 and |rho|≥ 0.4) with the inclusion levels of the significantly differentially included exons (when comparing postnatal to fetal muscles) (Fig. 7A-J). Our results suggest that the differential inclusion of the OBSCN exons during skeletal muscle development may be regulated by BUB3. In fact, Bub3 and BuGZ are two essential mitotic regulators that together interact with the splicing machinery in the interphase nucleus [13]. Silencing of either Bub3 or BuGZ has previously shown to enhance exon skipping in Human foreskin fibroblast (i.e. HFF) and ovarian carcinoma TOV21G cell lines [13]. Furthermore, BugZ was not differentially expressed (in postnatal skeletal muscles compared to prenatal skeletal muscles), therefore it is likely that the exon inclusion effects that we report are caused by the differential expression of Bub3 (Table S10). However, it is worth noting that due to scarcity of prior knock-down studies, especially in human muscle samples, a thorough analysis (beyond the scope of this study) is necessary before the precise role of Bub3 in RNA splicing, in muscles of human and other species can be concluded.

In addition to exon skipping, we discovered an alternative last exon usage event that is associated with skeletal muscle development. We discovered that the OBSCN isoform that ends with exon 98 (e.g. obscurin-A) is expressed much higher in fetal skeletal muscles, whereas the larger isoform that skip exons 97 and 98 and end with either exon 125 or exon 126 (e.g. obscurin-B) is more expressed in postnatal skeletal muscles (Fig. 4T and 6). To our knowledge the higher abundance of the longer OBSCN isoform in postnatal muscles compared to fetal muscles has neither been reported earlier.

Almost all the significantly differentially included exons detected in our study code for immunoglobulin domains. Similar to that in titin, obscurin feature long repeats of Ig domains (Table 2) and these tandem Ig domains are mostly coded by individual exons (Table 2) [6, 18]. Therefore, not surprisingly, the most affected (i.e. skipped) exons in our samples code for a complete (not partial) Ig domain (Table 2, Fig. 3). Even though the function of the repeated Ig domains in sarcomeric genes (e.g. OBSCN and TTN) has not thoroughly been studied previously, the extended tandem Ig domains are known to associate with increased elasticity in the isoforms [19]. Furthermore, this has been described as the reason that the sarcomere in skeletal muscle is more elastic comparted to the sarcomere in cardiac muscle [20].

The N-terminus of obscurin interacts with several proteins such as titin, slow myosin binding protein C, and myomesin. Furthermore, the 58th and 59th Ig domains in obscurin (coded by exons 67 and 68 of OBSCN, with ~ 100% of exon inclusion rate) are known to interact with the Z-band of titin, signifying the essential role of OBSCN in myofibrillogenesis [1, 21]. Further structural studies are needed to reveal the precise effects of the significantly differential exon inclusion events that we have described above. However, as the exons 48–58 are distanced from the 5’ end and the titin interacting sites (i.e. exons 67 and 68), the skipping of these exons neither is expected, nor has previously been reported to directly affect the interaction of the obscurin N-terminal with other proteins or to affect the obscurin-titin interaction. The C-terminus of obscurin-A interacts with small Ankyrin 1 (sAnk1) and Ankyrin-B [22, 23]. These interactions are essential to the Ca2 + homeostasis and the assembly of the dystrophin complex, respectively [21]. Furthermore, changes in calcium homeostasis and reduction of dystrophin have both been reported in aged skeletal muscles [24, 25]. However, even though it can be speculated that the downregulation of obscurin-A in adult skeletal muscles (compared to fetal skeletal muscles) that we have described above may contribute to these phenotypes, the connection of these phenotypes to obscurin splicing has not specifically been studied.

Finally, RNA splicing regulation information can assist the researchers and the clinicians to understand the clinical impacts of the exonic variants. As for instance we have recently shown how similar exon usage information for TTN can be used to explain the disease course in nine titinopathy patients [26]. Remarkably, the exon usage information has also been useful in ruling out titinopathy diagnosis for a prenatal case [26]. Therefore, we believe that information related to the OBSCN exon usage and splicing regulation during skeletal/cardiac muscle development, that we have described in detail here, is potentially useful for clinical interpretation of the exonic mutations in OBSCN.

Conclusion

In this study we have described several novel exon skipping events that are associated with human cardiac and skeletal muscle development. Additionally, we discovered an alternative final exon usage event associated with human skeletal muscle development. This information allows us to understand the regulation of OBSCN splicing during human muscle and heart development. Furthermore, the data is essential for clinical and prognostic interpretation of the OBSCN exonic variants and understanding the effects of these variants on the protein expression in different stages of life.

Our study is strengthened by the thoroughness of the analysis and the support for P-values (and FDRs) that describe how significantly differential the alternative splicing events are during human skeletal and cardiac muscle development. The study is however reliant on the analysis of RNAseq data from a limited selection of muscle types (e.g. tibialis anterior and vastus lateralis from the studied postnatal individuals). Therefore, it can be improved by including RNAseq data from more samples and from a more diverse types of muscle tissues in the analysis. Finally, we have developed an interactive visualization tool (using the shiny R package) that can easily be used by the clinicians to check the inclusion level of each OBSCN exon during skeletal and cardiac muscle development. The interactive R shiny application is available at http://psivis.it.helsinki.fi:3838/OBSCN_PSIVIS/.

Methods

In-house data

For prenatal analysis, a trained fetal pathologist collected fetal skeletal muscles (n = 20) and fetal cardiac muscles (n = 2) from 2 different fetuses, without muscle pathology, obtained from voluntary termination of pregnancy (TOP).

For postnatal analysis, we collected sample biopsies from an internal cohort of 44 individuals (Table 1). RNA was extracted with the Qiagen RNeasy Plus Universal Mini Kit (Qiagen, Hilden, Germany) according to the instructions provided by the manufacturer. Total RNA-Seq libraries were prepared using the Illumina Ribo-Zero Plus rRNA Depletion Kit (Illumina, Palo Alto, CA, USA) at the Oxford Genomics Center, Welcome Trust Institute, Oxford, United Kingdom and Novogene. Sequencing was performed using NovaSeq 6000 (Illumina), generating over 80 million 150 bp-long reads per sample.

External data

In addition to our in-house data, we studied four samples from ENCODE with accession IDs: ENCBS067RNA (fetal skeletal muscle tissue, 19 weeks female), ENCBS068RNA (fetal skeletal muscle tissues, 22 weeks male), ENCBS055RNA (fetal heart tissue, 28 weeks female) and ENCBS056RNA (fetal heart tissue, 19 weeks female) [27, 28]. Furthermore, for a view on OBSCN expression across different human tissues, we obtained data from the GTEx Portal (Analysis Release V8) and dbGaP accession number phs000424.v8.p2 on 09/13/23.

RNA-Seq read alignment

The paired RNA-Seq reads were mapped to the Human Genome (GRCh38.p13) using the splice-aware alignment software STAR (V2.7.7a) [29]. The software was run in 2-pass mode and most parameters were set to their default values. For the gene annotation, Gencode.v39 was used (further details available in supplemental methods, Additional file 1).

Splicing analysis and exon inclusion level estimation

The inclusion levels (i.e. PSI or Ψ values) of all unique exons in human genome, including those of OBSCN gene, were measured using the Intron Exon Retention Estimator (IntEREst) R/Bioconductor package (V1.26.1) [30]. IntEREst is a comprehensive RNA-Seq read summarization, differential intron retention and splicing analysis software. It supports tools that measure suitable Ψ values and run statistical differential test for splicing analysis. The inclusion Ψ values were measured for every OBSCN exon. The statistical significance of the increase or decrease of the inclusion levels of OBSCN exons was preformed genome-wide for all exons (except for the first and last exons), however later the results for the OBSCN exons were extracted [31]. The statistical test compared the variation of inclusion of each exon relative to the genome-wide variation observed for inclusion of the studied exons. The analysis was adjusted for possible biases introduced by the different sequencing batches by including this parameter as a covariate in the design model of the statistical tests (further details available in supplemental methods, Additional file 1).

In addition to exon skipping/inclusion, we ran a similar analysis for the inclusion of the alternative first and last exons, as well as the only case of alternative 3’ splicing in OBSCN (i.e. affecting exons 122 and 123, or 119 and 119a from the meta transcript) (Table 2). All the first/last exons of the studied genes (including OBSCN) were extracted using biomaRt [32]. All P-values were adjusted for multiple testing, using the Benjamini–Hochberg method [33]. An FDR < 0.05 cutoff was used to extract the significant results.

Finally, differential gene expression analysis was performed using DESeq2. The analysis was adjusted for the batch effects. The rank correlation of the VST normalized expression levels of the splicing factors (whose IDs were extracted from other studies [34,35,36]) with the exon inclusion $\Psi$ values was also measured using the Spearman method (further details available in supplemental methods, Additional file 1).

Real-time polymerase chain reaction (RT-qPCR) validation

RNA was extracted from two adult muscle samples and two fetal muscle samples using the Qiagen RNeasy Plus Universal Mini Kit (Qiagen, Hilden, Germany) and according to the instructions provided by the manufacturer. The cDNA synthesis was performed using SuperScript III Reverse Transcriptase (Invitrogen TM) and random primers, according to the protocol provided by the manufacturer. The UCSC In-Silico PCR tool and Primer3web v4.1.0 were used to design primers to target either exon-exon junctions or other regions near the junctions (Table S11). The RT-qPCR assays were performed using the iQ SYBR Green Supermix (BIO-RAD) and 25 nM of each specific primer. Furthermore, three technical replicates were taken into consideration. For the normalization, 18S was used as the reference gene. The final results were calculated using the ΔΔCt method and the relative quantification values were plotted (Fig. 4S-T).