Eight-week-old male wild-type C57BL/6 mice were purchased from Guangdong Yao Kang Biotechnology for subsequent verification experiments. The total number of animals used for the experiments is 75. In which, 48 were used for hematoxylin and eosin (H&E) staining and immunofluorescence for embryonic MyHC (eMyHC) in controls and injured muscles covering 12 h, 24 h, 36 h, 48 h, 60 h, 72 h and 84 h post injury. H&E staining and immunofluorescence were both in triplicates for each time point. Additionally, 3 mice were just used for immunofluorescence for eMyHC in 96 h injured muscle. 24 mice (controls, 12 h, 24 h, 84 h post injury) were used for RNA-seq and quantitative polymerase chain reaction (qPCR; the total RNA samples of three replicates for each time point were mixed into a pool and each time point is in duplicates). Right side muscles were administered with 100 μL normal saline or 100 μL of 10 μM CTX (Sigma-Aldrich) directly into the middle of anterior tibial muscle (TA). Nitrogen inhalation was used for mice euthanasia. TA muscles were rapidly dissected from injured legs and healthy muscle was used as negative control [2].

Total RNA of TA muscles was extracted by using TRIzol reagent (Life Technologies) according to the manufacturer’s instructions. The RNA integrity and concentration were assessed by NanoDrop ND-1000 (Thermo Fisher Scientific). All the RNA samples were extracted according to the experimental requirements of the Illumina sequencing platform (RNA integrity numbers [RIN]≥7.0 and 28S/18S≥1.0). The quality assessment report of RNA samples was showed in Supplement 1. RNA libraries were constructed and RNA sequencing was performed by Beijing Genomics Institute (BGI) using an Illumina Genome Analyzer. Sequencing reads were generated using Illumina’s Digital Gene Expression Reads Profiling Kit according to the manufacturer’s instructions. These reads were aligned and identified, depended on two reference genomes, Ensemble mouse (000001635.20) and the UniGene database. Clean reads were obtained after removing low quality sequences (there were more than 50% bases with quality lower than 20 in one sequence), reads with more than 5% N bases (bases unknown) and reads containing adaptor sequences, and then they were aligned to the reference database and annotated.

High-throughput sequencing was conducted to produce transcriptome profiles of skeletal muscle at the early stages of muscle regeneration in mice. After filtering adaptor reads, empty reads and low quality reads from the raw data of 8 muscle samples, 5.9±0.06×107 clean reads were remained (Supplement 2). In addition, the base along clean reads had a balanced composition and high quality (Supplements 3A, 3B). Through alignment of clean reads, 82.5±0.81% of clean reads were aligned to reference genome and 74.4±4.7% of clean reads were mapped to UniGene database (Supplement 4). In addition, the reads of 8 muscle samples were evenly distributed on reference genes (Supplement 3C).

Reads that could be uniquely mapped to a gene were used to calculate the expression level. Using RSEM tool, fragments Per Kilobase of exon region per Million mapped reads (FPKM) was calculated, whose value can be directly used as the normalized gene expression level. The fold changes (log2 ratio) were estimated according to the normalized gene expression level in each sample. We employed the absolute value of log2 ratio≥1 and probability≥0.8 as the threshold to screen the significantly differentially expressed gene (DEG).

To determine the main biological functions, DEGs were annotated by terms of the Gene Ontology (GO) database (http://www.geneontology.org/). GO enrichment analysis of functional significance terms in the GO database was applied using hyper geometric test to find significantly enriched GO terms in DEGs comparing to the genome background. The calculated p-value of GO terms goes through Bonferroni Correction, taking a false discovery rate (FDR)<0.05 as a threshold to screen significantly enriched GO terms in DEGs. In addition, GO functional analysis of candidate gene was carried out by DAVID analysis (http://david.abcc.ncifcrf.gov/). The molecular interaction network of candidate genes was constructed through STRING (http://string-db.org/) and Cytoscape 3.1.0 (http://www.cytoscape.org/). High confidence (0.7) was used in STRNG analysis.

Genes usually interact with each other to take part in certain biological functions. Pathway functional analysis of up-regulated DEGs was carried out by DAVID analysis (http://david.abcc.ncifcrf.gov/). Up-regulated DEGs of our comparison groups were classified into many significant enrichment pathways. Significantly enriched pathway was defined by p value<0.05.

11 genes were chosen to validate the accuracy of RNA sequencing. The remaining RNA from sequencing was used. cDNA was generated from 1 μg of total RNA by reverse transcription system (Promega) according to the manufacturer’s instructions. The qPCR reaction was carried out in the LightCycler 480 II system (Roche), using 2×SYBR Green Realtime PCR Master Mix (Takara) according to the manufacturer’s instructions. The analytic method of 2^−ΔΔCt was used. The Spearman’s correlation coefficient was calculated for the gene to assess the consistency between qPCR and RNA sequencing. Supplement 5 lists the gene-specific primers for gene expression.

TA was dissected from injured muscle, fixed in 4% paraformaldehyde (PFA) for 24 h at 4°C, then embedded with paraffin. Paraffin-embedded muscle was sectioned at 4 μm. After deparaffinating, muscle section was stained by H&E. Cell nucleus was stained by hematoxylin in 1 min, washed by 95% ethyl alcohol. Cytoplasm was stain by eosin in 20 seconds, washed by 95% ethyl alcohol, finally washed in PBS and mounted by neutral balsam. H&E staining images of muscle sections were taken by a fluorescence inverse microscope (ZEISS).

For ICC, OCT-embedded TA muscles were sectioned at 6 μm, fixed for 10 minutes in acetone at −20°C, permeabilized with 0.5% of Triton X-100 in PBS for 15 minutes followed by being blocked with Hydrogen Peroxide Block for 15 minutes, washed by PBS, blocked with Rodent Block (ab127055; Abcam) for 30 minutes, then incubated with primary antibodies against eMyHC (F1.652-b; DSHB) and LaminA2 (ab11576; Abcam) for 1 hour at room temperature. After removing primary antibodies and washing, muscle sections were incubated by anti-Mouse IgG secondary antibodies (#4409, CST) for eMyHC and anti-Rat IgG secondary antibodies (#4418S, CST) for LaminA2 in 10 min at room temperature. After being washed by PBS, sections were stained by DAPI for 2 minutes. Lastly the muscle sections were mounted by anti-fade Mounting Medium. Images were taken by confocal microscopy (Leica).

For screening the DEGs of each comparison group, we employed the absolute value of log2 ratio≥1 and probability≥0.8 as the threshold. FDR<0.05 was applied to screen significantly enriched GO terms in DEGs. In addition, p value<0.05 was used to define significantly enriched pathways. The Spearman’s correlation coefficient was calculated for the gene to assess the consistency between qPCR and RNA sequencing. Significant correlation was accepted when p<0.01.

To characterize the morphology of CTX-injured muscles, we performed H&E staining to characterize the morphologies of negative controls (NC) and injured muscles covering 12 h, 24 h, 36 h, 48 h, 60 h, 72 h and 84 h post injury (Figure 1A). At CTX12, the number of nucleus was obviously more than that in NC sample, which suggested that inflammatory cells had already migrated and proliferated. Subsequently, myofibers underwent necrosis and edematous (Figure 1B). Notably, the H&E results at CTX36, CTX48, and CTX60 closely resembled those observed at CTX24, with persistent necrosis and inflammatory infiltration. 8 muscle samples in total were grouped into four categories, of which NC sample was far away from injured muscle samples, the CTX12 sample was closely next to CTX24 sample and far away from CTX84 sample (Figure 1C).

To uncover the genes involved in early muscle regeneration, we totally obtained 4,774 DEGs with the absolute value of log2 fold change≥1 and probability≥0.8 as the threshold. In Figure 1D, totally 1,771 DEGs were up-regulated and 3,003 DEGs were down-regulated in three comparison groups. Additionally, in NC-CTX12 comparison group, 2,395 up-regulated and 1187 down-regulated genes were identified. In CTX12-CTX24 comparison group, we found 268 up-regulated and 309 down-regulated genes while 448 up-regulated and 394 down-regulated genes in CTX24-CTX84 comparison group (Figure 1E). These results indicated that the number of DEGs in NC-CTX12 is far more than that in CTX12-CTX24 or CTX24-CTX84, suggesting that the injured-activated genes at CTX12 are the most robust during early muscle regeneration.

To identify signaling pathways involved in early muscle regeneration, DEGs were mapped to the KEGG database. Significantly enriched functional pathways of up-regulated DEGs in three comparison groups were presented in Table 1. In NC-CTX12 comparison group, MAPK, Notch and chemokine signaling pathways were dramatically activated, which have previously been identified to regulate satellite cell proliferation [10,11]. DNA replication and cell cycle related genes were enriched in CTX12-CTX24 comparison group. Focal adhesion and extracellular matrix (ECM)-receptor interaction pathways were enriched in CTX24-CTX84 comparison group. More importantly, p53 signaling pathway was always significantly enriched in three comparison groups, which indicated that p53 pathway was involved in entire early process of muscle regeneration process and plays important role (Figure 2G).

To further study the biological function of DEGs during early muscle regeneration, the DEGs were classified by GO analysis. Bar graph shows that DEGs were mostly enriched in IS, cell proliferation and cell activation in NC-CTX12 group (Figure 2A). In CTX12-CTX24 group, DEGs were enriched in cell cycle, cell proliferation and IS processes (Figure 2B). In CTX24-CTX84 group, DEGs were enriched in IS, cell differentiation, muscle system (MS) and cell adhesion (Figure 2C). Interestingly, DEGs involved in IS and hormone response (HR) were shared in three comparison groups. Further, we totally found 516 immune related and 177 HR related up-regulated DEGs during early muscle regeneration. These implied that IS and HR were activated for regulating early muscle regeneration. Furthermore, up- and down-regulated DEGs of MS, IS and HR were classified. Statistical analysis showed that the number of up-regulated IS and HR related DEGs was more than that of down-regulated DEGs while the number of up-regulated MS related DEGs was less than that of down-regulated DEGs in NC-CTX12 group. However, the number of up-regulated IS and HR related DEGs was less than that of down-regulated DEGs in CTX12-CTX24 and CTX24-CTX84 groups while the number of up-regulated MS related DEGs was more than that of down-regulated DEGs in CTX24-CTX84 group (Figure 2D). These results indicated that the difference of DEG number between up- and down-regulated IS was similar with that of HR, but contrary to MS (Figures 2E, 2F). As IS is known to trigger muscle regeneration, we assumed that hormone regulation may also play an important role in early muscle regeneration.

To comprehensively explore how IS regulates satellite cell growth, 6 key immune factors regulating satellite cell proliferation were selected and these gene expression patterns were presented in heat maps. The heat maps illustrated that they displayed highly similar expression patterns (Figure 3A). Genes with similar expression patterns usually have similar biological function. Therefore, according to aforementioned expression patterns, 326 candidate genes were obtained. Then enrichment analysis was performed to inspect whether these genes shared common features. In Figure 3B, these candidate genes were mostly associated with IS, cell proliferation, activation and migration. The interaction network of these genes was predicted using STRING website and illustrated by Cytoscape software. It showed that 41 immune factors were identified and closely interconnected with genes involved in cell proliferation, activation and migration (Figure 3C), which implied that these immune factors potentially regulated satellite cell activation, migration and proliferation. Further, we calculated the fold change of 41 candidate immune factors expression level in injured muscle compared to NC (Supplement 6) and found that fold changes of S100a9, Csf3r, Cxcl3, Ppbp, Ccl3 and Il1rn were over 30 times, indicating that these immune factors may play more important role in satellite cell regulation during muscle regeneration.

To characterize the formation of newly formed myofibers during early muscle regeneration, eMyHC, a marker of newly formed myofibers, was detected by ICC. The results showed that eMyHC was extensively expressed in skeletal muscle at CTX84 (Figure 4A), suggesting that genes involved in myofiber formation were activated at CTX84. To further investigate candidate genes governing myofiber formation, key myogenic factors were chosen and their mRNA expression patterns during muscle early regeneration were presented by heat maps. It showed that these genes were significantly up-regulated at CTX84 (Figure 4B). Further, 320 candidate genes that shared similar expression patterns with myogenic factors were obtained. GO analysis indicated that most of these genes were associated with cell adhesion and muscle development (Figure 4C). STRING prediction result showed that these interacting proteins were classified into three kinds of functional proteins (myofibrillar proteins, myogenic proteins and cell adhesion proteins) (Figure 4D). Obviously, these cell adhesion proteins influence the formation of muscle fibers by regulating myogenic proteins during muscle regeneration process.

To validate RNA sequencing data, relative mRNA expression of 11 DEGs were detected by qPCR (Figure 5). Specifically, Myod1, Myf5, Pax7 and MyoG are crucial regulators in myogenesis and muscle regeneration [6]. TNFa, IL1β, IL6, Ccl2, Ccl3 and Ccl4 play important roles in regulation of satellite cell proliferation [7]. In addition, eMyHC involved in regeneration of myofibers upon injury is encoded by MyH3 [12]. The correlation between qPCR results and RNA sequencing data of DEGs was carried out by Spearman’s correlation coefficient analysis and the result showed that the correlation coefficient (R) was greater than 0.8 in all the genes, indicating that expression profiles of these genes have high consistency between two methods.

In our study, GO analysis of significantly up-regulated DEGs during early muscle regeneration indicated that several pathways known to be closely involved in regulation of satellite cell activity were enriched (Table 1). Specifically, chemokine and notch signaling have been reported to promote satellite cell activation, migration and proliferation [15–17]. Previous reports have shown that p53 acts as suppressor to regulate cell proliferation and apoptosis [18]. Moreover, p53 signaling pathway is activated during myogenesis and p53 interference prevents C2C12 myoblast differentiation [19]. Furthermore, recent findings have shown that mechano-signaling via Piezo1 prevents the activation of p53 and its associated senescence program in muscle stem cells, thereby promoting regenerative capacity following injury [20]. In our data, p53 pathway was enriched throughout the stages of early muscle regeneration, as evidenced by the sustained upregulation of key target genes (Figure 2G), indicating that it is very important during the entire early muscle regeneration.

Hormones, a class of secreted factors including adrenal cortical hormone, sex hormone, thymin, oxytocin and so on, plays important roles in development and physiological metabolism [21]. Recently, several studies have focused on the effects of hormones on muscle regeneration. Systemic administration of oxytocin promotes satellite cell proliferation for muscle regeneration, while satellite cell proliferation is impaired by inhibiting oxytocin signaling [22]. In addition, inactivation of thyroid hormone (T3) by type 3-deiodinase (D3) improves muscle regeneration through surviving satellite cell proliferation [23]. Decreased androgen impaired muscle regeneration [24]. Moreover, Liao et al [25] reported that estrogen promotes the polarization of macrophages toward an anti-inflammatory M2 phenotype and modulates T cell responses in injured muscle, thereby facilitating repair. Similarly, insulin and amino acids cooperatively activate the mTORC1 pathway in skeletal muscle, supporting regeneration through converging hormonal and nutritional signals [26]. Based on our transcriptome analysis, the number of DEGs between up- and down-regulated IS was similar with that of HR (Figure 2D). It was well known that IS is a key process in early muscle regeneration. Therefore, it allowed us to think that hormones may play an analogical role in early muscle regeneration, which could be a novel strategy for the treatment of muscle disease.

After muscle injury, the key inflammatory cells are rapidly activated and secrete immune factors to govern the microenvironment where satellite cell locates [27]. It has been determined that several important immune factors are able to regulate satellite cell activity for muscle regeneration. In particular, chemokines, including Ccl2, Ccl3 and Ccl4, directly promote myoblast proliferation by regulating the phosphorylation of ERK1/2 [16]. In parallel, the granulocyte colony-stimulating factor receptor (Csf3r) facilitates muscle repair by mediating neutrophil recruitment [28], while the interleukin-1 receptor antagonist (Il1rn) fine-tunes the inflammatory response to prevent excessive tissue damage [29]. IL1β and TNFα directly activate NF-κb to enhance the proliferation of satellite cells through promoting expression of cyclinD [30]. In addition, knocking out IL6 impairs satellite cell proliferation by decreasing STAT3 phosphorylation and cyclinD1 expression [31]. During early muscle regeneration, these immune factors had similar expression patterns (Figure 3A), which suggested that immune factors that have similar patterns may regulate satellite cell growth during early muscle regeneration.

To explore more candidate immune genes regulating satellite cell growth, we screened 58 immune genes that shared similar expression patterns with immune factors mentioned above (Figure 3B). Interaction network presented those 41 immune factors showed strong relationship with cell migration, activation and proliferation factors (Figure 3C). Further, we focused on S100a9, Csf3r, Cxcl3, Ppbp, Ccl3, Il1rn that may be crucial to migration, activation and proliferation of satellite cells because of their high expression level during early muscle regeneration. S100a9 and Ccl3 have been identified to regulate C2C12 cell proliferation [16,32]. In addition, Csf3r mediates hematopoietic progenitor mobilization in mice [33]. Cxcl3 regulates smooth muscle cell migration [34], and Ppbp facilitates intravascular leukocyte migration directly [35]. Therefore, Csf3r, Cxcl3, Ppbp may be also involved in regulation of satellite cell migration. Interleukin 1 receptor antagonist (Il1rn) inhibits the activation of NF-κb which is a crucial factor in satellite cell proliferation, suggesting that Il1rn could regulate satellite cell proliferation through NF-κb. Together with our study, these candidate immune genes may be involved in satellite cell activity during early muscle regeneration.

In our study, newly formed myofibers were widely formed at CTX84, meanwhile, the mRNA level of muscular regulators was significantly up-regulated (Figures 4A, 4B). Based on the expression patterns of muscular regulators in Figure 4B, we screened 320 genes, of which 27 genes were cell adhesion factors (Figure 4C), meanwhile ECM-receptor interaction and focal adhesion signaling were activated at CTX84 (Table 1), implying that 27 cell adhesion factors were potentially involved in myofiber formation. Moreover, molecular network showed that 16 cell adhesion proteins were related to myogenic regulators and myofibrillar proteins (Figure 4D). Previous report has found that cadherin regulates phosphorylation of β-catenin to influence myoblast differentiation [36], which implied that cdh2 and cdh15 may regulate myoblast differentiation during muscle regeneration. In addition, ablation of fibroblast in skeletal muscle promotes satellite cell premature differentiation, leading to impair muscle regeneration [37]. Collagen is mostly secreted by fibroblasts and acts as framework for guiding satellite cells to facilitate muscle regeneration [9,38,39]. Further, knocking out Col6a1 in mice attenuates satellite cell proliferation, resulting in muscle regeneration deficiency [5]. These findings suggested that collagen plays crucial roles in muscle regeneration. In our molecular network, 9 collagen proteins were found, which may play important roles in early muscle regeneration. In addition, Thbs2, a member of ECM proteins, can mediate collagen fibrillogenesis, which hinted that thbs2 could indirectly influence muscle regeneration through regulating collagen formation. Our results indicated that candidate cell adhesion factors may be involved in muscle regeneration (Figure 6).