Global alternative splicing landscape of skeletal muscle atrophy induced by hindlimb unloading.

BACKGROUND: Long-term exposure to microgravity will cause skeletal muscle atrophy, which can cause serious harm to astronauts in space travel. Therefore, it is important to explore skeletal muscle atrophy's molecular mechanism for its prevention and treatment. However, as an important regulatory approach of skeletal muscle physiology, the role of alternative splicing in skeletal muscle atrophy, especially skeletal muscle atrophy caused by disuse, is unclear.
METHODS: We established a rat hindlimb unloading model and performed RNA sequencing on soleus muscle, which was seriously atrophied during unloading. Several bioinformatics methods were used to identify alternative splicing events and determine their gene functions.
RESULTS: Many alternative splicing events were found to occur at different time points (12 h, 24 h, 36 h, 3 days, and 7 days) after hindlimb unloading. These differential alternative splicing events mainly occurred in the gene's coding domain sequence region, and 59% of the alternative splicing events caused open reading frameshift. Bioinformatics analysis results showed that genes with different alternative splicing events were enriched in multiple pathways related to muscle atrophy, including the insulin signaling pathway, endocytosis, mitophagy, and ubiquitin-proteasome pathway. Moreover, alternative splicing of several deubiquitinase genes persisted during skeletal muscle atrophy induced by unloading. Additionally, we identified 10 differentially expressed RNA binding proteins during skeletal muscle atrophy induced by unloading, mainly containing Xpo4, Eif4e2, P4ha1, Lrrfip1, Zc3h14, Emg1, Hnrnp h1, Mbnl2, RBfox1, and Mbnl1. Hnrnp h1 and Mbnl2 were significantly downregulated, and RBfox1 and Mbnl1 were significantly upregulated during skeletal muscle atrophy caused by unloading.
CONCLUSIONS: To the best of our knowledge, the present study is the first to propose alternative splicing alterations related to disuse-induced muscle atrophy, emphasizing that alternative splicing is a new focus of attention in the occurrence of muscle atrophy. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

After entering the space orbit, the spacecraft will be in a state of weightlessness. Extended periods of weightlessness or in a microgravity environment restrict astronauts’ movement, resulting in skeletal muscle atrophy. It has been found that after bed rest, physical restriction, or space flight for about 5 days, skeletal muscles will atrophy significantly (1,2). The specific characteristics of skeletal muscle atrophy are as follows: decreased muscle mass, muscle fiber cross-sectional area, and muscle contraction capacity. Further studies have shown that soleus muscle atrophy is the most severe type of atrophic muscle, accompanied by the transition from slow-switch fiber to fast-switch fiber (3,4). The molecular mechanism of skeletal muscle atrophy is complex and still unclear. Therefore, there are no effective drugs or methods to prevent skeletal muscle atrophy caused by weightlessness.

Although the causes of skeletal muscle atrophy are diverse, such as denervation, physical damage, cancer, aging, and disuse, they have similar pathological characteristics and some common activation pathways (5-7). The cause of the atrophic phenotype is the imbalance of protein metabolism, that is, protein synthesis is inhibited and protein degradation is activated (8). Three major protein degradation pathways are activated during muscle atrophy as follows: the ubiquitin-proteasome pathway, autophagy-lysosome pathway, and calpain pathway (9-11). Two studies using proteomics in denervated muscle atrophy model found that a large number of protein expression changes were accompanied by muscle atrophy, ubiquitin proteasome pathway was widely activated (12-14). Our study showed that the ubiquitin-proteasome pathway is activated in the middle and late stages of skeletal muscle atrophy, suggesting that this pathway does not trigger muscle atrophy. In the early stage of skeletal muscle atrophy, oxidative stress, and inflammation-related pathways are relatively active, which may be an important trigger of skeletal muscle atrophy (15,16). Also, numerous studies have shown that a large number of non-coding RNAs are differentially expressed during muscle atrophy. They form a complex competing endogenous RNA network with target genes and participate in the process of muscle atrophy (17-20). Therefore, skeletal muscle atrophy is regulated by diverse events.

However, as an indispensable part of gene expression regulation, the role of alternative splicing in the process of muscle atrophy has attracted less attention. A single gene produces multiple transcripts through alternative splicing, which greatly expands the proteome’s complexity and the diversity of protein functions (21). Alternative splicing can also change mRNA’s stability and the production of non-coding RNA and regulate the expression of a large number of genes (22-24). It has been found that 95% of human multiple exon genes are alternative spliced, and they are specifically regulated under different physiological or pathological conditions (25). There are extensive alternative splicing phenomena in skeletal muscle, which are involved in various physiological processes of muscle cells, such as contraction, aging, development, and regeneration (21,26-30).

Moreover, the imbalance of RNA splicing can also cause pathological changes in skeletal muscle (31), and RNA binding protein (RBP) deletion can directly lead to the loss of skeletal muscle mass (32-34). This suggests that alternative splicing may be a new regulatory mechanism of muscle atrophy. However, there are few studies on alternative splicing events during muscle atrophy.

Hind limbs suspension can simulate weightlessness, and hanging for extended periods can cause a phenotype similar to the clinical muscle atrophy caused by disuse. In the present study, we established a rat hindlimb unloading model. RNA sequencing (RNA-seq) was performed on the soleus muscle, which was severely atrophic during unloading (12 h, 24 h, 36 h, 3 days, and 7 days). We analyzed the dynamic changes and biological functions of the genes that underwent alternative splicing during skeletal muscle atrophy caused by hindlimb unloading and found several genes involved in the ubiquitin-proteasome pathway have alternative splicing phenomenon. We also analyzed the expression of RBPs, which drove the changes of alternative splicing. This can enrich the molecular mechanism of skeletal muscle atrophy and provide a scientific basis for the targeted therapy of muscle atrophy.

We present the following article following the ARRIVE reporting checklist (available at http://dx.doi.org/10.21037/atm-20-5388).

Methods

Establishment of rat models of hindlimb unloading-induced atrophy

Eighteen 8-week-old male Sprague-Dawley rats were used in the present study. The experiments were done following the recommendations of and approved by the Institutional Animal Care and Use Committee of Nantong University, China (No. S20200312-003). The rats’ hindlimbs were just off the ground, and the tilt angle between the body and the ground was approximately 30 degrees, following the internationally recognized tail unloading method of rodents (35). At this time, the rats could move freely with access to food in the cage. The animals were raised in individual cages at 23±2 °C, and with 12 h of light per day. According to the randomization principle, the rats were killed by cervical dislocation at 6 time points (0 h, 12 h, 24 h, 36 h, 3 days, 7 days), and the soleus muscles were taken. After cutting off excess fat and fascia, the samples were quickly immersed in liquid nitrogen for storage. Normal (0 h) rat soleus muscle was used as the control group, and samples at each time point were obtained from 3 independent repeated experiments (15 suspended and 3 control animals in total).

RNA extraction

Total RNA was extracted using TRIzol (Vazyme, Nanjing, Jiangsu, China). Briefly, the muscle tissue was soaked in 1 mL TRIzol and completely crushed by a homogenizer. Using 1/5 volume of chloroform, RNA was extracted at 12,000 rpm at 4 °C. An equal volume of isopropanol was used to precipitate the RNA, and then the precipitate was washed with 75% ethanol.

cDNA synthesis and library construction

In total, 4 µg RNA was purified and enriched by magnetic beads with oligo(dT). The mRNA was broken into short fragments by adding interruption of reagent. Using interrupted mRNA as a template, the first-strand cDNA was synthesized by using 6-base random primers. Double-stranded cDNA was synthesized and purified. Purified double-stranded cDNA was then subjected to end repair. Poly(A) tailing was added and connected to a sequencing adapter. Finally, polymerase chain reaction amplification was performed. After the Agilent 2100 Bioanalyzer qualified the constructed library, it was sequenced using the Illumina HiSeq X Ten sequencer to generate 125 or 150 bp double-ended data.

Filtering and mapping of sequencing data

The raw data (raw reads) generated by high-throughput sequencing in FASTQ format. To obtain high-quality reads that could be used for subsequent analyses, further quality filtering of raw reads was required. First, Trimmomatic software was used to control quality. Low-quality bases and N bases were filtered out, and finally, high-quality clean reads were obtained. The clean reads were compared using HISAT2 as the reference genome of the species. The software parameters were the default parameters, and the genome mapping rate evaluated samples.

Identification of alternative splicing events

To analyze the differences in alternative splicing events between different samples, rMATS software was used to classify alternative splicing events based on RNA-seq data. Alternative splicing events were divided into the following 5 categories using rMATS: skipped exon (SE), alternative 5' splice site (A5SS), alternative 3' splice site (A3SS), mutually exclusive exon (MEX), and retained intron (RI). In the analysis, junction counts were used for the alternative splicing event detection. The number of reads across the inclusion junction site and across the skipping junction site was calculated. The value of the inclusion level (IncLevel) was used to define the degree of alternative splicing. The IncLevel was equal to the number of reads across the inclusion junction site/(number of reads across the inclusion junction site + number of reads across the skipping junction site). The difference in the value of the IncLevel difference between the experimental and control groups was >0.01, and the false discovery rate (FDR) was not >0.05 as the criterion for screening differential alternative splicing events. The ΔIncLevel was equal to IncLevel2-IncLevel1.

Gene function analysis

We performed Gene Ontology (GO) enrichment analysis on differentially alternative spliced genes to describe gene function. The method for the enrichment analysis of the GO function was as follows. All protein coding genes/transcripts were used as the background list. A list of differential protein encoding genes/transcripts was utilized as a candidate list screened from the background list. The hypergeometric distribution test was used to calculate the P value, representing whether the GO function library was significantly enriched in the differential protein coding genes/transcripts list. The P value was corrected by Benjamini-Hochberg multiple tests to obtain FDR. For the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (http://www.genome.jp/kegg/), a hypergeometric distribution test was also used to calculate the significance of differential gene enrichment in each pathway entry. The calculated results could return a significant P value, and a low P value indicated that the differential genes were enriched in the pathway. The STRING database (https://string-db.org/) was used to analyze the interaction network between RBPs.

Statistical analysis

The values of gene expression and the IncLevel alternative splicing events in the present study were calculated from the average of 3 repeated experiments. The relative change represented the difference between each experimental group and the control group, and the negative binomial distribution test was used to test the significance (P value) of the difference between the 2 groups of data.

Results

Differential alternative splicing events during muscle atrophy caused by hindlimb unloading

To identify the alternative splicing events during muscle atrophy caused by hindlimb unloading, we acquired the soleus muscle at 12 h, 24 h, 36 h, 3 days, and 7 days and performed RNA-seq of the poly(A) method to construct the library. We obtained the original base number >7.4 G/sample, the effective base percentage >93.46%, and Q30 >94.2%. At least 97.15% of reads matched the rat genome, suggesting higher reliability of sequencing (Table S1).

Mammals contained the following 5 common alternative splicing modes: SE, A5SS, A3SS, MEX, and RI. Using rMATS software for the data analysis (36,37), we identified a total of 33,505, 33,547, 33,179, 31,966, and 32,906 alternative splicing events at 12 h, 24 h, 36 h, 3 days, and 7 days after hindlimb unloading, which matched the 13,616, 13,707, 13,650, 13,374, and 13,543 genes, respectively (, B). Compared with the control group, according to the absolute value of the IncLevel difference >0.01 and FDR ≤0.05, we identified 2,118, 1,975, 2,132, 2,242, and 2,478 differential alternative splicing events at 5 time points, which matched the 1,762, 1,655, 1,775, 1,849, and 1,946 genes, respectively (, D). The number of differentially spliced genes was almost the same at every time point. This suggests that the responses were early and constant and did not change much over the time course of 7 days. Of these, SE was the main alternative splicing mode, regardless of the event or gene, and accounted for >70%. We then analyzed the changing trends of these differential splicing events at 5 time points. Compared with the control group, the IncLevel difference >0 was defined as splicing activation, representing a longer transcript, and the IncLevel difference <0 was defined as splicing inhibition. The results showed that for the alternative splicing events in the MEX, A5SS, and A3SS modes, the proportions of activation and inhibition were not significantly different. The inhibited alternative splicing events accounted for a relatively large number in the late stage of atrophy for the RI model, reaching 75% (3 days) and 69% (7 days). As the most common alternative splicing mode, SE increased the number of alternative splicing events activated during muscle atrophy ().

Figure 1

Differential alternative splicing events during muscle atrophy caused by hindlimb unloading. Each group of data contained 3 independent samples; normal rat soleus muscle was used as the control group. (A) Total number of alternative splicing events identified during atrophy. (B) Total number of genes corresponding to alternative splicing events. (C) Number of differential alternative splicing events at each time point compared with the control group. (D) Number of differential alternative splicing genes at different time points. (E) Number and proportion of the 5 alternative splicing forms that activated and inhibited alternative splicing events during muscle atrophy. White box represents splicing inhibition, and black box represents splicing activation. A3SS, alternative 3' splicing site; A5SS, alternative 5' splicing site; MEX, mutually exclusive exons; RI, retention intron; SE, skipped exon.

Distribution of differential alternative splicing events on genes and their influence on open reading frames

Alternative splicing occurred in the process of pre-mRNA excision of introns to form mature mRNA. The mRNA contained 3 regions, 5' untranslated region (5'UTR), coding sequence (CDS), and 3' untranslated region (3'UTR). Alternative splicing in different regions may affect gene expression through different mechanisms, such as transcriptional activation, protein-coding, and mRNA stability. Therefore, we analyzed the location of SE form alternative splicing events on the mRNA at 5 time points. The results showed that most alternative splicing occurred in the CDS region. Taking 7 days as an example, the alternative splicing events that occurred in the 5'UTR, 5'UTR-CDS, CDS, CDS-3'UTR, and 3'UTR regions were 183, 233, 985, 170, and 148, respectively, of which the CDS region accounted for 57.3% ().

Figure 2

Distribution of alternative splicing events on genes and their impact on open reading frames. (A) Distribution of differential alternative splicing events on genes at 5 time points: Gene is divided into 5 regions: 5'UTR, CDS-5'UTR, CDS, CDS-3'UTR, and 3'UTR. (B) Number and proportion of frameshift and non-frameshift alternative splicing of genes corresponding to alternative splicing events in the CDS region. White box represents the genes that can be spliced with frameshift alternative splicing, and the black box represents the genes with non-frameshift alternative splicing. CDS, coding sequence; SJ, splicing junction; 3'UTR, 3' untranslated region; 5'UTR, 5' untranslated region.

If the alternative exon's base number was not a multiple of 3, alternative splicing could cause protein frameshift or activate the nonsense-mediated mRNA decay pathway to degrade mRNA (38). The impact on the protein was devastating. If the base number of the alternative exon was a multiple of 3, it would only cause changes in the number and type of several amino acids. In most cases, the impact on protein function was relatively minimal. The number and proportion of genes containing frameshift alternative splicing events in the CDS region are shown in . The findings indicated that the proportion of frameshift alternative splicing genes was higher than that of non-frameshift alternative splicing genes. At the 5 time points, we identified 458, 398, 420, 435, and 489 frameshift alternative splicing genes, respectively ().

Functional analysis of differential alternative splicing genes

To obtain a functional view of differential alternative splicing genes, we selected 2 representative time points as follows: the intermediate time point of 36 h and the late time point of atrophy of 7 days. Because SE is the main alternative splicing mode, we performed a KEGG analysis of genes with alternative differential splicing in the SE model. According to the different biological functions, these differential alternative splicing genes were classified into 6 processes as follows: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. The results showed that the number of genes enriched in each pathway at 7 days was generally greater than those at 36 h, suggesting that skeletal muscle had more active molecular changes at 7 days. The following 4 pathways had the most significant increase in the number of genes: endocrine system, infectious diseases, endocrine, and metabolic diseases, and signal transduction (). We analyzed the Top 20 most significant KEGG pathways, and the findings indicated that 9 pathways were significantly activated at 36 h and 7 days as follows: pancreatic cancer, bacterial invasion of epithelial cells, glucagon signaling pathway, insulin signaling pathway, neurotropin signaling pathway, endocytosis, mitophagy-animal, mRNA surveillance pathway, and propanoate metabolism. Of these, the insulin signaling pathway, endocytosis pathway, and mitophagy-animal pathway have been reported to be associated with muscle atrophy () (39,40). It is worth noting that the mammalian target of rapamycin (mTOR) signaling pathway was enriched to Top20 at 36 h instead of 7 days (by comparing the P value, it is 0.0003 in 36 h and 0.0072 in 7 days, by comparing the enrichment score, it is 2.2257 in 36 h and 1.7689 in 7 days), suggesting that alternative splicing may affect protein synthesis at the late stage of atrophy (41).

Figure 3

Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of differential alternative splicing genes in skipped exon form at 36 h and 7 days after hindlimb unloading. Genes that were differentially spliced at 36 h (A) and 7 days (B) are clustered into 5 biological processes as follows: cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems. At 36 h (C) and 7 days (D) the TOP20 most significant KEGG pathway. Color of the dots represents the P-value, and the size of the dots represents the number of differential alternative splicing genes in the pathway.

Next, we analyzed the GO of the alternative splicing genes with SE mode. The most significant GO terms of Top30 were divided into the following 3 aspects: biological process, molecular function, and cellular component. We first compared the functions of differentially alternative splicing genes at 36 h and 7 days. Overall, the degree of activation at 7 days was generally stronger than at 36 h, which was mainly reflected in the negative regulation of the establishment of the protein localization, regulation of peroxisome organization, protein binding, actin binding, and cadherin binding pathways (). It is worth noting that 1 ubiquitin-related pathway was enriched at 3 days, and 4 ubiquitin-related pathways were enriched at 7 days, suggesting that changes in alternative splicing events may promote significant activation of the ubiquitin-proteasome pathway at 7 days (, B). Next, we analyzed the function of SE mode splicing activation gene and splicing repression gene; the results are shown in (36-h activation), (7-day activation), (36-h repression), and (7-day repression). After analysis, we found that the ubiquitin-proteasome pathway was mainly represented by splicing repressing genes, rather than splicing activating genes. Interestingly, the regulation of RNA splicing and the regulation of transcription by RNA polymerase II were widely activated, suggesting that alternative splicing may cause expression changes in a large number of downstream genes ().

Figure 4

Gene Ontology (GO) analysis of the differential alternative splicing genes in the skipped exon form at 36 h and 7 days after hindlimb unloading. (A) GO analysis of all differential alternative splicing genes at 36 h; (B) GO analysis of all differential alternative splicing genes at 7 days; (C) GO analysis of alternative splicing activator genes at 36 h; (D) GO analysis of alternative splicing activator genes at 7 days; (E) GO analysis of alternative splicing repression genes at 36 h; (F) GO analysis of alternative splicing repression genes at 7 days. The TOP30 most significant term is shown. Green histogram represents the biological process, red represents the molecular function, and blue represents the cellular component.

Alternative splicing continues to affect the ubiquitin–proteasome pathway during muscle atrophy caused by hindlimb unloading

Muscle atrophy is a gradual process and may require certain molecules to continue to function. Therefore, we took the intersection of the differential alternative splicing genes at different time points, and 201 continuously changing genes were identified (). GO analysis of these genes showed that the biological processes of cell-cell adhesion, transcription, DNA-templated, protein deubiquitination, skeletal muscle contraction, skeletal muscle tissue development, actin filament organization, cytoskeleton organization, microtubule cytoskeleton organization, regulation of transcription, DNA-templated, and regulation of muscle contraction were enriched (). Molecular functions of poly(A) binding, protein binding, actin-binding, structural constituent of muscle, cadherin binding involved in cell-cell adhesion, N6-methyladenosine-containing RNA binding, thiol-dependent ubiquitinyl hydrolase activity, tropomyosin binding, lamin binding, and actin filament binding were enriched (). As the ubiquitin-proteasome pathway plays an important role in muscle atrophy, we focused on the splicing changes of 6 genes in protein deubiquitination and thiol-dependent ubiquitinyl hydrolase activity pathway and found that Brca1 associated protein 1 (Bap1), ubiquitin specific peptidase 15 (Usp15), and ubiquitin specific peptidase 54 (Usp54) splicing was continuously suppressed; ubiquitin specific peptidase 2 (Usp2) splicing was continuously activated; ubiquitin specific peptidase 36 (Usp36) and ubiquitin specific peptidase 9x (Usp9x) splicing showed fluctuating changes (). Our data suggest that AS changes may represent a new ubiquitin-proteasome activation mechanism associated with skeletal muscle atrophy.

Figure 5

Alternative splicing continues to affect the ubiquitin–proteasome pathway during muscle atrophy caused by hindlimb unloading. (A) In the Venn diagram of the differential alternative splicing genes of the skipped exon form at 5 time points, 201 genes overlapped. Functional analysis of overlapping genes. On the basis of biological processes (B) and molecular functions (C), these genes are enriched. (D) Dynamic changes in alternative splicing of 6 ubiquitin proteasome pathway genes during muscle atrophy. IncLevel on the y-axis represents the inclusion level of exons.

Identification of differentially expressed RBPs during skeletal muscle atrophy

To explore the mechanism of alternative splicing changes during skeletal muscle atrophy, we analyzed the expression of RBPs, which are considered splicing regulators. In the sequencing results, we identified a total of 15,781 genes that corresponded with 2,269 identified rat RBPs (42), and 1,639 RBP genes were obtained (). We identified the differentially expressed RBPs at 5 time points (P<0.01) and calculated the proportion of their numbers in 1,639 genes. The results showed that ~20% of RBP expression changed during skeletal muscle atrophy (). We took the intersection of RBPs that were differentially expressed at 5 time points to obtain the following 6 RBP genes: exportin 4 (Xpo4), eukaryotic translation initiation factor 4E family member 2 (Eif4e2), prolyl 4-hydroxylase subunit alpha 1 (P4ha1), LRR binding FLII interacting protein 1 (Lrrfip1), zinc finger CCCH type containing 14 (Zc3h14), and EMG1 N1-specific pseudouridine methyltransferase (Emg1). Analysis of the dynamic changes of these 6 genes during atrophy demonstrated continuous change; only Zc3h14 showed an upregulated trend, and the rest were all downregulated. Of these genes, Xpo4 demonstrated a continuous downward trend (). It has been reported that alternative splicing plays an important role in skeletal muscle biology, which may be strongly associated with the role of certain RBP family proteins (Hnrnp, Rbfox, Celf, Mbnl, Pabpc, and Pabpn) (21,43). Our results confirmed that Hnrnp h1, and Mbnl2 were significantly downregulated during skeletal muscle atrophy caused by unloading, and Rbfox1 and Mbnl1 were significantly upregulated during skeletal muscle atrophy caused by unloading (). The interaction network between these RBPs is shown in . The findings indicated that RBPs have complex interactions and connections, making the molecular mechanism of muscle atrophy more complex.

Figure 6

Identification of differentially expressed RNA binding proteins (RBPs) during skeletal muscle atrophy. (A) Recognition of RBPs in sequencing results. Blue circle represents the total number of RBPs, and orange circle represents the total number of genes identified in the sequencing results. A total of 1639 RBPs overlapped. (B) Histogram represents the number of differentially expressed RBPs at 5 time points during atrophy, and the percentage value at the top of the histogram represents the ratio of differentially expressed RBPs to the total number of RBPs identified in the study. (C) Heat map of the RBPs that continuously changed at the 5 time points in panel (B). (D) Heat map of RBPs reported in previous research that play an important role in the physiology of skeletal muscle. Genes with Fragments Per Kilobase of exon model per Million mapped fragments value >10 are shown, and the 4 RBPs with significant changes in expression are in red. Expression level of each gene in the control group is defined as 1, and the expression level at other time points refers to a multiple relative to 1. (E) STRING database was used to analyze the interaction network between RBPs. Line thickness indicates the strength of data support.

Discussion

In the present study, we conducted an alternative splicing analysis of the process of rat skeletal muscle atrophy caused by hindlimb unloading. We found that many differential splicing events occurred at 12 h, 24 h, 36 h, 3 days, and 7 days, and their main locations were the CDS region of the gene. Most of these alternative splicing events occurred in the CDS region, resulting in the gene’s frameshift. We found that alternative splicing activated the insulin signaling, endocytosis, mitophagy-animal, and ubiquitin-proteasome pathways at 36 h and 7 days through functional analysis. The alternative splicing of several USP-related genes continues to be altered during muscle atrophy. In terms of molecular mechanism, we found that 6 RBP genes’ expressions continued to change during atrophy, and 4 common RBPs showed significantly differential expressions. Our results emphasize that alternative splicing is a new focus of muscle atrophy, and in the present study, we preliminarily explored the molecular mechanism of alternative splicing changes.

The splicing of mammalian precursor mRNA is the process of intron removal and exon interconnection and is an important link during gene expression. In a specific time and space (such as different developmental stages or pathological conditions), cells usually produce proteins with different or even opposite functions by alternative splicing to adapt to environmental changes (44). Although alternative splicing is strongly associated with skeletal muscle physiology (26,45), the role of alternative splicing during skeletal muscle atrophy caused by hindlimb unloading has not been studied at the whole-gene level. In the present study, we found that a large number of differential alternative splicing events (genes) occurred during skeletal muscle atrophy caused by hindlimb unloading, which expanded our understanding of the process of skeletal muscle atrophy. It has been found that splicing disorders of some genes may cause skeletal muscle lesions, leading to muscle atrophy (46,47). Antisense nucleic acid technology is a potential therapy to correct and splice neuromuscular diseases. At present, antisense nucleic acid drugs for spinal muscular atrophy and Duchenne muscular dystrophy have been approved by the American Food and Drug Administration (48). The present study provides several new alternative splicing targets for the development of anti-atrophic drugs.

In the present study, we identified the differential alternative splicing events (genes) during muscle atrophy and analyzed their location on the gene and their impact on the open reading frame. The findings indicated that 59% of these differential alternative splicing events occurred in the CDS region and caused protein frameshifts. In addition to changing the reading frame, frameshifted transcripts were also one of the conditions for the activation of nonsense-mediated mRNA decay (49), and the transcripts that activated nonsense-mediated mRNA decay would eventually be degraded. Several recent studies have found that skeletal muscle atrophy was regulated by miRNAs, such as miR-29b (50), miR-125b-5p (18), and miR-351 (51). Although our results show that most AS occurs in the CDS region, a considerable portion still occurs in the UTR area. Alternative splicing events at 3'UTR will increase or decrease the binding site of miRNA.

Moreover, alternative splicing that occurs in 3'UTR can reduce the 3'UTR length to avoid miRNA or RBP-mediated regulation, it is also called alternative cleavage and polyadenylation (APA). APA also greatly contributes to the proteome’s selectivity and determines the ultimate fate of mRNA (52). Therefore, our results suggest that skeletal muscle atrophy is an extremely complex pathological process, and there are synergistic effects of multiple mechanisms.

The activation of the ubiquitin–proteasome pathway is one of the characteristics of skeletal muscle atrophy. In addition to the 5 classic molecules MAFbx, MuRF1, Traf6, UBR5, TRIM62, many other E3 ubiquitin ligases have been shown to regulate skeletal muscle atrophy (53-56). In the present study, 6 genes were found to be involved in the ubiquitin-proteasome pathway, including Usp15 and Usp2, and were found to have splicing differences. These splicing differences are widespread in muscle atrophy and provide a new perspective for the activation of ubiquitin-proteasome. Seven Usp2 gene alternative splicing products have been reported, 5 of which can encode proteins (57). The functions of the 3 main Usp2 alternative splicing products, Usp2-201, Usp2-202, and Usp2-204, have been extensively studied (51). Usp2-201 is involved in recycling epidermal growth factor receptor and plays a vital role in the cell cycle, Usp2-202 can activate the apoptosis signaling pathway and participate in apoptosis, and Usp2-204 can reduce the antiviral response of cells (58). Kotani et al. confirmed that the changes in inclusion (activation) or skipping (inhibition) of exon 7 of Usp15 could affect its substrate specificity and interactome (59). In our results, Usp15 splicing was inhibited during muscle atrophy, and Usp2 splicing was activated, implying that these changes may have a broad impact on substrate selection and downstream signaling pathways in the ubiquitin response.

RBPs are direct regulators of alternative splicing. The present study’s findings indicated that approximately 20% of RBPs were differentially expressed during skeletal muscle atrophy caused by hindlimb unloading, which could explain why a large number of alternative splicing events occur during atrophy caused by hindlimb unloading. We found that 6 RBP genes continued to change during muscle atrophy. Zc3h14 is generally upregulated during muscle atrophy and is involved in the formation of the poly(A) tail of mRNA, which affects the stability of mRNA (60); Xpo4 continues to be downregulated during muscle atrophy and encodes a nuclear transport protein that participates in the nuclear output of eIF5A and Smad3 (61). Although there is currently no evidence that Xpo4 are involved in alternative splicing regulation, they may still play an important role in muscle atrophy.

Multiple families of RBPs have been reported to participate in the physiological process of skeletal muscle (43), they may contain key alternative splicing regulators that we focus on. Singh et al. found that Rbfox1/2 is necessary to maintain skeletal muscle quality and proteostasis. Rbfox1/2 knockout in skeletal muscle has been found to cause changes in alternative splicing of many genes (32). Our results showed that the expression of Rbfox1 was disordered, suggesting that it may play an important role in skeletal muscle atrophy by regulating alternative splicing. Mbnl1 and Mbnl2 are also RBPs with large changes in expression during muscle atrophy. It has been reported that Mbnl1/2 is a skeletal muscle-specific RBP and participates in the tissue-specific alternative splicing regulation of skeletal muscle. Mbnl expression’s dysregulation will affect the severity of hereditary muscle atrophy, such as myotonic dystrophy (62). Interestingly, Rbfox1 and Mbnl1 have been reported to bind the same cis-acting element competitively and coordinately regulate alternative splicing (63), suggesting that the differential alternative splicing event at each time point may be the result of the interaction of many differentially expressed RBPs at that time point.

Splicing regulation depends on the interaction of cis-acting elements and trans-acting factors. Although we found several potentially differentially expressed RBPs, it is unclear why those RBPs are differentially expressed in atrophic skeletal muscle. Therefore, it is important to determine the complex regulatory relationship between them and the differential alternative splicing events in future studies so that the molecular mechanism of muscle atrophy can be further elucidated. This can also provide a scientific basis for the prevention and treatment of muscle atrophy.

The article’s supplementary files as