Longitudinal RNA-Seq analysis of acute and chronic neurogenic skeletal muscle atrophy

Skeletal muscle is a highly adaptable tissue capable of changes in size, contractility, and metabolism according to functional demands. Atrophy is a decline in mass and strength caused by pathologic loss of myofibrillar proteins, and can result from disuse, aging, or denervation caused by injury or peripheral nerve disorders. We provide a high-quality longitudinal RNA-Seq dataset of skeletal muscle from a cohort of adult C57BL/6J male mice subjected to tibial nerve denervation for 0 (baseline), 1, 3, 7, 14, 30, or 90 days. Using an unbiased genomics approach to identify gene expression changes across the entire longitudinal course of muscle atrophy affords the opportunity to (1) establish acute responses to denervation, (2) detect pathways that mediate rapid loss of muscle mass within the first week after denervation, and (3) capture the molecular phenotype of chronically atrophied muscle at a stage when it is largely resistant to recovery.


Background & Summary
Skeletal muscle atrophy is the loss of muscle mass and function that occurs in response to diverse stimuli including disuse/immobility, glucocorticoid treatment, cancer, aging, and denervation [1][2][3][4][5] . Biologically, atrophy reflects the active loss of skeletal muscle contractile proteins, leading to loss of strength and functional impairment with substantial impact on quality of life and, in some cases, reduced survival [6][7][8] . In addition, chronically denervated, atrophied muscle shows impaired capacity for reinnervation and functional recovery, which significantly limits prospects for recovery in settings of chronic neuromuscular disease, delayed repair, or large nerve lesions 9-12 . Nerve-evoked contraction is the most important factor for maintaining or regaining muscle mass and force 13 . Neurogenic atrophy refers specifically to skeletal muscle atrophy resulting from denervation, as may occur in traumatic injury or diseases that affect the peripheral nervous system, such as amyotrophic lateral sclerosis (ALS) [14][15][16][17] . A number of "atrogenes" are induced as a result of denervation and in response to various triggers of muscle atrophy; among these are specific ubiquitin ligases targeting components of the sarcomere [18][19][20][21][22][23][24][25][26][27][28][29] . A comprehensive analysis of the global gene pathways that change in response to denervation and during atrophy may offer an optimal chance of identifying means to pharmacologically maintain or increase muscle mass and function in atrophy-associated disease states.
We provide here a comprehensive RNA-Seq dataset 30 to identify gene expression changes across the entire longitudinal course of muscle atrophy, affording the opportunity to (1) establish acute responses to denervation within the first day, (2) detect pathways that mediate rapid proteolysis and loss of muscle mass within the first week after denervation, and (3) capture the molecular phenotype of chronically atrophied muscle (weeks to months after denervation) at a stage when it is largely resistant to reinnervation and recovery.
We generated a longitudinal RNA-Seq dataset from a cohort of adult (8-week-old) wild type C57BL/6 J male mice denervated for 0 (baseline), 1, 3, 7, 14, 30, or 90 days (n = 4 for each timepoint) 30 . We elected to use tibial nerve transection as a model for muscle denervation, as this approach is physiologically meaningful while limiting the morbidity (i.e., pain and immobility) associated with complete sciatic nerve transection 31 . The tibial nerve is the largest branch of the sciatic nerve that supplies skeletal muscles of the posterior compartment of the lower limb, including the gastrocnemius and soleus. In brief, we identified and separated the tibial nerve from other branches of the sciatic nerve, then ligated, cut distally, and sutured the proximal stump in place to prevent muscle reinnervation during chronic studies. We have established that this model reliably induces significant gastrocnemius atrophy within one week after denervation, with atrophy becoming progressively more severe over time (Fig. 1c,d).
The samples collected and described in this manuscript include transcriptional profiles from a total of 28 denervated gastrocnemii and 28 contralateral (paired) intact gastrocnemii, comprising 4 denervated and 4 contralateral (paired) intact gastrocnemii for each of 7 denervation durations [0 (baseline), 1, 3, 7, 14, 30, and 90 days] 30 . All specimens were generated from a cohort of male C57BL/6 J mice that were 8 weeks of age at the start of the study. These data provide a comprehensive description of baseline gene expression in adult mouse skeletal muscle and a broad assessment of the acute and longitudinal gene expression changes in atrophying muscle associated with denervation.

Methods
Animal husbandry. 8-week-old C57BL/6 J male mice (Stock #000664) were obtained from the Jackson Laboratory (Bar Harbor, ME) and randomized into 7 groups of n = 4 mice per group for the following denervation timepoints: 0, 1, 3, 7, 14, 30, and 90 days. Animal subjects were housed in a controlled environment with a 12:12-h light-dark cycle with ad libitum access to water and food (Envigo 2018 SX). All mouse experiments were carried out under protocols approved by the JHU Animal Care and Use Committee.
Tibial nerve denervation surgery. Mice were anesthetized with 1.5% isoflurane/2% oxygen using a VetEquip inhalation system (Livermore, CA). The left hindlimb was shaved and sterilized, and a 1 cm incision was introduced in the skin overlying the dorsal thigh. Myofascial planes were gently separated to reveal the sciatic nerve. The tibial nerve branch was identified at its distal branch point and gently separated from the sciatic and peroneal nerves, then ligated proximally and distally using a 10-0 polyamide monofilament suture. The tibial nerve was then transected, the nerve length between ligatures carefully resected, and the proximal stump sutured to the biceps femoris muscle to prevent distal reinnervation. The incision was then closed using stainless steel wound clips. Mice were monitored for recovery from anesthesia and then returned to their home cages.
Myofiber morphometry. Gastrocnemii   Overview of the experimental procedure. The tibial nerve, the largest branch of the sciatic nerve, supplies the gastrocnemius muscle and other muscles of the lower limb posterior compartment. In our mouse model of denervation atrophy, the sciatic nerve is identified, and its branches separated to isolate the tibial nerve (a; nerve identities are as follows: 1, sural nerve; 2, tibial nerve; 3, common peroneal/fibular nerve; 4, sciatic nerve). We generated a cohort of C57BL/6 J male mice denervated for 0, 1, 3, 7, 14, 30, or 90 days (b,c). Significant atrophy is apparent by 7 days after denervation, with consistent decline in mass during chronic denervation (d); ***P < 0.001 compared to baseline.
www.nature.com/scientificdata www.nature.com/scientificdata/  Gastrocnemius myofiber morphometry. Atrophy of type I, IIa, and IIb myofibers was analyzed by assessment of minimum Feret diameter at baseline (t = 0 days) and 7, 14, 30, and 72 days post-denervation. All three myofiber types showed significant atrophy within the first week after denervation, with the greatest change in magnitude observed for type IIb myofibers overall. Scale bar, 100 μm. RNA-Seq library preparation, sequencing, and bioinformatics analysis. RNA-sequencing was carried out using TrueSeq RiboZero gold (stranded) kit (Illumina, catalogue #20020597). Libraries were indexed and sequenced over 18 lanes using HiSeq4000 (Illumina) with 69-bp paired end reads. Quality control (QC) was performed on base qualities and nucleotide composition of sequences using FastQC version 0.11.5 34 , to identify problems in library preparation or sequencing. Sequence quality for the dataset described here was sufficient that no reads were trimmed or filtered before input to the alignment stage. Paired-end reads were aligned to the most recent Mus musculus mm10 reference genome (GRCm38.75) using the STAR spliced read aligner (version 2.4.0) 35 . Average input read counts were 58.0 M per sample (range 39.1 M to 91.0 M) and average percentage of uniquely aligned reads was 81.9% (range 72.3% to 88.6%). Total counts of read-fragments aligned to known gene regions within the mouse (mm10) refSeq (refFlat version 07.24.14) reference annotation were used as the basis for quantification of gene expression. Fragment counts were derived using HTSeq (version 0.6.0) and the mm10 refSeq transcript model 36 . Low count transcripts were filtered, and count data were normalized using the method of trimmed mean of M-values (TMM) 37 followed by removing unwanted variation using Bioconductor package RUVseq 38 with k value of 1. Differentially expressed genes (FDR < 0.1) were then identified using the Bioconductor package limma with voom function to estimate mean-variance relationship, followed by empirical www.nature.com/scientificdata www.nature.com/scientificdata/ Bayes moderation [39][40][41] . Pairwise comparisons between denervated and contralateral intact muscle at each timepoint were used as the basis for model contrasts. All bioinformatics analyses were conducted using R version 3.5.1 42 .

Data Records
Sequencing data in the fastq format have been deposited in NCBI Sequence Read Archive (SRA) 30 . A metadata table (Supplementary Table S1) is available with details for each sample.

technical Validation
Reproducible skeletal muscle atrophy using tibial nerve denervation model. Tibial nerve denervation resulted in a reliable time-dependent loss of skeletal muscle mass, with a significant difference in mass between denervated and contralateral intact gastrocnemii detected by day 7 post-denervation (Fig. 1c,d). All mice used in this study entered the cohort at the same time, with sequential denervation according to the designated timepoints, to remove age as a potential confounding variable. Mouse gastrocnemius contains a mixed population of myofiber types including so-called slow twitch myofibers (type I) and fast twitch myofibers (type IIa and IIb). After muscle denervation, all three of these myofiber populations showed a significant reduction in size as measured by minimum Feret diameter, with the most substantial rate of individual myofiber atrophy occurring within the first two weeks post-denervation (Fig. 2). Type IIb myofibers, the most abundant myofiber www.nature.com/scientificdata www.nature.com/scientificdata/ type in mouse gastrocnemius, showed the largest magnitude of atrophy (Fig. 2f). Multiple linear regression with myofiber type, myofiber type-time interactions, and time modeled with a spline at t = 14 days was used to model rates of atrophy among type I, IIa, and IIb myofibers; bootstrapping was used to estimate standard errors. Results are presented in Table 1.
RNA quality control. RNA integrity was analyzed using an Agilent 2100 Bioanalyzer (Fig. 3). The mean RNA Integrity Number (RIN) for RNA isolated from denervated and contralateral intact gastrocnemii was 7.8 ± 0.3 and 8.3 ± 0.1 (mean ± SEM), respectively, with no significant difference in RIN by denervation status.
Read quality and base-calling accuracy. Read quality was high with Phred quality score >70 for the majority of the cycles, and lower quartile base qualities were generally high (Fig. 4). No reads or samples necessitated exclusion based on read quality. The nucleotide composition patterns (proportions of A/C/G/T) of all samples were as expected, with nearly uniform proportions of each nucleotide across sequencing cycles (with the exception of a non-random pattern of nucleotide proportions in the first 13 sequencing cycles as a result of random hexamer priming) (Fig. 5). No read trimming or filtering was required because the quality distribution and variance appeared normal across all reads and samples. Tables 2-9. Similar sequencing depths and mapping rates were observed for the denervated and contralateral intact skeletal muscle samples.     Table 3. Day 0 (baseline) alignments.
Unsupervised clustering analysis of longitudinally denervated samples. Multidimensional scaling using expression levels of all genes demonstrated temporal clustering based on denervation status, with replicates within each denervation timepoint clustering closer to each other than to other denervation timepoints (Fig. 7).   www.nature.com/scientificdata www.nature.com/scientificdata/ time-dependent comparison of denervated and contralateral intact skeletal muscle transcriptomes. Normalized gene counts from denervated and contralateral intact skeletal muscle at each timepoint are compared in scatter plots (Fig. 8).
Differential expression analysis. MA-plots showing the log-fold change (M-values, the log of the ratio of counts for each gene across the two samples being compared) against the normalized log-average (A-values, the average counts for each gene across the two samples being compared) indicate substantial differences in gene expression in skeletal muscle during acute and chronic neurogenic atrophy (Fig. 9a-g). Volcano plots indicate minimal differences in gene expression at baseline (intact muscle) (Fig. 9h), but demonstrate that thousands of genes are significantly differentially expressed (FDR < 0.1) within the first day after denervation (Fig. 9i) and beyond (Fig. 9j-n). A summary of the number of differentially expressed genes at each timepoint is shown in Fig. 9o.  www.nature.com/scientificdata www.nature.com/scientificdata/ Fig. 9 Differential expression analysis. MA-plots comparing the log2 fold change of gene expression for denervated vs. contralateral intact skeletal muscle at each timepoint plotted against the normalized average of the counts (a-g). Volcano plots showing the -log10 FDR for difference in expression between denervated and contralateral intact skeletal muscle for each gene detected, plotted against the log2 fold-change (h-n). Genes with FDR < 0.1 are depicted in red. The total number of significantly differentially expressed genes (FDR < 0.1) at each timepoint is summarized in panel (o).

Usage Notes
The RNA-Seq dataset presented in this study provides a detailed view of the acute and chronic gene expression changes that occur in denervated, atrophying skeletal muscle. These data may provide insight into the early events associated with acute loss of neuronal input that trigger rapid atrophy, as well as the gene expression changes in chronically denervated and severely atrophied skeletal muscle associated with impaired capacity for reinnervation. Defining these changes may afford opportunities to limit the rate and severity of skeletal muscle atrophy, and to enhance functional reinnervation.

Code availability
Scripts used in the RNA sequencing analyses are available at https://github.com/icnn/RNAseq-PIPELINE.git.