First, we built an early state skeletal muscle injury animal model in 0–48 h using rats (see METHODS section). Then, we examined the injured rats in groups at seven time points after SMI (Figure 1B). To confirm the solid of animal model and the processes of muscle tissue injury to repair, we performed a preliminary histological analysis with hematoxylin and eosin (H&E) staining. H&E staining revealed the changes in muscle morphology after injury (Figure 1A). In comparison with the normal control group, the transverse sections of skeletal muscles exhibited breakage and the destruction of large numbers of muscle fibers at 4 h after injury. At 8 h after injury, interstitial blood vessels were congested with the infiltration of large numbers of inflammatory cells. At 12 h after injury, muscle fibers exhibited obvious hyperemia and edema, which reflected a strong inflammatory response at this time point. At 16 and 20 h after injury, the inflammation was reduced and the amount of fibrous tissue in the injured area gradually increased. At 24 and 48 h after injury, the numbers of fibroblasts increased, and muscle fibers were partially rebuilt. This was observed in longitudinal sections of muscle fibers at 4 and 48 h after injury (Supplementary Figure 1).

To objectively describe and confirm the processes of pathological changes in the early period after SMI based on transcription profiles, we examined genes with obvious changes in expression, Gene Ontology (GO) terms, and pathways during the early phases of the post-injury response. In this study, we obtained raw next-generation sequencing data from 60 samples (in total 20,662 genes) (Supplementary Table 1). Then we subjected the standardized data (Supplementary Figure 2) to principal components analysis (PCA). After removing abnormal samples, the results of PCA indicated that the samples from each injury group were clearly discriminated from the normal samples and also that the groups at 4, 8, 24, and 48 h after injury were clearly distinct from one another. The groups at 12, 16, and 20 h after injury were not clearly differentiated (Supplementary Figure 3). These divisions were also reflected in the clustering correlation heat map (Supplementary Figure 4). Therefore, after removing genes with low expression, 15,901 common genes expressed in 53 samples (Table 1) were subjected to weighted gene co-expression network analysis (WGCNA) (Supplementary Table 2) (Zhang and Horvath, 2005). A soft power of 22 was selected to construct a scale-free network (Figure 1C). We then used dynamic tree cutting to divide the gene co-expression modules according to the common gene expression file. The clearest gene co-expression module distribution obtained at early post-SMI period using the Dynamic Merge tool included 24 modules comprising 23 meaningful modules and one undefined gray module (Figure 1D) (Zhang et al., 2019b).

To understand the correlations between each co-expression module and specific time points after injury as well as the potential biological significance of these modules, we used Spearman correlation analysis to calculate R values and p-values. Correlations with time and the biological significance of the 23 meaningful modules were visualized in a module-trait relationship heat map (Figure 1E). To analyze the expression of related genes in each meaningful module, we selected modules associated with positive time-point correlations (R > 0) (Figure 1E), setting p < 0.05 (Jin et al., 2019) as the threshold value. For each time point, we only considered the module with the highest R value.

Six co-expression modules (represented by different colors in Figure 1E, in this study we named the modules using the name of corresponding time point) were selected based on the results of WGCNA analysis (Table 1; gene annotation information for each selected module listed in Supplementary Table 2). However, we did not find an eligible module in 16 h group (Table 1 lists the correspondence between module color and time. The modules corresponding to time point in the following text are directly named by time). Finally, the most relevant module in each group, referred to as a highly correlated (HCr) module, was selected for gene expression analysis. Next, to confirm that the module division according to WGCNA and the selection of HCr modules based on gene co-expression profiles as outlined above (Table 1) were reasonable, we explored the relevance of gene co-expression in the HCr modules using heat maps of all genes that clustered in a given module. The results indicated that the modules identified by this analysis largely corresponded to the variation in gene expression at specific time points after SMI (Supplementary Figure 5). In addition, expression levels (corresponding to the feature vectors) (Zhang and Horvath, 2005) in the module of eigengene (represented by sparse loads of the most important feature vectors and variance directions (variable genes with maximum variable percentage) corresponding to the major contributions of the expression profile) (Shen et al., 2006) were plotted in histograms based on transcriptome analysis results for each HCr module in all time-point groups (Supplementary Figure 5). In detail, in the “4,” “8,” “12,” and “48” modules, eigengene expression levels were highest at 4, 8, 12, and 48 h, respectively, among all the time points after injury examined in this study. In the “20” and “24 modules, eigengene expression levels were relatively high at 20 and 24 h, respectively (Supplementary Figure 5). Based on eigengene expression in each selected HCr module, we could use these modules for further analyses.

After selecting the most relevant meaningful module in each group, we analyzed gene expression profiles to reflect the transcriptional kinetic characteristics of each module. We performed functional analysis to determine the transcriptional dynamics changes during the early post-SMI period.

We focused on biological processes among the GO terms that were enriched in the selected modules and selected the GO terms of interest based on their associated p-values (p < 0.01). Due to the unequal numbers of genes in each selected module, the numbers of biologically meaningful GO terms chosen differed among modules (Table 1). The screening process of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways for each selected module was similar to that described for GO terms. According to the above principles, we focused on several phases of pathological changes in the early period after SMI, Based on this, the GO terms and KEGG pathways screened are shown in Supplementary Figure 6.

As outlined above, the entire course from SMI to repair can be roughly divided into five interconnected biological stages: response to stimuli and stress, inflammation, immune process, cellular process, and repair process (Figure 2A). In the stage associated with stimuli and stress, most GO terms were enriched in the “4” and “8” modules. As shown in Figure 2A, almost all GO terms enriched in the “12” module were associated with inflammation and the immune process. Especially the 48 h module, most of the enriched GO terms were related to the repair process.

During the entire course of the SMI repair process, based on the enrichment of GO terms in all HCr modules, each stage could be further divided into some sub-classes (Figure 2A). In the stage associated with stimuli and stress, the response to external stimuli and the stress began immediately after injury. The oxidative stress responses also began in the early period of the SMI repair process. We divided the inflammation stage into sub-stages associated with the release of inflammatory mediators, inflammatory processes, and the inflammatory cell response.

In particular, among the sub-classes pertaining to inflammatory cell response, the migration of various inflammatory cells (neutrophil, macrophages, lymphocyte, and mononuclear) exhibited a time-varying pattern. From the number of enriched GO terms, we inferred that inflammation predominated from 4 to 12 h after injury. Similarly, in the immune process stage, most immune-related GO terms were enriched in the “12” module. By contrast, autophagy was only enriched in the “24” module (Figure 2A). In addition, the GO terms of DNA repair, recombination, and regulation of DNA-dependent DNA replication were enriched in the “48” module, suggesting the involvement of DNA-related intracellular mechanisms in the course of tissue repair (Supplementary Figure 6A).

Taken together, the findings outlined above demonstrated the temporal variation in GO terms enriched in each HCr module during the entire course of the SMI repair process. The changes in gene expression profiles revealed the basic biological processes and pathological changes occurring during the SMI repair process.

Similar to GO term enrichment analysis, the entire course of the SMI repair process could be roughly divided into stages associated with damage, the repair process, and signal regulation based on KEGG pathway analysis (Figure 2B). We also classified the KEGG pathways into 2–3 sub-classes for each stage of the SMI repair process (Figure 2B). The distribution of KEGG enriched pathways varied in a time-dependent manner. For example, the IL-17 signaling pathway, RIG-I-like receptor signaling pathway, and Toll-like receptor signaling pathway, which are related to the production of inflammatory factors, were enriched in the “4” and “8” modules, whereas pathways related to inflammation and the immune response, such as natural killer cell-mediated cytotoxicity, were enriched in the “12” module (Figure 2B). Pathways related to DNA replication and cell cycle in the repair stage were enriched in the “48” module (Figure 2B). These observations were consistent with the time-dependent variation exhibited by GO terms and the mechanisms underlying the muscle damage-repair process. Taken together, GO and KEGG enrichment analyses for all HCr modules indicated the temporal pathological variation and regulatory processes occurring during the early phase of SMI repair in rats.

In this turn, temporal variation in the active biological processes and pathways after injury were precisely regulated by some aspects of regulation-related KEGG pathways, especially during the initial period of the SMI repair process (i.e., in the “4” module; Figure 2B). The enriched regulatory KEGG pathways included those related to cascade reactions (such as the NOD-like receptor and MAPK signaling pathways); transcription factors (such as the FoxO signaling pathway); pathway switching, including protein phosphorylation regulation and histone modification (such as the Jak-STAT signaling pathway and ubiquitin-mediated proteolysis); and second messenger regulation (such as the cAMP signaling pathway) (Yang et al., 2007; Sharma et al., 2012; Icli et al., 2019).

In summary, throughout the entire SMI repair process, temporal variation in the distribution of KEGG pathways and biological processes (as indicated by GO term enrichment) in each HCr module objectively reflected the transcriptional dynamics as observed in gene expression profiles and the biological significance of the HCr modules (Figure 2).

To further confirm the HCr module selected previously to represent each period, and to study the expression of representative functional genes in each HCr module, we screened the top 30 high-connectivity genes in each HCr module based on the results outlined above. First, we entered all co-expressed genes in each module (Table 1) into the web-based analysis program Cytoscape (Kohl et al., 2011) to screen for the top 30 high-connectivity genes in homologous HCr modules. Due to the differences in numbers of co-expressed genes in each HCr module (Table 1), the screening process yielded a total of 161 genes distributed in the HCr modules (Figure 3A).

We performed detailed enrichment analyses of GO terms representing biological processes and KEGG pathways for these 161 high-connectivity co-expressed genes (Supplementary Figure 7). The bubble diagram of GO terms in Supplementary Figure 7A lists the details of the top three GO terms in each HCr module discussed above (in total, 50 GO terms representing biological processes were screened). As shown in Supplementary Figure 7B, seven KEGG pathways related to the immune response, proteasome, and repair-related process were enriched.

Next, to better show the significance of these genes in each HCr module, a heat map of the top high-connectivity co-expressed genes was produced, which clearly revealed some genes were more highly expressed at specific time points after injury, indicating the relevance of the corresponding HCr module to each time point (Figure 3B). The expression time-series graphs to the right of the heat map showed the relative mean expression of genes in each HCr module, reflecting the significance of these high-connectivity co-expressed genes. Next, GO term enrichment analysis was performed to determine the functional representativeness, of the top high-connectivity co-expressed genes in each HCr module. Among them, we screened the top three representative GO terms in each HCr module to determine the functions of the high-connectivity co-expressed gene cluster and to confirm the transcriptional dynamics in the SMI repair process (Figure 3B). The dominant functions of the co-expressed genes (eigengene in part) in each HCr module matched those expected from the classical biological process of SMI repair illustrated in Figure 2A. Moreover, the HCr modules screened from the results of WGCNA and the high-connectivity co-expressed genes in each module represent the early period of the SMI repair process and reflect the transcriptional dynamics of the gene expression profiles.

Having analyzed the HCr modules identified by WGCNA and expanded on the transcriptional dynamics of the gene expression profiles in all time periods covering the early stages of the SMI repair process, we focused on the time-dependent variation in gene expression. Based on the data analysis workflow (Supplementary Figure 2), we identified time-dependent genes using the maSigPro (microarray Significant Profiles) R package (version 2.0.2) (Nueda et al., 2014).

As described above, raw clean data on the experiment-wide gene expression profiles were obtained from 53 samples and 20,662 genes (Supplementary Table 1). The results of maSigPro analysis indicated that during the early period after SMI, 1512 genes exhibited marked temporal changes in expression at the 95% confidence level. Using hierarchical clustering and maSigPro default parameters, these genes exhibiting time-dependent changes in expression were divided into six clusters with different expression profiles during the experiment (Figure 4A and Supplementary Table 3) (Nueda et al., 2014). The raw profiles for each group in each cluster are presented in Supplementary Figure 8A, and the median profile of gene expression and profile differences between the time-dependent experimental and control groups are presented in Figure 4A. In comparison to the control group, the experimental groups exhibited marked downregulation of genes in Cluster 2 and upregulation of genes in Cluster 5. Thereafter, the expression levels of these genes mostly increased in Cluster 2 and decreased in Cluster 5 throughout the whole period after injury, coming close to the levels in the control group. The genes in Cluster 1, Cluster 3, and Cluster 6 were upregulated from 4 to 48 h after injury, whereas the genes in Cluster 4 exhibited the opposite trend (Figure 4A). Due to the different numbers of genes in each cluster (divided according to maSigPro), we chose the high-connectivity genes in the top 10% from the six clusters using a protein–protein interaction (PPI) network visualized with Cytoscape software (version 3.7.2), and a total of 152 genes were selected (Supplementary Figure 8B).

Next, to study the temporal variation in the expression of DEGs at different time points, we performed pattern analysis using STEM software (Ernst and Bar-Joseph, 2006), which has been widely used for in-depth analyses of gene-expression time-series clustering data.

Comparison of the gene expression profiles between the control group and the experimental groups revealed that there were 540 DEGs expressed in common in the seven experimental groups. The distribution of all DEGs in 7 experiment groups (|log₂FC|≥1) was shown in Supplementary Figure 9A.

These genes are common differential expression gene sets that at each stage of 0–48 h after skeletal muscle injury compared with the control group. According to the GO and KEGG enrichment results of this part of genes (Supplementary Figure 10), the biological functions of the common differential genes at 7 time points after injury enriched mainly included immune performance, inflammation, cell cycle and stress response. From this point of view, according to the literature (Crow et al., 2019), the set of 540 common differential genes can be regarded as systemic global prior differential genes (DE genes).

STEM software was used to identify dynamic gene expression clusters in the common 540 DEGs using default parameters (Supplementary Table 4). The results revealed six expression patterns in the DEGs (p < 0.001; Supplementary Figure 9B); however, no similar gene expression trends among these patterns were found (Figure 4B). The analysis also identified 439 common DEGs exhibiting dynamic expression at different time points after injury (Supplementary Table 5). Similar to the results of maSigPro analysis, we then chose high-connectivity genes in the top 10% representing the six expression patterns using a PPI network that was visualized with Cytoscape software, and a total of 41 genes were selected (Supplementary Figure 9C). Furthermore, 41 time-related highly linked genes screened by post-injury common differential genes (global prior differential genes) may play a role in global regulation.

Taken together, the maSigPro analysis based on raw RNA-seq data and STEM analysis based on common DEGs revealed time-dependent marker genes and genes exhibiting dynamic expression during this phase from two perspectives.

To comprehensively study the function of the time-dependent marker genes selected as outlined above using an in-depth analysis, we integrated the 152 genes screened by maSigPro and the 41 genes (global prior differential genes) selected by STEM. To obtain marker genes fully reflecting temporal variation during the SMI repair process, we extracted the union set from the total of genes and identified 188 genes for further in-depth analysis (Figure 4C).

This set of 188 genes included marker genes exhibiting time-dependent variation in expression and genes with dynamic expression during the early period after SMI. The circular heat map in Figure 4D indicates that throughout the whole period after SMI, nearly one-third of these genes were downregulated, whereas the remaining two-thirds were upregulated. This observation indicated that these 188 genes exhibited time-dependent variation in expression and could be used as time-dependent marker genes.

We then examined the expression of the 188 time-dependent marker genes using GO term enrichment (biological processes) and KEGG pathway analyses. These 188 genes were associated with 737 GO terms and 48 KEGG pathways. In the GO term enrichment analysis, we selected GO terms based on the p-value (p < 0.01) and identified 52 GO terms (Supplementary Figure 11A). The selected GO terms indicated that the 188 genes reflected the mechanisms underlying the early stages of the SMI repair process. Then, we identified 42 KEGG pathways, which are listed in Supplementary Figure 11B. The KEGG pathways associated with these 188 genes also reflected the precise regulation of multiple processes after injury, e.g., the TNF signaling pathway, IL-17 signaling pathway, osteoclast differentiation, cell cycle, and chemokine signaling pathway (Supplementary Figure 11B).

In summary, the results suggested that these genes not only reflected time-dependent changes but also included hub genes and core regulatory genes participating in the many complex processes involved in SMI repair.

Having illustrated the characteristics of the transcriptional dynamics and time-dependent gene expression variation, in this section, we further illustrate the possible connection between transcriptional dynamics that reflect biological processes and time-dependent marker genes to identify potential biomarkers that will allow us to elucidate the correlations between temporal changes in gene expression and the biological processes that occur in the early period after SMI. In addition, we verified our results based on morphological and immunohistochemical analyses.

First, we identified the genes in common between the 161 high-connectivity co-expressed genes (eigengenes) in the HCr modules and the 188 time-dependent marker genes. The results indicated that there were 18 genes common to the two gene sets (Figure 5A). Next, we analyzed the positions of these 18 genes in a network in which the set of 161 high-connectivity co-expressed genes and the set of 188 time-dependent marker genes were combined, as well as the functions of the 18 genes during the early period after SMI. In Figure 5A, the 18 genes are listed in a red-outlined box in the middle of the two gene networks. Genes in these two networks are connected by lines of different thickness. The functions of the 18 genes also corresponded to the main biological processes discussed above, and they may be hub genes in these complex processes after injury (Figure 5A and Table 2). In addition, we produced a heat map of the 18 genes to visualize changes in their expression levels during the early period after SMI. The expression levels of the bottom one thirds of genes from the “4” and “8” modules were initially upregulated, and then gradually downregulated. By contrast, the expression levels of the top two thirds of genes from the “48” module exhibited the opposite trends (Figure 5B).

These observations indicated that the transcriptional dynamics of the co-expressed genes in the HCr modules were closely related to the time-dependent marker genes through these 18 genes. In addition, based on the integration of transcriptional dynamics and the time-dependent marker genes, some of the hub genes may reflect temporal gene-expression changes associated with biological processes occurring in the early stage following SMI.

Next, to intend confirm the process of inflammation during the early period after SMI, we labeled neutrophils using an anti-MPO primary antibody. The results of immune-histochemical analysis suggested that, in comparison with the normal control group, with the development of inflammation, the degree of neutrophil infiltration and the neutrophil distribution varied in a characteristic manner during the early period (0–24 h) after SMI (Figure 5C).

Taken together, the results of the histo-morphological and immune-histochemical analyses reflected the characteristics of the SMI repair process.

Male 6–7-week-old Sprague Dawley rats were obtained from the Experimental Animal Center of Shanxi Medical University (Taiyuan, Shanxi, China) and housed in environmentally enriched ventilated cages under specific pathogen-free conditions at the Provincial Key Laboratory of Forensic Medicine, Shanxi Medical University (Jinzhong, Shanxi, China) under a 12h light/dark cycle with ad libitum access to water and food [RM1 (P); Regular Diet Services]. Sixty-four male Sprague Dawley rats were divided randomly into eight groups (seven experimental groups and one control group, n = 8 rats/group) once their body weight reached 240 g (±20 g) (Figure 1B). The animals in the experimental groups were anesthetized with 10% chloral hydrate (2.5 ml/kg body weight, intraperitoneal injection) and the hair on their right posterior limb was removed using a depilatory agent. A 500 g weight was dropped from a height of 50 cm onto the right hind limb with a free fall motion. The rats in the experimental groups were sacrificed at 4, 8, 12, 16, 20, 24, or 48 h after injury using a lethal dose of 10% chloral hydrate (4 ml/kg body weight, intraperitoneal injection). The animals in the control group were directly sacrificed using the same method as described above in the middle of the 0–48 h post-injury period defined in this study. A muscle sample of ~200 mg was dissected from the wound site of each rat and divided into two parts: (1) samples from the central area of damage (~100 mg) were snap-frozen with liquid nitrogen and cut into small squares for histological analysis after fixation with formalin solution; and (2) the remaining samples (~100 mg) were frozen immediately in liquid nitrogen and saved at −80°C until RNA-seq (50 mg) analysis after the fascia was removed.

Frozen muscle samples (~50 mg each) were pulverized in liquid nitrogen with an RNase-free mortar and pestle, and then dissolved using RNAiso Plus 9108 (Takara Bio) in accordance with the manufacturer's instructions. The concentration (ng/mL) and purity of the freshly extracted total RNA were measured using a microplate reader (Infinite M200 Pro; Tecan). The quality of the extracted total RNA was determined using spectrophotometry (NanoDrop 2000; Thermo Fisher Scientific). The concentration of total RNA was determined using fluorometry (Qubit 3.0; Life Technologies) and the integrity of total RNA was determined with an Agilent 2100 Nano 6000 Assay kit on an Agilent 2100 Bioanalyzer system. Only RNA samples with OD₂₆₀/OD₂₈₀ ratios within the range 1.8–2.2 and an RNA integrity number > 7.0 were used in the following experiments. Samples not meeting these criteria were excluded from the analysis (n = 4).

Sixty samples were submitted to Illumina for library preparation using a NEBNext Ultra RNA Library Prep Kit (poly-A selection; NEB) and fluorometer (Qubit 3.0; Life Technologies) to build and preliminarily quantify the library. The library was then quantified accurately (IQ SYBR GRN Kit; Bio-Rad) and the insert size of the library was determined (Agilent 2100 Bioanalyzer). The library of 60 samples was then subjected to pair-end 150-bp sequencing, aiming for coverage of 40–60 M reads, using a HiSeq PE Cluster Kit v4-cBot-HS on the HiSeq-2500 platform (Illumina). The off-machine data were converted into raw sequence reads after base recognition using bcl2fastq2 software (Love et al., 2014), and the results were stored in the FASTQ file format. The read quality was then inspected using FastQC and MultiQC (Zhou et al., 2018), and trimmed with Trimmomatic (version 3) (Bolger et al., 2014). These procedures yielded clean data for bioinformatics analysis.

Further bioinformatics data analysis was performed using the high-throughput data analysis software Chipster (version 3.16) (Kallio et al., 2011). The clean paired-end RNA-seq data obtained with Trimmomatic were aligned to the rn6 rat genome (HISAT2 version 2.1.0) (Kim et al., 2015) and annotated based on the Rattus_norvegicus.Rnor_6.0.95¹ file. In this study, the alignment files contained paired-end data. The mapped and aligned reads were quantified to determine gene-level counts using uniquely aligned unambiguous reads in HTSeq (version 0.6.1) (Anders et al., 2015) with the default settings and reverse strandedness. Raw counts were processed using the R bioconductor package DESeq2 (version 1.12.4 in R Studio version 3.6.3) (Love et al., 2014) and normalized using the DESeq method to remove library-specific artifacts. Among the 53 samples, total absolute read counts of <5 genes were considered to have low expression and were filtered out. DEGs between the experimental and control groups were calculated using the Wald test in DESeq2. Genes with log₂-fold change (FC; Injured/Control) > 1 or < −1 and adjusted p < 0.05 corrected for multiple testing using the Benjamini–Hochberg method (Benjamini and Hochberg, 1995) were considered significant and subjected to further downstream analysis.

WGCNA was performed to identify the gene co-expression modules at each post-injury time point using the WGCNA package in R Studio (version 3.6.3) (Langfelder and Horvath, 2008). For subsequent analysis to generate modules, the same parameters were used to construct all modules for each post-injury time point in independent analyses. A signed weighted correlation matrix containing pairwise Pearson correlation coefficients between all genes across all samples was computed using a soft threshold of β = 22 to attain a scale-free topology (Bogenpohl et al., 2019). Using this adjacency matrix, the topological overlap measure, which measures the network interconnectedness (Yip and Horvath, 2007), was calculated and used as an input to group highly correlated genes together using average linkage hierarchical clustering. The WGCNA dynamic hybrid tree-cut algorithm was used to detect network modules of co-expressed genes. After we have selected the aforementioned HCr modules (Jin et al., 2019), functional enrichment analysis of genes in the HCr modules was performed to detect enriched KEGG pathways and GO terms representing biological processes, and statistically significant clusters were identified using the Metascape² database (Zhou et al., 2019).

In gene co-expression networks, high-connectivity genes, referred to as hub genes (high degrees of connectivity), are critical for the maintenance of overall network stability. In WGCNA, intra-module connectivity and correlations with module eigengenes were used to select hub genes without any statistical criteria (Bi et al., 2015). STRING was used to construct a PPI network, with closely related genes located closer together. In the present study, the degrees of hub genes in each HCr module were calculated and visualized using Cytoscape version 3.7.2³. The top 30 genes according to degree of connectivity were screened as high-connectivity co-expressed genes in each HCr module.

Time-dependent marker gene analysis was performed using normalized log₂ gene expression values after DESeq analysis using the Next maSigPro R package (version 1.6.0) (Nueda et al., 2014). Statistical analysis of time-series data identified genes that exhibited changes in their expression over time and/or followed a specific expression pattern. Briefly, the p vector function was used to compute a regression fit for each gene in both control and separate experimental groups. Temporally DEGs were detected using the generalized linear model setting with p ≤ 0.05 after false discovery rate correction. The final selection of temporally DEGs was performed by filtering the results of the second regression model using the get siggenes function, with the R² parameter set to 0.7 and the vars parameter set to groups (Grilli et al., 2018). Finally, significant genes were grouped into k = 6 groups (set value).

As mentioned above, STRING was used to construct a PPI network to select the top 10% of genes in each cluster according to degree of connectivity, and the results were visualized using Cytoscape.

STEM (Short Time Series Expression Miner, version 1.3.12) (Ernst and Bar-Joseph, 2006) was used to identify the dynamic gene-expression clusters among the DEGs common to all seven experimental groups, and significantly enriched gene families with similar expression patterns were clustered according to the default parameters. Similar to how genes were clustered, we used STRING to construct a PPI network to select the top 10% of genes in each cluster according to degree of connectivity, and the results were visualized using Cytoscape. Functional GO term and KEGG pathway enrichment analyses of the significant clusters were performed using the Metascape database (see above).

The injured limbs from the experimental groups (n = 56) and control group (n = 8) were harvested and placed in 10% neutral buffered formalin for 24 h. The muscle tissue was rinsed and flushed with phosphate-buffered saline, gently squeezed to remove non-adherent bacteria, and embedded in paraffin. Then, longitudinal sections 5 μm thick were cut for histological observations and immune-histochemical processing. Histological sections for each time point were stained with H&E to observe the progression from damage to repair in the early period following SMI.

A primary antibody against MPO (1:500 dilution, ab9535; Abcam) was used with a Metal Enhanced DAB Chromogenic Kit (AR1026; Boster) to stain infiltrating neutrophils. A secondary anti-mouse/rabbit horseradish peroxidase immune-globin G antibody (SA1020; Boster) was used to identify sites of inflammation. Sections were counterstained with hematoxylin.

Following staining, the slides were imaged using a Tissue Fax Plus 2000 slide scanner (Tissuen Gnostics) at a magnification of 40× (H&E staining) or 20× (immune-histochemical staining).