Original Research ARTICLE
Microarray Analysis of Gene Expression Provides New Insights Into Denervation-Induced Skeletal Muscle Atrophy
- 1Key Laboratory of Neuroregeneration of Jiangsu, Ministry of Education, Jiangsu Clinical Medicine Center of Tissue Engineering and Nerve Injury Repair, Co-Innovation Center of Neuroregeneration, Nantong University, Nantong, China
- 2The Second Affiliated Hospital of Nantong University, Nantong University, Nantong, China
- 3Department of Surgery, Changshu Affiliated Hospital of Nanjing University of Chinese Medicine, Changshu Traditional Chinese Medicine Hospital, Changshu, China
- 4Department of Medical Laboratory, School of Public Health, Nantong University, Nantong, China
- 5Department of Orthopedics, Affiliated Hospital of Nantong University, Nantong, China
Denervation induces skeletal muscle atrophy, accompanied by complex biochemical and physiological changes, with potentially devastating outcomes even an increased mortality. Currently, however, there remains a paucity of effective strategies to treat skeletal muscle atrophy. Therefore, it is required to understand the molecular mechanisms of skeletal muscle atrophy and formulate new treatment strategies. In this study, we investigated the transcriptional profile of denervated skeletal muscle after peripheral nerve injury in rats. The cDNA microarray analysis showed that a huge number of genes in tibialis anterior (TA) muscles were differentially expressed at different times after sciatic nerve transection. Notably, the 24 h of denervation might be a critical time point for triggering TA muscle atrophy. According to the data from self-organizing map (SOM), Pearson correlation heatmap, principal component analysis (PCA), and hierarchical clustering analysis, three nodal transitions in gene expression profile of the denervated TA muscle might partition the period of 0.25 h–28 days post nerve injury into four distinct transcriptional phases. Moreover, the four transcriptional phases were designated as “oxidative stress stage”, “inflammation stage”, “atrophy stage” and “atrophic fibrosis stage”, respectively, which was concluded from Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene ontology (GO) analyses at each transcriptional phase. Importantly, the differentially expressed genes at 24 h post sciatic nerve transection seemed to be mainly involved in inflammation, which might be a critical process in denervation-induced muscle atrophy. Overall, our study would contribute to the understanding of molecular aspects of denervation-induced muscle atrophy, and may also provide a new insight into the time window for targeted therapy.
Skeletal muscle, as a highly plastic organ, is generally capable of adapting to mild external changes through secretion of myokines and myometabolites, thus maintaining whole body homeostasis (Aversa et al., 2019; Oost et al., 2019). However, a diverse array of pathological stimuli, including disuse, nutrition deprivation, denervation, and systematic disease (e.g., AIDS, sepsis, renal and cardiac failure, muscular dystrophies, or cancer), are likely to induce skeletal muscle atrophy, characterized by muscle weakness and dysfunction (Bonaldo and Sandri, 2013; Szentesi et al., 2019). Skeletal muscle atrophy may cause adverse impacts on the quality of life and life prognosis in patients. Unfortunately, currently available treatments are often unable to achieve desirable therapeutic outcomes (Qiu et al., 2018). In this sense, it is urgently required to better understand molecular mechanisms of skeletal muscle atrophy and actively seek new therapeutic strategies.
During skeletal muscle atrophy, a series of biochemical and physiological alterations appear in atrophic muscle, triggering the gene expression changes. Bodine et al. (2001) and Bodine and Baehr (2014) demonstrate that muscle RING finger 1 (MuRF1) and muscle atrophy F-box (MAFbx), two muscle-specific E3 ubiquitin ligases, are up-regulated in different models of muscle atrophy, playing a pivotal role in atrophy programs; Wu et al. (2014) indicate that NF-kB sites are necessary for MuRF1 promoter activation in disuse-induced muscle atrophy. Recently, an increasing number of studies have adopted cDNA microarray to reveal that many genes are differentially expressed in atrophic muscles (Gomes et al., 2001; Kunkel et al., 2011; Fedorov et al., 2014; Kruger et al., 2015; Jeong et al., 2017; Weng et al., 2018). Although these findings provide extensive gene expression data, their design of time points is not comprehensive enough to fully cover the entire process of muscle atrophy, especially not including an early stage. A global analysis is necessary for obtaining more valuable information on gene expression profiles of skeletal muscle in response to different atrophy-inducing stimuli.
Animals with sciatic nerve transection have been commonly chosen as a study model for denervation-induced skeletal muscle atrophy (Madaro et al., 2018; Lala-Tabbert et al., 2019). Once sciatic nerve is transected, target muscles will lose their function of “muscle pump” due to the loss of nerve innervation, which leads to relatively reduced perfusion of target muscles, facilitating skeletal muscle atrophy (Laughlin, 1987; Marra et al., 2015). This type of skeletal muscle atrophy is regulated by different molecular mediators (Tang et al., 2014; Li et al., 2017; Choi et al., 2018; Pigna et al., 2018; Lala-Tabbert et al., 2019), but the issue of how molecular mediators are deliberately orchestrated to activate atrophic programs is to be systematically addressed. The purpose of this study was to furnish a global perspective of transcriptional regulation for denervation-induced skeletal muscle atrophy. We constructed a model of denervated skeletal muscle atrophy by using the tibialis anterior (TA) muscles of rats which underwent sciatic nerve transection, and performed cDNA microarray to identified thousands of differentially expressed genes (DEGs) in TA muscles at different times after sciatic nerve transection. The microarray data were then analyzed by clustering and bioinformatic methods to visualize some interesting aspects of gene regulation in atrophic skeletal muscles during the denervation period.
Materials and Methods
A total of 72 adult male Sprague-Dawley (SD) rats (weighing ∼200 g) were provided by the institutional Experimental Animal Center (Nantong, China), and were maintained under temperature 22°C and a 12-h light/dark cycle with free access to food and water. The experiments involving animals were carried out in accordance with the animal care guidelines of Nantong University and ethically approved by Jiangsu Administration Committee of Experimental Animals. The animals were randomly divided into 11 experimental groups and a sham group (n = 6 each group). For experimental groups, the rats were anesthetized with intraperitoneal injection of mixed narcotics (100 mg/kg ketamine plus 10 mg/kg xylazine) before the sciatic nerve was exposed through an incision at the mid-thigh of the left hind limb for transection to leave a 10-mm long defect. For relieving post-surgery pain, the rats were given buprenorphine (0.2 mg/kg subcutaneously) when performing the operation. The rats were sacrificed by cervical decapitation to harvest TA muscles at 0.25, 0.5, 3, 6, 12, and 24 h, 3, 7, 14, 21, and 28 days post nerve injury respectively. For a sham group, the rats were subjected to similar surgical procedures but without sciatic nerve transection, and TA muscles were harvested immediately after rat sacrifice (referred to as at 0 h post surgery).
The TA muscle sample was homogenized, and RNA was extracted with RNasey Mini Kit (Qiagen, San Francisco, CA, United States) as per kit guidelines. The quantity of total RNA was measured by using NanoDrop (ND-2000 Thermo Scientific, Delaware, ME, United States), and RNA integrity was evaluated by using Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, United States). Microarray analysis consisted of several successive steps, including sample labeling, hybridization, scanning, and image analysis, according to standard protocols. In brief, total RNA was transcribed to double strand cDNA and synthesized into cRNA. After labeled with Cyanine-3-CTP (Cy3), the cRNA was hybridized onto Agilent SurePrint G3 Rat GE microarray (8 × 60,000). Agilent Scanner (G2505C) was used to scan microarray slides, and Agilent Feature Extraction Software (version 10.7.1.1) was used to read the scanned images. The raw data were quantile-normalized and log2 transformed for processing with Genespring Software (version13.1). The probes for which at least 60% samples have flags in “Detected” were chosen for further data analysis.
The expression level of mRNAs at each time was compared to that at 0 h post surgery (control). The gene whose mRNA expression showed a fold change of greater or lower than 2 and a P-value of less than 0.05 was considered a DEG.
Time-Specific Pattern Identification
Microarray data mining was performed by a combination of self-organizing map (SOM), Pearson correlation heatmap, principal component analysis (PCA), and hierarchical clustering analysis to identify time-specific patterns in the mRNA expression of DEGs.
SOMs, as a class of unsupervised neural networks, is used to build up 1-, 2-, or 3-dimensional representation of high dimensional information (Chavez-Alvarez et al., 2014). Here SOMs were created by using Matlab software R2008a (MathWorks, Natick, MA, United States). Pearson correlation coefficients were calculated by the cor function of R “stats’ package1 for ensuring the reliability of expression data, and the results were viewed as a heatmap by ‘pheatmap’ R package (Kolde, 2015). PCA is an unsupervised learning method for constructing high-dimension data patterns without reference to prior knowledge, and the transformation in PCA is achieved by geometrically projecting data onto special dimensions called principal components (PCs) (Lever et al., 2017). Here PCA were performed by R ‘prcomp’ function of R ‘stats’ package1, where the first 3 PCs were selected to view the expression pattern of samples. Hierarchical clustering analysis was performed from Euclidean distance matrix data by using the complete-linkage cluster in the R ‘dendextend’ package (Galili, 2015), and Euclidean distance measure was used to calculate the similarity in gene expressions between samples, and to group samples into clusters by the Ward.D method in the R ‘stats’ package1.
Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Release 91.0) and Gene Ontology (GO) category database were applied for functional annotation of DEGs. KEGG pathway enrichment analysis was tested upon hypergeometric distribution by R ‘phyper’ function (Khatri et al., 2012). Enrichment analysis of GO categories was performed by DAVID tools2. Up- or down- regulated genes were separated for functional enrichment analysis. KEGG or GO categories with a P-value of less than 0.05 were considered the significantly enriched ones, which were viewed in a heatmap. The expression profile for each enriched category was shown based on the average z-score transformed expression value. Z score (x) = (x−μ)/δ, where μ is the mean value and δ is the corresponding standard deviation (Parikh et al., 2010).
Quantitative Real Time PCR (qRT-PCR)
For qRT-PCR validation, the primers, designed by Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, United States), were purchased from Generay Biotech (Shanghai, China). They are listed in Supplementary Table S1. The RNA samples were reverse transcribed and synthesized to cDNA with an iScript cDNA synthesis kit as per the kit manual (Bio-Rad, Hercules, CA, United States). The qRT-PCR reaction was conducted on an ABI7500 Real Time System (Applied Biosystems, Foster City, CA, United States) by using SYBR Green I mix (Roche, Lewes, United Kingdom). The thermal cycle program was as follows: denaturation at 95°C for 10 s, annealing at 60°C for 30 s, extension at 70°C for 10 s (He et al., 2016). All reactions were performed in triplicates, and the specificity of PCR amplification was determined by melting point curve analysis. Quantification of the mRNA expression was performed on a 7500 real-time PCR system (Applied Biosystems). After the cycle threshold (Ct) values were determined based on the cycle number of PCR at which a threshold in fluorescence emission reached, the relative mRNA expression was calculated using comparative 2–Δ Δ Ct method (Livak and Schmittgen, 2001), where beta actin (ACTB) was used as an internal control.
Differential Expression of Genes in Denervated TA Muscles
The RNA samples were extracted from TA muscles and were subjected to cDNA microarray analysis. In each assay, at least three biologically independent samples were used to make data reproducible. A huge number of genes in experimental groups were differentially expressed as compared to those in control (Supplementary Table S2). Out of the DEGs, 14, 282, 327, 255, 310, 641, 1383, 1661, 1891, 1854, and 1578 genes were up-regulated, and 69, 104, 317, 366, 341, 976, 1488, 1951, 2348, 2212, and 1905 genes were down-regulated at 0.25, 0.5, 3, 6, 12, and 24 h, 3, 7, 14, 21, and 28 days post nerve injury, respectively (Figure 1A). Interestingly, the number of DEGs displayed a marked increase at 24 h post nerve injury, suggesting that the atrophic process might be extensively activated at this time point.
Figure 1. Four distinct transcriptional phases during denervation-induced muscle atrophy. (A) The number of DEGs at different times in TA muscles post nerve injury as compared to control (0 h post sham surgery). (B) SOMs: Color-coded output maps representing the final weight of denervated muscle atrophy. The cluster centroid value is coded with a range of colors from blue (for the lowest expression level) to red (for the highest expression level). (C) Pearson correlation heatmap between different times post nerve injury, as labeled on down and right margins. Arrows indicate three nodal transitions. (D) PCA of DEGs at different times post nerve injury. The percentage of variance captured by principal components (PCs) is shown. (E) Dendrogram of gene expression profiles, created by hierarchical clustering analysis. The three nodal transitions in temporal gene expressions segregate the four distinct transcriptional phases within the period of 28 days post nerve injury, and they were designated as “oxidative stress stage,” “inflammation stage,” “atrophy stage” and “atrophic fibrosis stage,” respectively.
Based on SOM analysis, the gene expression profile in TA muscles might be considered to be composed of three main portions: (I) covering control (0 h), 0.25, 0.5, 3, 6 and 12 h; (II) covering 24 h; (III) covering 3, 7, 14, 21, and 28 days (Figure 1B), implying that 24 h post nerve injury was a critical time point in gene expression changes. Pearson correlation heatmap of microarray data indicated that three nodal transitions between 12 and 24 h, 24 h and 3 days, and 7 and 14 days post nerve injury were observable in the expression profile of DEGs (Figure 1C). PCA suggested that the period of 28 days post nerve injury could be partitioned into four distinct transcriptional phases by taking the three nodal transitions as boundaries between phases (Figure 1D). A dendrogram, created by hierarchical cluster (Figure 1E), provided further evidence for the existence of four transcriptional phases.
KEGG Pathway Enrichment Analysis
KEGG pathways related to DEGs were identified by using KEGG database. The enriched categories of KEGG pathways for the up-regulated genes at different times after sciatic nerve transection were labeled in Figure 2. During transcriptional phase I (0.25–12 h post nerve injury), the pathways, such as IL-17 signaling pathway, p53 signaling pathway, PPAR signaling pathway, drug metabolism-cytochrome P450, metabolism of xenobiotics by cytochrome P450, and HIF-1 signaling pathway, were enriched; also enriched were calcium signaling pathway, apelin signaling pathway, circadian entrainment, and so on. The data suggested that oxidative stress was mainly activated at this stage. During transcriptional phase II (24 h post nerve injury), the pathways, such as IL-17 signaling pathway, TNF signaling pathway, JAK-STAT signaling pathway, cytokine-cytokine receptor interaction, TGF-beta signaling pathway, chemokine signaling pathway, and NOD-like receptor signaling pathway, were all enriched, suggesting that inflammation was widely activated at this stage. During transcriptional phase III and IV (3–28 days post nerve injury), the pathways, such as lysosome, Fc gamma R-mediated phagocytosis, phagosome, endocytosis, apoptosis, FoxO signaling pathway, cell adhesion molecules (CAMs), extracellular matrix (ECM)-receptor interaction, cellular senescence, hippo signaling and Wnt signaling pathway, were enriched, suggesting that muscle atrophy and fibrosis were activated at these two stages.
Figure 2. KEGG pathway analysis of up-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of pathways at different times post nerve injury. (B) The average expression profiles of DEGs involved in major pathways. The y-axis represents the averaged z-score of genes in the specific KEGG pathway.
The enriched categories of KEGG pathways for the down-regulated genes at different times after sciatic nerve transection were labeled in Figure 3. During transcriptional phase I, the pathways, such as cell cycle, DNA replication, rap1 signaling pathway, TGF-beta signaling pathway, insulin resistance, MAPK signaling pathway and Ras signaling pathway, were enriched, suggesting that proliferation was inhibited and oxidative stress was activated at this stage. During transcriptional phase II, the pathways, such as MAPK signaling pathway, AMPK signaling pathway, insulin signaling pathway and adipocytokine signaling pathway, were enriched. During transcriptional phases III and IV, the pathways, such as metabolic pathway, oxidative phosphorylation, citrate cycle, glycolysis were enriched; also enriched were MAPK signaling pathway, AMPK signaling pathway, insulin signaling pathway, adipocytokine signaling pathway, Rap1 signaling pathway, PI3K-AKT signaling pathway, gap junction, and so on. These results collectively suggested that energy metabolism was significantly inhibited at the last three stages.
Figure 3. KEGG pathway analysis of down-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of pathways at different times post nerve injury. (B) The average expression profiles of DEGs involved in major pathways. The y-axis represents the averaged z-score of genes in the specific KEGG pathway.
GO Biological Process Analysis
Biological processes related to DEGs were identified using GO analysis. The enriched categories of biological processes for the up-regulated genes at different times after sciatic nerve transection were labeled in Figure 4. During transcriptional phase I, the categories, such as response to calcium ion, response to cytokine, response to hypoxia, and oxidative-reduction process, were enriched; also enriched were MAPK cascade, negative regulation of cell proliferation, and positive regulation of cell death, and so on. These data suggested that stress response progress was activated at this stage. During transcriptional phase II, the categories, such as response to lipopolysaccharide, neutrophil chemotaxis, immune response, positive regulation of inflammatory response, inflammatory response, cellular response to interleukin-1, cellular response to tumor necrosis factor, and lipopolysaccharide-mediated signaling pathway, were enriched, suggesting that inflammatory response was widely activated at this stage. During transcriptional phase III and IV, the categories, such as endocytosis, phagocytosis, extracellular matrix organization, and collagen fibril organization, were enriched, suggesting that biological processes of muscle atrophy and fibrosis were activated at these two stages.
Figure 4. GO biological process analysis of up-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of biological processes at different times post nerve injury. (B) The average expression profiles of DEGs involved in major biological processes. The y-axis represents the averaged z-score of genes in the specific biological process terms.
The enriched categories of biological processes for the down-regulated genes at different times after sciatic nerve transection were labeled in Figure 5. During transcriptional phase I, the categories, such as regulation of cell proliferation, DNA replication, cell cycle, positive regulation of cell proliferation, regulation of cell cycle, cell proliferation and positive regulation of growth, were enriched, suggesting that the cell proliferation was suppressed at this stage. During transcriptional phase II, the categories, such as response to ischemia, cellular response to hypoxia, mitochondrion organization, and ATP metabolic process, were enriched. During transcriptional phases III and IV, the categories, such as cellular respiration, tricarboxylic acid cycle (TCA cycle), fatty acid beta-oxidation, response to oxidative stress, and glycolytic process, were enriched. These results collectively indicated that biological processes of energy metabolism were significantly suppressed at the last three stages.
Figure 5. GO biological process analysis of down-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of biological processes at different times post nerve injury. (B) The average expression profiles of DEGs involved in major biological processes. The y-axis represents the averaged z-score of genes in the specific biological process terms.
GO Molecular Function Analysis
Molecular functions related to DEGs were also examined by GO analysis, and the enriched categories of molecular functions for the up-regulated genes at different times after sciatic nerve transection were labeled in Figure 6. During transcriptional phase I, the categories, such as calcium-induced calcium release activity, monooxygenase activity, N,N-dimethylaniline monooxygenase activity, and oxidoreductase activity, were enriched, suggesting that oxygenase activity was significantly enhanced at this stage. During transcriptional phase II, the categories, such as CCR1 chemokine receptor binding, CCR5 chemokine receptor binding, cytokine activity, chemokine activity and CCR2 chemokine receptor binding, were enriched, suggesting that inflammatory factor activity was significantly increased at this stage. During transcriptional phases III and IV, the categories, such as extracellular matrix structural constituent, cell adhesion molecule binding, heparin sulfate proteoglycan binding, collagen binding, and structural constituent of ribosome, were enriched, suggesting that structural components of ECM were significantly changed at these two stages.
Figure 6. GO molecular function analysis of up-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of molecular functions at different times post nerve injury. (B) The average expression profiles of DEGs involved in major molecular functions. The y-axis represents the averaged z-score of genes in the specific molecular function terms.
The enriched categories of molecular functions for the down-regulated genes at different times after sciatic nerve transection were labeled in Figure 7. During transcriptional phase I, the categories, such as microtubule motor activity, kinase activity, DNA replication origin binding, DNA helicase activity, RNA polymerase ii transcription factor activity, and transcription factor binding, were enriched, suggesting that the activities of DNA replication and transcription were reduced at this stage. During transcriptional phase II, the categories, such as ATP binding, ubiquitin-protein transferase activity, heat shock protein binding, and steroid hormone receptor activity, were enriched. During transcriptional phases III and IV, the categories, such as NADH dehydrogenase activity, electron carrier activity, chaperone binding, and GTPase activator activity, were enriched. These results collectively indicated that the energy metabolism was significantly inhibited at the last three stages.
Figure 7. GO molecular function analysis of down-regulated expressed genes in TA muscles after sciatic nerve transection. (A) The enriched categories (red color) of molecular functions at different times post nerve injury. (B) The average expression profiles of DEGs involved in major molecular functions. The y-axis represents the averaged z-score of genes in the specific molecular function terms.
Four Transcriptional Phases
Based on KEGG and GO analyses, oxidative stress was activated and cell proliferation was inhibited mainly during transcriptional phase I; inflammation was widely activated during transcriptional phase II; muscle atrophy (phagosome, endocytosis, and lysosome) was activated and energy metabolism was significantly suppressed during transcriptional phase III; muscle atrophy and fibrosis (structural components of extracellular matrix) was activated during transcriptional phase IV. In this way, transcriptional phases I-IV might be sequentially referred to as “oxidative stress stage,” “inflammation stage,” “atrophy stage” and “atrophic fibrosis stage” to concisely specify major functions related to DEGs at different time stages.
To validate microarray data, 12 genes, including Ccnb1, Ccnb2, Cdk1, Kif11, Prc1, Kif23, Pttg1, Plk1, Tnfrsf1a, Cebpb, Myd88, and Serpina3, were subjected to qRT-PCR detection for their differential expressions. As to Ccnb1, Ccnb2, CDK1, Kif11, Prc1, Kif23, and Pttg1, qRT-PCR data were quite consistent with microarray results. As to other five genes, qRT-PCR data were roughly consistent with microarray results as far as the change tendency of gene expression was concerned, despite more or less deviations occurring at some time points (Figure 8). In general, however, qRT-PCR and microarray data were correlated to each other in this study.
Figure 8. qRT-PCR validation of expression levels of DEGs in TA muscles after sciatic nerve transection. The curve plots comparing qRT-PCR and microarray data on the time-dependent expression changes of 12 genes.
Various pathological stimuli provoke skeletal muscle atrophy, which, in its nature, results from an imbalance of protein synthesis and protein degradation. The orchestrated cellular events and complex interactions between signaling pathways, mediated by an enormous number of functional biomolecules, are critically involved in the action mechanisms underlying skeletal muscle atrophy (Bonaldo and Sandri, 2013; Cohen et al., 2015; Bilodeau et al., 2016). Setting out from individual genes or proteins has proven to be a disappointing approach to acquiring a global knowledge about the mechanisms of skeletal muscle atrophy. In recent years, therefore, many studies have examined a large number of DEGs in atrophic muscles by microarray analysis (Gomes et al., 2001; Kunkel et al., 2011; Fedorov et al., 2014; Kruger et al., 2015; Jeong et al., 2017; Weng et al., 2018). Despite high-throughput data obtained in these pioneering studies, this study focused on microarray analysis over a broader time window of denervated skeletal muscle atrophy, and applied various clustering methods and bioinformatic tools to mine microarray data with a hope to obtain the detailed mechanistic information on denervated skeletal muscle atrophy.
We identified a huge number of DEGs in TA muscles after sciatic nerve transection. The number of either up-regulated genes or down-regulated genes in TA muscles showed a time-dependent change during the period of 0.25 h to 28 days post nerve injury. Interestingly, the number of DEGs displays a marked increase at 24 h post nerve injury, and SOM analysis also indicated that 24 h post nerve injury was a critical time point of the gene expression pattern in TA muscles, and that the most number of genes began to be activated at this time point, thus initiating the atrophic process.
Pearson correlation heatmap, PCA, and hierarchical cluster analysis were combined to analyze the gene expression profile in denervated TA muscles. The four transcriptional phases were divided within the period of 28 days post nerve injury, and different transcriptional phases might represent the activation of different functional genes, different signaling pathways, different biological processes, and different molecular functions (Sacheck et al., 2007; Glass and Roubenoff, 2010). This division provided a scientific basis for specific therapeutic interventions at different time stages of denervated skeletal muscle atrophy.
KEGG analysis and GO annotation indicated that oxidative stress was activated and cell proliferation was inhibited during transcriptional phase I, inflammation was widely activated during transcriptional phase II, muscle atrophy was activated and energy metabolism was significantly suppressed during transcriptional phase III, muscle atrophy and fibrosis was activated during transcriptional phase IV. Depending on these features, transcriptional phase I, II, III, and IV are, in sequence, designated as “oxidative stress stage,” “inflammation stage,” “atrophy stage,” and “atrophic fibrosis stage.”
At oxidative stress stage, cytochrome P450-related signaling pathways were remarkably activated. As is well known, cytochrome P450 enzymes catalyze the oxidative transformation of toxic metabolites into reactive oxygen species (ROS, such as hydroxyl radical and peroxide and superoxide anions), and mediate the mass production of ROS from estrogen metabolites (Khan et al., 2015; He L. et al., 2017), thereby becoming the marker of oxidative stress as well as the source of cellular ROS (Sato et al., 2011; Tomasi et al., 2018). On the other hand, monooxygenase activity, as a category of molecular function, was significantly enhanced also at oxidative stress stage. Actually, cytochrome P450 belongs to a class of monooxygenase enzymes, and is the main source of ROS due to its nature of terminal oxidases (Bhattacharyya et al., 2014). Our observations implied that a large amount of ROS was generated in response to oxidative stress during transcriptional phase I. Meanwhile, PPAR and HIF-1 signaling pathways were also activated at oxidative stress stage. PPARs (PPAR-α, -β/δ, and -γ), belonging to transcription factors, form a tightly connected triad, and they are involved in regulation of ROS production and degradation (Aleshin and Reiser, 2013; Morris et al., 2018). HIF-1, with a full name of hypoxia-inducible factor-1, plays a key role in genetic regulation of cellular adaptation to hypoxia. The excessive production of ROS under oxidative stress will activate HIF-1 signaling pathway to eliminate excess ROS and to prevent cellular damage (Jalouli et al., 2017; Eyrich et al., 2019). In a word, we found that oxidative stress was present at an early stage of denervated skeletal muscle atrophy. Although the involvement of oxidative stress in muscle atrophy (e.g., due to cancer cachexia or aging) has been reported (Faitg et al., 2017; Zhao et al., 2018), the genetic elucidation of oxidative stress impacts on denervated muscle atrophy, especially at the early stage, was really one of the highlights in this study. It is believable that this highlight may perhaps be extended to establishing an early antioxidant therapy to delay the progress of denervated skeletal muscle atrophy.
At inflammation stage, the up-regulated genes were related to activation of TNF, JAK-STAT, TGF-beta, NOD-like receptor, and NF-kappa B signaling pathways, as well as cytokine-cytokine receptor interaction, thus triggering a persistent inflammatory response (Philpott et al., 2014; Banerjee et al., 2017; He Y. et al., 2017). The inflammatory pathways are involved in different types of skeletal muscle atrophy caused by cachexia, aging, and/or chronic diseases (Perry et al., 2016; Zimmers et al., 2016; Sakurai and Harashima, 2017; Ma et al., 2018). STAT3 and NF-kappa B respond to the same upstream signal and cooperate to promote the expression of pro-cachectic genes, thus contributing to muscle atrophy (Ma et al., 2017). Interventions with JAK-STAT3 related inhibitors in a cancer cachexia model decrease IL-6 level and alleviate skeletal muscle atrophy (Bonetto et al., 2012). TGF-beta 1 mediates muscle atrophy (Zhu et al., 2017) and promotes fibroblast proliferation, collagen synthesis, and muscle fibrosis (Mendias et al., 2012). Moreover, proinflammatory cytokines (e.g., IL-6 and TNF-alpha) affect aging-induced skeletal muscle atrophy by regulating the related signaling pathways (Fan et al., 2016). Besides inflammatory signaling pathways, several biological processes, including immune response, inflammatory response, cellular response to IL-1 and cellular response to TNF, were also activated at inflammation stage. As to the enhanced molecular function at the same stage, they were CCR1 chemokine receptor binding, CCR2 chemokine receptor binding, CCR5 chemokine receptor binding, and cytokine activity. Our observation about the presence of inflammation stage not only provided new evidence for the significance of inflammation in skeletal muscle atrophy, but also inspired a new idea of using anti-inflammatory therapy for treating denervated skeletal muscle atrophy.
At atrophy stage, the signaling pathways related to the up-regulated genes were mainly phagosomes, lysosomes, endocytosis, and P53 signaling pathways, which were essential in proteolysis (Bonetto et al., 2012; Fanzani et al., 2012; Bonaldo and Sandri, 2013). There is overwhelming evidence that activation of autophagy/lysosomal and proteasome signaling pathways promotes protein degradation and causes muscle atrophy (Zhao et al., 2007). Further studies revealed that insulin activation inhibits the ubiquitin proteolytic system and the lysosomal hydrolysis system, promoting skeletal muscle hypertrophy or relieving skeletal muscle atrophy (O’Neill et al., 2016). Our results showed that the insulin signaling pathway was continuously inhibited to aggravate skeletal muscle atrophy. At atrophic fibrosis stage, the up-regulated genes were mainly involved in ribosomes, CAMs and ECM, phagosomes, and endocytosis. It has been reported that denervation for 14 days induces skeletal muscle fibrosis (Huang et al., 2016) and that synthesis of CAMs and ECM aggravates the fibrosis of injured skeletal muscle (Wynn and Ramalingam, 2012). Activation of ribosome-related signaling pathways might promote the synthesis of enzymes involved in the process of collagenation and fibrosis. These data suggested that at atrophic fibrosis stage, both the atrophic process (induced by protein proteolysis) and the fibrosis process (induced ECM overproduction) coexisted in denervated skeletal muscle.
At oxidative stress and inflammation stages, the down-regulated genes were involved in many molecular functions, including cell proliferation, cell cycle, DNA replication initiation, and mitotic sister chromatid segregation, indicating that cell proliferation was significantly inhibited at the two earlier stages of denervated skeletal muscle atrophy. At atrophy phase and atrophic fibrosis stages, the down-regulated genes were involved in many molecular functions, including cellular respiration, tricarboxylic acid cycle, fatty acid beta-oxidation, response to oxidative stress, and glycolytic process, indicating that the energy metabolism was significantly inhibited at the middle and later stages of denervated skeletal muscle atrophy. This phenomenon has also been noted in other models of skeletal muscle atrophy (Ibebunjo et al., 2013; Rao et al., 2018). It follows that the down-regulated genes during denervated skeletal muscle atrophy may be mainly involved in two molecular functions: cell proliferation and energy metabolism.
It should be mentioned that there are some limitations in this study. The rats used in the experiments grew considerably during the period of 0.25 h–28 days post nerve injury, and the response of skeletal muscle to denervation stimulus was accompanied by the body growth. It seems more suitable that the TA muscle samples harvested from the uninjured, contralateral side of animals at distinct times serve as multiple controls, which may be better than a single control (from sham-operated rats and designated as 0 h post surgery). Another drawback of our study was that lncRNA and microRNA were not investigated. Perhaps, the detection of mRNA expressions is not enough to support the whole transcriptome survey for obtaining more comprehensive information.
In terms of four stages occurring within the period of 0.25 h–28 days post nerve injury, we assume that denervated skeletal muscle atrophy represents a cascade of gene expression changes, and it also consists of a sequence of pathological processes with each being triggered by a previous one. In brief, denervation (loss of innervations) causes the disuse of skeletal muscle, the disuse leads to ischemia and hypoxia for generating a large amount of ROS, ROS production initiates inflammation, and uncontrolled inflammation activates the downstream proteolysis process to induce skeletal muscle atrophy (see schematic in Figure 9). Overall, our results in this study may help understand the molecular mechanisms of denervated skeletal muscle atrophy and contribute to developing targeted therapeutic approaches.
Figure 9. Schematic diagram illustrating a sequence of pathological processes in denervation-induced skeletal muscle atrophy. The disuse of skeletal muscle due to loss of innervation leads to ischemia and hypoxia of skeletal muscle, which contributes the production of ROS. This occurs at oxidative stress stage. A large amount of ROS may induce inflammation at inflammation stage. If inflammation is not controlled in time, it will activate the downstream proteolysis process, contributing to muscle atrophy. This occurs at atrophy stage and atrophic fibrosis stage.
Data Availability Statement
The gene expression dataset is available on EBI and the accession number is E-MTAB-8009.
The animal study was reviewed and approved by Nantong University Administration Committee of Experimental Animals, Jiangsu Province, China.
HS and JG designed the study. YS, RZ, LX, ZH, QW, WM, JZ, and MS performed the experiments. YS, ZH, WM, QW, JG, JZ, and LX collected and assembled the data. ZH, JZ, and WM performed the data analysis. FD provided the scientific expertise. HS wrote the manuscript.
This work was supported by the National Key Research and Development Program of China (Grant No. 2017YFA0104703), National Natural Science Foundation of China (Grant Nos. 81671230, 31730031, and 81871554), the 973 Program (Grant Nos. 2014CB542202 and 2014CB542203), a project funded by Jiangsu Provincial Key Medical Center and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), and 2019 Jiangsu College Students Innovation and Entrepreneurship Training Program.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2019.01298/full#supplementary-material
TABLE S1 | The primers used in the qPCR reactions are listed.
TABLE S2 | Differentially expressed genes (DEGs) are identified from experimental groups compared to sham group.
Aleshin, S., and Reiser, G. (2013). Role of the peroxisome proliferator-activated receptors (PPAR)-alpha, beta/delta and gamma triad in regulation of reactive oxygen species signaling in brain. Biol. Chem. 394, 1553–1570. doi: 10.1515/hsz-2013-2215
Aversa, Z., Zhang, X., Fielding, R. A., Lanza, I., and LeBrasseur, N. K. (2019). The clinical impact and biological mechanisms of skeletal muscle aging. Bone 127, 26–36. doi: 10.1016/j.bone.2019.05.021
Banerjee, S., Biehl, A., Gadina, M., Hasni, S., and Schwartz, D. M. (2017). Erratum to: JAK-STAT signaling as a target for inflammatory and autoimmune diseases: current and future prospects. Drugs 77, 521–546. doi: 10.1007/s40265-017-0772-777
Bhattacharyya, S., Sinha, K., and Sil, P. C. (2014). Cytochrome P450s: mechanisms and biological implications in drug metabolism and its interaction with oxidative stress. Curr. Drug Metab. 15, 719–742. doi: 10.2174/1389200215666141125121659
Bilodeau, P. A., Coyne, E. S., and Wing, S. S. (2016). The ubiquitin proteasome system in atrophying skeletal muscle: roles and regulation. Am. J. Physiol. Cell Physiol. 311, C392–C403. doi: 10.1152/ajpcell.00125.2016
Bodine, S. C., and Baehr, L. M. (2014). Skeletal muscle atrophy and the E3 ubiquitin ligases MuRF1 and MAFbx/atrogin-1. Am. J. Physiol. Endocrinol. Metab. 307, E469–E484. doi: 10.1152/ajpendo.00204.2014
Bodine, S. C., Latres, E., Baumhueter, S., Lai, V. K., Nunez, L., Clarke, B. A., et al. (2001). Identification of ubiquitin ligases required for skeletal muscle atrophy. Science 294, 1704–1708. doi: 10.1126/science.1065874
Bonetto, A., Aydogdu, T., Jin, X., Zhang, Z., Zhan, R., Puzis, L., et al. (2012). JAK/STAT3 pathway inhibition blocks skeletal muscle wasting downstream of IL-6 and in experimental cancer cachexia. Am. J. Physiol. Endocrinol. Metab. 303, E410–E421. doi: 10.1152/ajpendo.00039.2012
Chavez-Alvarez, R., Chavoya, A., and Mendez-Vazquez, A. (2014). Discovery of possible gene relationships through the application of self-organizing maps to DNA microarray databases. PLoS One 9:e93233. doi: 10.1371/journal.pone.0093233
Choi, W., Lee, J., Lee, J., Ko, K. R., and Kim, S. (2018). Hepatocyte growth factor regulates the miR-206-HDAC4 cascade to control neurogenic muscle atrophy following surgical denervation in mice. Mol. Ther. Nucleic Acids 12, 568–577. doi: 10.1016/j.omtn.2018.06.013
Eyrich, N. W., Potts, C. R., Robinson, M. H., Maximov, V., and Kenney, A. M. (2019). Reactive oxygen species signaling promotes hypoxia-inducible factor 1α stabilization in sonic hedgehog-driven cerebellar progenitor cell proliferation. Mol. Cell. Biol. 39:e00268-18.
Fanzani, A., Conraads, V. M., Penna, F., and Martinet, W. (2012). Molecular and cellular mechanisms of skeletal muscle atrophy: an update. J. Cachexia Sarcopenia Muscle 3, 163–179. doi: 10.1007/s13539-012-0074-6
Fedorov, V. B., Goropashnaya, A. V., Stewart, N. C., Toien, O., Chang, C., Wang, H., et al. (2014). Comparative functional genomics of adaptation to muscular disuse in hibernating mammals. Mol. Ecol. 23, 5524–5537. doi: 10.1111/mec.12963
Gomes, M. D., Lecker, S. H., Jagoe, R. T., Navon, A., and Goldberg, A. L. (2001). Atrogin-1, a muscle-specific F-box protein highly expressed during muscle atrophy. Proc. Natl. Acad. Sci. U.S.A. 98, 14440–14445. doi: 10.1073/pnas.251541198251541198
He, L., He, T., Farrar, S., Ji, L., Liu, T., and Ma, X. (2017). Antioxidants maintain cellular redox homeostasis by elimination of reactive oxygen species. Cell Physiol. Biochem. 44, 532–553. doi: 10.1159/000485089
He, Y., Xiao, Y., Yang, X., Li, Y., Wang, B., Yao, F., et al. (2017). SIRT6 inhibits TNF-alpha-induced inflammation of vascular adventitial fibroblasts through ROS and Akt signaling pathway. Exp. Cell Res. 357, 88–97. doi: 10.1016/j.yexcr.2017.05.001
He, Q., Qiu, J., Dai, M., Fang, Q., Sun, X., Gong, Y., et al. (2016). MicroRNA-351 inhibits denervation-induced muscle atrophy by targeting TRAF6. Exp. Ther. Med. 12, 4029–4034. doi: 10.3892/etm.2016.3856
Huang, Q. K., Qiao, H. Y., Fu, M. H., Li, G., Li, W. B., Chen, Z., et al. (2016). Mir-206 attenuates denervation-induced skeletal muscle atrophy in rats through regulation of satellite cell differentiation via tgf-beta1, smad3, and hdac4 signaling. Med. Sci. Monit. 22, 1161–1170. doi: 10.12659/msm.897909
Ibebunjo, C., Chick, J. M., Kendall, T., Eash, J. K., Li, C., Zhang, Y., et al. (2013). Genomic and proteomic profiling reveals reduced mitochondrial function and disruption of the neuromuscular junction driving rat sarcopenia. Mol. Cell. Biol. 33, 194–212. doi: 10.1128/MCB.01036-1012
Jalouli, M., Mokas, S., Turgeon, C. A., Lamalice, L., and Richard, D. E. (2017). Selective HIF-1 regulation under nonhypoxic conditions by the p42/p44 MAP kinase inhibitor PD184161. Mol. Pharmacol. 92, 510–518. doi: 10.1124/mol.117.108654
Jeong, H. O., Park, D., Im, E., Lee, J., Im, D. S., and Chung, H. Y. (2017). Determination of the mechanisms that cause sarcopenia through cDNA microarray. J. Frailty Aging 6, 97–102. doi: 10.14283/jfa.2017.13
Khan, M. A., Zubair, H., Srivastava, S. K., Singh, S., and Singh, A. P. (2015). Insights into the role of microRNAs in pancreatic cancer pathogenesis: potential for diagnosis, prognosis, and therapy. Adv. Exp. Med. Biol. 889, 71–87. doi: 10.1007/978-3-319-23730-5_5
Kruger, K., Dischereit, G., Seimetz, M., Wilhelm, J., Weissmann, N., and Mooren, F. C. (2015). Time course of cigarette smoke-induced changes of systemic inflammation and muscle structure. Am. J. Physiol. Lung Cell Mol. Physiol. 309, L119–L128. doi: 10.1152/ajplung.00074.2015
Kunkel, S. D., Suneja, M., Ebert, S. M., Bongers, K. S., Fox, D. K., Malmberg, S. E., et al. (2011). mRNA expression signatures of human skeletal muscle atrophy identify a natural compound that increases muscle mass. Cell Metab. 13, 627–638. doi: 10.1016/j.cmet.2011.03.020
Lala-Tabbert, N., Lejmi-Mrad, R., Timusk, K., Fukano, M., Holbrook, J., St-Jean, M., et al. (2019). Targeted ablation of the cellular inhibitor of apoptosis 1 (cIAP1) attenuates denervation-induced skeletal muscle atrophy. Skelet Muscle 9:13. doi: 10.1186/s13395-019-0201-206
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Ma, J. F., Sanchez, B. J., Hall, D. T., Tremblay, A. K., Di Marco, S., and Gallouzi, I. E. (2017). STAT3 promotes IFNgamma/TNFalpha-induced muscle wasting in an NF-kappaB-dependent and IL-6-independent manner. EMBO Mol. Med. 9, 622–637. doi: 10.15252/emmm.201607052
Madaro, L., Passafaro, M., Sala, D., Etxaniz, U., Lugarini, F., Proietti, D., et al. (2018). Denervation-activated STAT3-IL-6 signalling in fibro-adipogenic progenitors promotes myofibres atrophy and fibrosis. Nat. Cell Biol. 20, 917–927. doi: 10.1038/s41556-018-0151-y
Marra, A. M., Arcopinto, M., Bossone, E., Ehlken, N., Cittadini, A., and Grunig, E. (2015). Pulmonary arterial hypertension-related myopathy: an overview of current data and future perspectives. Nutr. Metab. Cardiovasc. Dis. 25, 131–139. doi: 10.1016/j.numecd.2014.10.005
Mendias, C. L., Gumucio, J. P., Davis, M. E., Bromley, C. W., Davis, C. S., and Brooks, S. V. (2012). Transforming growth factor-beta induces skeletal muscle atrophy and fibrosis through the induction of atrogin-1 and scleraxis. Muscle Nerve 45, 55–59. doi: 10.1002/mus.22232
Morris, C. J., Kameny, R. J., Boehme, J., Gong, W., He, Y., Zhu, T., et al. (2018). KLF2-mediated disruption of PPAR-gamma signaling in lymphatic endothelial cells exposed to chronically increased pulmonary lymph flow. Am. J. Physiol. Heart Circ. Physiol. 315, H173–H181. doi: 10.1152/ajpheart.00635.2017
O’Neill, B. T., Lee, K. Y., Klaus, K., Softic, S., Krumpoch, M. T., Fentz, J., et al. (2016). Insulin and IGF-1 receptors regulate FoxO-mediated signaling in muscle proteostasis. J. Clin. Invest. 126, 3433–3446. doi: 10.1172/jci86522
Oost, L. J., Kustermann, M., Armani, A., Blaauw, B., and Romanello, V. (2019). Fibroblast growth factor 21 controls mitophagy and muscle mass. J Cachexia Sarcopenia Muscle 10, 630–642. doi: 10.1002/jcsm.12409
Parikh, J. R., Klinger, B., Xia, Y., Marto, J. A., and Bluthgen, N. (2010). Discovering causal signaling pathways through gene-expression patterns. Nucleic Acids Res. 38, W109–W117. doi: 10.1093/nar/gkq424
Perry, B. D., Caldow, M. K., Brennan-Speranza, T. C., Sbaraglia, M., Jerums, G., Garnham, A., et al. (2016). Muscle atrophy in patients with Type 2 Diabetes Mellitus: roles of inflammatory pathways, physical activity and exercise. Exerc. Immunol. Rev. 22, 94–109.
Philpott, D. J., Sorbara, M. T., Robertson, S. J., Croitoru, K., and Girardin, S. E. (2014). NOD proteins: regulators of inflammation in health and disease. Nat. Rev. Immunol. 14, 9–23. doi: 10.1038/nri3565
Pigna, E., Renzini, A., Greco, E., Simonazzi, E., Fulle, S., Mancinelli, R., et al. (2018). HDAC4 preserves skeletal muscle structure following long-term denervation by mediating distinct cellular responses. Skelet Muscle 8:6. doi: 10.1186/s13395-018-0153-152
Qiu, J., Fang, Q., Xu, T., Wu, C., Xu, L., Wang, L., et al. (2018). Mechanistic role of reactive oxygen species and therapeutic potential of antioxidants in denervation- or fasting-induced skeletal muscle atrophy. Front. Physiol. 9:215. doi: 10.3389/fphys.2018.00215
Sacheck, J. M., Hyatt, J. P., Raffaello, A., Jagoe, R. T., Roy, R. R., Edgerton, V. R., et al. (2007). Rapid disuse and denervation atrophy involve transcriptional changes similar to those of muscle wasting during systemic diseases. FASEB J. 21, 140–155. doi: 10.1096/fj.06-6604com
Tang, H., Inoki, K., Lee, M., Wright, E., Khuong, A., Khuong, A., et al. (2014). mTORC1 promotes denervation-induced muscle atrophy through a mechanism involving the activation of FoxO and E3 ubiquitin ligases. Sci. Signal. 7:ra18. doi: 10.1126/scisignal.2004809
Tomasi, M. L., Ramani, K., Ryoo, M., Cossu, C., Floris, A., Murray, B. J., et al. (2018). SUMOylation regulates cytochrome P450 2E1 expression and activity in alcoholic liver disease. FASEB J. 32, 3278–3288. doi: 10.1096/fj.201701124R
Weng, J., Zhang, P., Yin, X., and Jiang, B. (2018). The whole transcriptome involved in denervated muscle atrophy following peripheral nerve injury. Front. Mol. Neurosci. 11:69. doi: 10.3389/fnmol.2018.00069
Wu, C. L., Cornwell, E. W., Jackman, R. W., and Kandarian, S. C. (2014). NF-kappaB but not FoxO sites in the MuRF1 promoter are required for transcriptional activation in disuse muscle atrophy. Am. J. Physiol. Cell Physiol. 306, C762–C767. doi: 10.1152/ajpcell.00361.2013
Zhao, J., Brault, J. J., Schild, A., Cao, P., Sandri, M., Schiaffino, S., et al. (2007). FoxO3 coordinately activates protein degradation by the autophagic/lysosomal and proteasomal pathways in atrophying muscle cells. Cell Metab. 6, 472–483. doi: 10.1016/j.cmet.2007.11.004
Zhao, J., Tian, Z., Kadomatsu, T., Xie, P., Miyata, K., Sugizaki, T., et al. (2018). Age-dependent increase in angiopoietin-like protein 2 accelerates skeletal muscle loss in mice. J. Biol. Chem. 293, 1596–1609. doi: 10.1074/jbc.M117.814996
Zhu, X., Kny, M., Schmidt, F., Hahn, A., Wollersheim, T., Kleber, C., et al. (2017). Secreted frizzled-related protein 2 and inflammation-induced skeletal muscle atrophy. Crit. Care Med. 45, e169–e183. doi: 10.1097/CCM.0000000000002056
Keywords: denervation, skeletal muscle atrophy, microarray, bioinformatics, oxidative stress, inflammation
Citation: Shen Y, Zhang R, Xu L, Wan Q, Zhu J, Gu J, Huang Z, Ma W, Shen M, Ding F and Sun H (2019) Microarray Analysis of Gene Expression Provides New Insights Into Denervation-Induced Skeletal Muscle Atrophy. Front. Physiol. 10:1298. doi: 10.3389/fphys.2019.01298
Received: 27 July 2019; Accepted: 27 September 2019;
Published: 11 October 2019.
Edited by:Jose Renato Pinto, Florida State University, United States
Reviewed by:Christopher P. Cardozo, United States Department of Veterans Affairs, United States
Sven Tågerud, Linnaeus University, Sweden
Copyright © 2019 Shen, Zhang, Xu, Wan, Zhu, Gu, Huang, Ma, Shen, Ding and Sun. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
†These authors have contributed equally to this work