The Whole Transcriptome Involved in Denervated Muscle Atrophy Following Peripheral Nerve Injury

Peripheral nerve injury (PNI) usually leads to progressive muscle atrophy and poor functional recovery. Previous studies have demonstrated that non-coding ribonucleic acid (ncRNA) is a key regulator of muscle atrophy and beneficial for the treatment of PNI. We aimed to analyze the whole transcriptome involved in denervated muscle atrophy after PNI. Animal models of sciatic nerve injury were assessed at 0 (control group), 1, 2, 4, and 8 weeks after injury. The expression patterns in the whole transcriptome in the gastrocnemius muscle were profiled using RNA sequencing at each time point and compared to that obtained in the control group. Six-hundred and sixty-four long non-coding RNAs, 671 microRNAs, 236 circular RNAs, and 12,768 messenger RNAs (mRNAs) were differentially expressed (DE) after injury. Changes in some of the DE ncRNAs and mRNAs were validated using quantitative polymerase chain reaction. Gene Ontology and Kyoko Encyclopedia of Genes and Genomes analysis revealed the potential functions of and relationships among the DE ncRNAs and mRNAs. To our knowledge, this is the first study to expound the whole transcriptome involved in denervated muscle atrophy, and provides a theoretical basis for further research targeting ncRNAs.


INTRODUCTION
Peripheral nerve injury (PNI) may result in developmental atrophy of the target skeletal muscle and poor functional recovery due to delayed surgery. Due to their slow rate of regeneration (1-3 mm/d), the injured axons usually require at least 3 months to regenerate to distal organs. However, the target muscle usually undergoes atrophy and cannot be reinnervated after the prolonged period required for axonal regeneration. This in turn leads to poor neurological functional outcomes (Höke and Brushart, 2010). The major effects of progression of this type of atrophy are a loss in muscle mass and decreased strength. In order to improve functional motor recovery after nerve injury, it is necessary to reduce the time required for innervation and to decelerate the progress of denervated muscle atrophy (Moimas et al., 2013). The progression of denervated muscle atrophy includes a series of pathological events and complex molecular mechanisms that are still largely unknown. Therefore, investigation of the underlying mechanisms of denervated muscular atrophy and the search for new therapeutics for this condition are particularly important. Our team has studied peripheral nerve repair and the protection of atrophic muscles for decades. The transcriptome has become an important field of study in the post-genomic era, and has unique advantages in the study of PNI (Phay et al., 2015). We therefore have great interest in the whole transcriptome involved in denervated muscle atrophy. Our aim is to better understand the genetic and neurobiological bases of denervated muscle atrophy in order to develop effective therapeutic strategies for this condition.
Significant molecular changes occur in the target muscle before it is reached by the regenerated axons (Sando and Cederna, 2016). Previous studies have shown that non-coding ribonucleic acids (ncRNA) plays a key role in regulating protein expression during muscle atrophy (Wilkins et al., 2004;Wang et al., 2011). In contrast to typical RNA, ncRNAs do not encode proteins, but instead functionally regulate the expression of proteins (Mattick and Makunin, 2006). ncRNAs may be subdivided based on length into small ncRNAs (<200 bp), which include long non-coding RNAs (lncRNAs) and microRNAs (miRNAs), and circular RNAs (circRNAs) >200 bp (Thum, 2014). The functional mechanisms of action of ncRNAs in denervated muscular atrophy have not been fully studied. To date, few cases have been reported in this field. Therefore, a systematic analysis of the expression of ncRNAs involved in the progression of denervated muscular atrophy is essential for the exploration of effective treatments and the prevention of this disorder. The transcriptome has a dynamic characteristic with highly specific features in space and time. Furthermore, the whole transcriptome is known to be a dynamic reflection of changes in ncRNAs and messenger RNAs (mRNAs). This can be used to further study the mechanisms underlying denervated muscle atrophy. Next-generation RNA sequencing (RNA-seq) is an extremely sensitive method for analyzing differentially expressed (DE) RNAs (Sîrbu et al., 2012). RNA-seq is superior to the use of gene microarrays, as it is not limited to only a subset of known RNAs, but can also be used to analyze the entire transcriptome in tissues and organs to provide biological information regarding possible functions of annotated or novel genes (Hurd and Nelson, 2009).
Here we first evaluated the expression profiles of ncRNAs and mRNAs in denervated muscle atrophy after PNI using RNA-seq. We detected significant expression changes in ncRNAs in the gastrocnemius muscle. Our findings provide a novel information platform that can be used to identify target genes involved in the mechanisms underlying denervated muscle atrophy. In addition, we utilized Gene Ontology (GO) and Kyoko Encyclopedia of Genes and Genomes (KEGG) analysis to identify the biological processes and pathways associated with the DE genes, such as the mitogen-activated protein kinase (MAPK) signaling pathway. These data will serve as useful resources for the development of therapeutic strategies or novel diagnostics in the future.

Animal Model and Tissue Collection
Fifteen healthy female C57BL/6 mice (weighing 22-25 g and aged 6-8 weeks) were randomly assigned to five experimental groups. Three mice were assigned to each group of mice assessed at 0 (control group), 1, 2, 4, and 8 weeks after sciatic nerve injury. The mice were deeply anesthetized with an intraperitoneal injection of 10% chloral hydrate (3 ml/kg). Once the mice were completely anesthetized, a 0.5 cm-long incision was made over the sciatic nerve of the left hind limb. The sciatic nerves were then exposed by gentle separation and cut approximately 1 cm proximal to the division of the tibial and common peroneal nerves. The incisions were closed using 4-0 sutures. All procedures were performed under sterile conditions. To evaluate the level of atrophy due to the injury, the animals were euthanized 0, 1, 2, 4, and 8 weeks post-injury. The gastrocnemius muscles were harvested from the operated experimental limb and processed for the next step.

RNA Isolation and Sequencing Library Preparation
Total RNA was isolated using the Trizol (Invitrogen, Carlsbad, CA, United States) according to the protocol. The quality of the purified RNA was assessed using a BioAnalyzer 2100 (Agilent Technologies; Santa Clara, CA, United States). A NanoPhotometer spectrophotometer (Implen; Westlake Village, CA, United States) was used to determine the quantity of purified RNA by measuring the absorbance at 260/280 nm. All of the RNA samples were stored at −80 • C and were sequenced by Novogene using the Illumina platform (Beijing, China).
A total of 8 µg of RNA per sample (3 µg for lncRNA and 5 µg for circRNA) was used to generate libraries for lncRNA and circRNA sequencing. After ribosomal RNAs were depleted using the Epicentre Ribo-zero TM rRNA Removal Kit (Epicentre, Madison, WI, United States), the rRNA-depleted RNAs were treated with RNase R (Epicenter) and extracted using Trizol (Invitrogen). Sequencing libraries were generated using the NEBNext R Ultra TM Directional RNA Library Prep Kit for Illumina R (NEB, Ipswich, MA, United States) according to the protocol. First, performed by the divalent cations in NEBNext First Strand Synthesis Reaction Buffer at high temperature. First strand cDNA was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (RNase H-). DNA Polymerase I and RNase H were then used to synthesize the second strand cDNA. The deoxythymidine triphosphates in the deoxynucleotide mixture in the reaction buffer were substituted with deoxyuridine triphosphates. The remaining overhangs were then made into blunt ends by exonuclease/polymerase activity. The NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization following the adenylation of the 3 ends of the DNA fragments. In order to select cDNA fragments 150-200 bp in length, the library fragments were purified using the AMPure XP system (Beckman Coulter; Beverly, CA, United States). After size selection using 3 µl of USER Enzyme (NEB), the cDNA was ligated at 37 • C for 15 min and then incubated at 95 • C for 5 min. Polymerase chain reaction (PCR) was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer.
A total of 3 µg RNA per sample was used to prepare the library for miRNA sequencing. Sequencing libraries were generated using the NEBNext R Multiplex Small RNA Library Prep Set for Illumina R (NEB) according to the manufacturer's instructions. Index codes were used to attribute sequences to the samples. First, the NEB 3 SR Adaptor was directly ligated specifically to the 3 ends of miRNAs. Following the 3 ligation reaction, some 3 SR Adaptor remained free. At this time, SR RT Primer was used to hybridize the excess 3 SR Adaptor. The single-stranded DNA adaptor was then transformed into a double-stranded DNA molecule in order to prevent the formation of adaptor dimers. In the subsequent ligation step, the 5 SR Adaptor did not ligate double-stranded DNAs, which do not undergo ligation mediated by T4 RNA Ligase 1. The 5 ends of the adapter were then ligated to the 5 ends of the miRNAs. M-MuLV Reverse Transcriptase (RNase H-) was used to synthesize the first strand cDNA. PCR amplification was performed using the LongAmp Taq 2X Master Mix, SR Primer for Illumina, and Index (X) Primer. An 8% polyacrylamide gel was used to purify the PCR products (100 V, 80 min), and 8 µL of elution buffer was used to dissolve DNA fragments 140-160 bp in length DNA fragments. Finally, the Agilent Bioanalyzer 2100 system was used to purify the products and evaluate the quality of the library.

Clustering and Sequencing
Index-coded sample clustering was performed using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) on a cBot Cluster Generation System following the manufacturer's protocol. An Illumina HiSeq X platform was used to sequence the library preparations and generate the 125/150 bp paired-end reads and the 50 bp single-end reads when the cluster generation was complete.

GO and KEGG Enrichment Analysis of DE Genes
Gene Ontology enrichment analysis of the DE genes was performed using the clusterProfiler R package to revise the gene length bias. DE genes with corrected P-values less than 0.05 were considered enriched. KEGG is a database resource for exploring high-level functions of biological systems at the molecular level using genome sequencing or other high-throughput experimental technologies. We used the clusterProfiler R package to assess the statistical enrichment of DE genes in KEGG pathways.

Statistical Analysis
All data are expressed as means ± standard errors of the mean. One-way analysis of variance was used to statistically analyze muscle weight. The qPCR results were analyzed using Student's t-tests. P-values less than 0.05 were considered statistically significant.

General Observation of Muscle Atrophy
The volumes of the gastrocnemius muscle in the injury groups were less than those in the control group after nerve surgery (Figure 1). The muscle weight on the experimental side is presented as a percentage of the muscle weight on the control side in each group. This ratio was significantly lower in the surgery groups than in the control group ( Table 1). The gastrocnemius muscle on the injured side therefore had obvious atrophy after PNI, indicating that our denervated muscle atrophy animal models were successfully established.

Quality Control and Correlations of the Data
General information regarding RNA-seq data quality was obtained by calculating the ratio of clean data to all data (Figure 2A). The ratio of clean reads as a percentage of all raw reads was greater than 93% for all samples. This suggests that the data were high-quality and suitable for subsequent analysis. Correlations were analyzed using the Pearson correlation   coefficient, as calculated in R ( Figure 2B). The R 2 values at each time point were larger than 0.92, suggesting that the biological replicates of the samples in each group met the requirements for subsequent analysis.

Validation of ncRNA and mRNA Expression Using qPCR
We confirmed the accuracy of sequencing data for selected ncRNAs and mRNAs (ENSMUSG00000087523.1, LNC_000280,  mmu-miR-34c-5p, mmu-miR-142a-3p, mmu_circ_0001068, novel_circ_0011051, Myl4, and Myom3) that were DE after injury using qPCR. The mRNA validation results were consistent with the RNA-seq data (Figure 6). These convincing sequencing data can be used as the basis for future research.

Functional Annotation Based on GO
To predict the functions of the mRNAs in PNI, genes with absolute correlation values >0.95 were selected and underwent GO and KEGG pathway analysis. The GO enrichment analysis results for the DE lncRNAs, miRNAs, circRNAs, and mRNAs are shown in Figures 7-10, respectively. Directed acyclic   Based on the co-expressed genes with the DE lncRNAs, the most significantly enriched BP involved the generation of precursor metabolites and energy at 1, 2, 4, and 8 weeks after the injury. The most involved BP were related to single-organism metabolic processes at 1, 2, 4, and 8 weeks after the injury. The most significantly enriched CC was the mitochondrion at 1, 2, 4, and 8 weeks after injury. The most involved CC were the cell at 1 and 4 weeks after injury, intraCCs at 2 weeks after the injury, and the cytoplasm at 8 weeks after the injury. The most significantly enriched MF were protein binding at 1 week after the injury, phosphoric ester hydrolase activity at 2 weeks after the injury, and cofactor binding at 4 and 8 weeks after injury. The most involved MF were binding at 1, 2, and 4 weeks after injury, and protein binding at 8 weeks after injury (Figure 7).
Based on the target genes of the DE miRNAs, the most significantly enriched BP were the positive regulation of biological processes at 1 and 2 weeks after injury, localization at 4 weeks after injury, and developmental processes at 8 weeks after injury. The most involved BP were metabolic processes at 1, 2, 4, and 8 weeks after the injury. The most significantly enriched CC were membrane-bounded organelles at 1 and 8 weeks after injury, and intracellular componenets at 2 and 4 weeks after injury. The most involved CC was the cell at 1, 2, 4, and 8 weeks after injury. The most significantly enriched MF were protein binding at 1, 4, and 8 weeks after injury, and binding at 2 weeks after injury. The most involved MF was binding at 1, 2, 4, and 8 weeks after injury (Figure 8).  Frontiers in Molecular Neuroscience | www.frontiersin.org FIGURE 6 | Validation of DE ncRNAs and mRNAs by qPCR. The results confirmed that the relative expression levels of ENSMUSG00000087523.1, LNC_000280, mmu-miR-34c-5p, mmu-miR-142a-3p, mmu_circ_0001068, novel_circ_0011051, Myl4, and Myom3 were significantly regulated after injury. * P < 0.05 compared to the control group. CON, control; INJ, injured.  Based on the target genes of the DE circRNAs, the most significantly enriched BP were muscle contraction at 1 and 4 weeks after injury, and muscle system processes at 2 and 8 weeks after the injury. The most involved BP were muscle contraction at 1 and 4 weeks after injury, and muscle system processes at 2 and 8 weeks after injury. The most significantly enriched CC were myosin filaments at 1 week after injury, and myofibrils at 2, 4, and 8 weeks after injury. The most involved CC were myofibrils at 1, 2, 4, and 8 weeks after injury. The most significantly enriched MF were calmodulin binding at 1 week after the injury, adenosine triphosphate binding at 2 weeks after injury, microfilament motor activity at 4 weeks after injury, and actin binding at 8 weeks after the injury. The most involved MF were hydrolase activity at 1 and 2 weeks after injury, adenosine triphosphate binding at 4 weeks after injury, and cytoskeletal protein binding at 8 weeks after injury (Figure 9).
Based on the target genes of the DE mRNAs, the most significantly enriched BP were the negative regulation of biological processes at 1 week after injury, positive regulation of biological processes at 2 and 4 weeks after injury, and single-organism metabolic processes at 8 weeks after injury. The most involved BP were cellular processes at 1, 2, and 4 weeks after injury, and developmental processes at 8 weeks after the injury. The most significantly enriched CC was the cytoplasm at 1, 2, 4, and 8 weeks after injury. The most involved CC were the cell at 1, 2, and 4 weeks, and intraCCs at 8 weeks after the injury. The most significantly enriched MF was protein binding at 1, 2, 4, and 8 weeks after surgery. The most involved MF was binding at 1, 2, 4, and 8 weeks (Figure 10).

Pathway Prediction Based on KEGG
lncRNAs have extensive regulatory functions that can directly regulate the structure of DNA and the transcription and translation of RNA. They can also inhibit target gene regulation by miRNAs to indirectly regulate gene expression. An enriched scatter diagram was used to display the KEGG enrichment analysis results for the intersection of genes co-localized and co-expressed with the DE lncRNAs, target genes of the DE miRNAs, and predicted mRNAs. The degree of KEGG enrichment was reflected by the Rich factor, the Q-value, and the number of genes. Based on the KEGG results, the top 5 significantly involved pathways in the pathogenesis of denervated muscle atrophy were osteoclast differentiation, MAPK signaling, carbon metabolism, the proteasome, and leishmaniasis 1 week after the injury; osteoclast differentiation, focal adhesion, carbon metabolism, MAPK signaling, and hypoxia-induced factor-1 signaling 2 weeks after the injury; carbon metabolism, citric acid cycle, leishmaniasis, nonalcoholic fatty liver disease, and oxytocin signaling 4 weeks after injury; and non-alcoholic fatty liver disease, carbon metabolism, Parkinson's disease, oxidative phosphorylation, and Alzheimer's disease 8 weeks after the injury (Figure 11). The MAPK signaling pathway was the most involved pathway in all groups based on the number of participating genes. The genes involved in the MAPK signaling pathway might play crucial roles in denervated muscle atrophy and are thus worthy of further research. These results illustrate the regulatory relationships between ncRNAs and mRNAs in the mechanisms underlying denervated muscle atrophy after PNI. Further in-depth studies by our team will be forthcoming.

DISCUSSION
Peripheral nerve injury, especially that involving long nerve segments, may result in poor functional recovery such as progressive skeletal muscle atrophy, which occurs until axons regenerate into the target tissue. Previous studies have shown that some ncRNAs and mRNAs can regulate the progress of nerve regeneration in response to muscle atrophy (Pan et al., 2015;Huang et al., 2016). This is the first study to investigate the expression profile of the whole transcriptome, including ncRNAs and mRNAs, during skeletal muscle atrophy after PNI using sequencing analysis. We also used GO and KEGG pathway analyses to predict the potential functions and mechanisms of the DE RNAs. Our results support the opinion that ncRNAs play key roles in denervated muscle atrophy pathogenesis, and may highlight a potential therapeutic target for the treatment of this condition.
Peripheral nerve injury is a common clinical disease. Incomplete nerve regeneration leads to poor recovery, such as loss of limb movement. Long periods of denervation lead to several progressive degenerative processes, resulting in atrophy of the denervated muscle. In the absence of axonal regeneration, the target muscles lose volume and undergo fibrosis. In addition, the number of myofibrils is reduced and they decrease in size, while the number of receptive motor endplates is decreased (Carlson et al., 1996;Ma et al., 2011). Previous studies from our group have been performed to accelerate axonal regeneration using experimental therapies (Wei et al., 2009;Kou et al., 2013). Today, the ideal therapy would also enable the maintenance of muscle functions during the regeneration process. Kobayashi et al. (1997) have reported significant changes in both muscle mass and whole motor function when the nerve regeneration interval is longer than 1 month. Thus, therapeutic prevention of muscle atrophy following PNI may be beneficial in these cases. Recent studies support the opinion that early intervention is appropriate for denervated muscles when reinnervation is delayed, or in the worst case, when there is no opportunity for reinnervation (Gao et al., 2005;Deshpande et al., 2006). Therefore, the identification of molecular mechanisms for the prevention of skeletal muscle atrophy is necessary.
ncRNAs have been shown to play important roles in the regulation of specific target mRNAs and proteins in many diseases, including muscle atrophy (Crist and Buckingham, 2010;Russell et al., 2013). Eisenberg et al. (2007) have identified 185 miRNAs that are significantly regulated in 10 major muscular disorders. Greco et al. (2009) have reported that some miRNAs related to regeneration and degeneration were modified in both patients with muscular dystrophy and in mouse models. The above studies indicate the important role of miRNAs in pathophysiological pathways in muscle atrophy. However, no studies have investigated ncRNAs, which include lncRNAs and circRNAs, during denervated muscle atrophy. We utilized RNA-seq to detect and further analyze the biological roles of DE ncRNAs and mRNAs in denervated gastrocnemius muscles. Six-hundred and sixty-four lncRNAs, 671 miRNAs, 236 circRNAs, and 12,768 mRNAs were significantly DE during muscle atrophy. Our findings regarding ncRNAs may provide a theoretical basis for the study of the mechanisms underlying denervated muscle atrophy. Due to the special regulatory functions of lncRNAs, it is worth mentioning that 73 DE lncRNAs were detected during the progression of atrophy after PNI. These lncRNAs should be considered potential therapeutic targets in further studies.
Gene Ontology analysis is currently the main bioinformatics tool used to unify the representation and product attributes of genes (Altshuler et al., 2008). It has been confirmed that GO terms and GO annotations can be used to predict gene functions and trends (Du et al., 2016). The KEGG pathway database is currently widely used as an enrichment analysis platform for gene functions due to the high levels of functional information that can be used for systematic analysis (Camon et al., 2004). We investigated gene functions and pathways related to the DE ncRNAs and mRNAs in the gastrocnemius muscle after PNI using GO and KEGG analyses. Our results indicated that the most significantly involved pathway in the pathogenesis of denervated muscle atrophy was the MAPK signaling pathway. Other pathways involved in denervated muscle atrophy, such as carbon metabolism, the proteasome, focal adhesion, the citric acid cycle, and oxidative phosphorylation are all fundamental biological processes present in all organisms. The results of the GO and KEGG analyses support investigation of the underlying mechanisms of denervated muscle atrophy after PNI. The accuracy of the sequencing data as the basis for future research was confirmed. It is well known that higher biological replicate numbers can minimize the false positive rate in RNA-seq results. A limitation of our study was that we tested three samples in parallel in each group to obtain the minimum number of biological replicates. However, we strictly controlled other variables (animal genetic background, feeding conditions, age, sex, and external shape). Furthermore, the correlation analysis results suggested that the individual differences in the samples were small in each group.
Our study is the first to investigate the expression profile of the whole transcriptome in denervated muscle atrophy using innovative RNA-seq analysis. Here we provide a novel bioinformatics platform for the study of the roles of ncRNAs and mRNAs in denervated muscle atrophy, which may be benefit future research. This study supports the opinion that ncRNAs, and especially lncRNAs, may have interactions and regulatory relationships with their co-localized and coexpressed genes, which may in turn affect the pathogenesis of denervated muscle atrophy. However, the regulatory roles of ncRNAs in the pathogenesis of denervated muscle atrophy are complicated and largely unknown. Our team will further investigate the related genes and corresponding signaling pathways for the predicted ncRNAs and mRNAs in a future study to elucidate the mechanisms underlying denervated muscle atrophy.

ETHICS STATEMENT
The study was approved by the Animal Experiment Ethics Committee of Peking University People's Hospital. All of the procedures were performed in accordance with the relevant guidelines and regulations.