Molecular Networks of Human Muscle Adaptation to Exercise and Age

Summary

Physical activity and molecular ageing presumably interact to precipitate musculoskeletal decline in humans with age. Herein, we have delineated molecular networks for these two major components of sarcopenic risk using multiple independent clinical cohorts. We generated genome-wide transcript profiles from individuals (n = 44) who then undertook 20 weeks of supervised resistance-exercise training (RET). Expectedly, our subjects exhibited a marked range of hypertrophic responses (3% to +28%), and when applying Ingenuity Pathway Analysis (IPA) up-stream analysis to ∼580 genes that co-varied with gain in lean mass, we identified rapamycin (mTOR) signaling associating with growth (P = 1.4×10⁻³⁰). Paradoxically, those displaying most hypertrophy exhibited an inhibited mTOR activation signature, including the striking down-regulation of 70 rRNAs. Differential analysis found networks mimicking developmental processes (activated all-trans-retinoic acid (ATRA, Z-score = 4.5; P = 6×10⁻¹³) and inhibited aryl-hydrocarbon receptor signaling (AhR, Z-score = −2.3; P = 3×10⁻⁷)) with RET. Intriguingly, as ATRA and AhR gene-sets were also a feature of endurance exercise training (EET), they appear to represent “generic” physical activity responsive gene-networks. For age, we found that differential gene-expression methods do not produce consistent molecular differences between young versus old individuals. Instead, utilizing two independent cohorts (n = 45 and n = 52), with a continuum of subject ages (18–78 y), the first reproducible set of age-related transcripts in human muscle was identified. This analysis identified ∼500 genes highly enriched in post-transcriptional processes (P = 1×10⁻⁶) and with negligible links to the aforementioned generic exercise regulated gene-sets and some overlap with ribosomal genes. The RNA signatures from multiple compounds all targeting serotonin, DNA topoisomerase antagonism, and RXR activation were significantly related to the muscle age-related genes. Finally, a number of specific chromosomal loci, including 1q12 and 13q21, contributed by more than chance to the age-related gene list (P = 0.01–0.005), implying possible epigenetic events. We conclude that human muscle age-related molecular processes appear distinct from the processes regulated by those of physical activity.

Introduction

Discovery of the biological determinants of muscle mass and functional molecular phenotypes has substantial bearing on the fields of human performance (e.g. hypertrophy, strength or endurance adaptations [1], [2]) and human health (countering muscle atrophy and deconditioning occurring in older age or with conditions such as cancer cachexia [3], [4], respiratory disease or medically enforced immobilization (e.g. hospitalized bed-rest, cast-immobilization [5], [6])). Resistance exercise (RE) training (RET) is an effective intervention to increase muscle mass in many, but not all people, and thus provides an excellent opportunity to study gene-network regulation during muscle hypertrophy and the proposed relationship to muscle aging. Many exogenous factors may influence RET-induced hypertrophy including manipulation of exercise volume [7], intensity [8] and adequate macronutrient availability [9] all of which interact with an individual's genotype to determine muscle growth. Establishing the molecular diagnostics that enable a personalized approach to tackle ageing has great appeal.

To date, the molecular regulation of muscle hypertrophy has mostly focused on aspects of post-genomic signaling, with early work concluding that canonical insulin-like growth factor (IGF-1) signals control muscle hypertrophy though a phosphatidyl-inositol-3 kinase/protein kinase B/mechanistic target of rapamycin (PI3K/AKT/mTORc1) pathway [10], abbreviated to mTOR. This protein complex can control cell growth through two mechanisms; firstly, mTORC1 regulates efficiency of translation through inducing phosphorylation of its substrates, ribosomal protein S6 kinase (S6K1) and 4E binding protein 1 (4EBP1) and, secondly, mTORC1 increases translational capacity through regulating ribosomal RNA (rRNA) production within the nucleolus [11]. There are however conflicting data regarding the importance of mTOR regulation (protein phosphorylation or target gene mRNA responses) or its up-stream regulators, and acute anabolic or chronic growth responses to resistance exercise [12]–[17] reported from the same laboratories, indicating that important biological rather than methodological issues remain to be identified.

More recent evidence indicates that the mechanisms regulating muscle hypertrophy go beyond the canonical IGF-1/PI3K/AKT/mTORc1 pathway. While circulating IGF-1 concentrations do not determine RET-induced hypertrophy in humans [18], hypertrophy has now been shown to potentially occur through both PI3K-AKT [19] and mTOR [20] independent pathways, even in pre-clinical models. Perhaps the most convincing observation in favor of a more divergent regulation of muscle growth, is the fact that disparate exercise modes (e.g. RET vs. endurance exercise training (EET)) can produce similar protein signaling patterns in humans [21]. This suggests that the molecules, so far studied, are pleiotropic and in our opinion probably important for any type of tissue remodeling, regardless of the final physiological phenotype [22]. Therefore, a more innovative approach is needed to define links between molecules and ensuing in vivo physiological adaptations, than can be achieved with targeted western-based molecular profiling.

Exercise training has also been postulated as a key tool to reverse the impact of ageing on human skeletal muscle phenotypes [23], [24]. Yet, while some ‘genomic features’ of ageing have been reported [25], we have noticed that the available global molecular profiles of human muscle [23], [24], [26], [27] do not identify consistent molecular features. Furthermore, our recent work has highlighted that physiological adaptations to exercise, whether that be hypertrophy [28] or aerobic function [29], are highly heterogeneous in humans, implying that exercise may not be able to “reverse” muscle ageing [23] for some people. For example, following >10-wk of supervised EET, ∼20% of subjects show no improvement in aerobic capacity while ∼30% demonstrate no improvement in insulin sensitivity [30], [31]. Similarly, we reported muscle hypertrophy ranging from 0.8 to 6.0 kg [28], while Raue et al reported changes in muscle cross-sectional area (CSA) ranging −1.2 to +10.4 cm² [24]. Both of these RET studies reported that the outcome of supervised progressive RET did not relate to pre-existing differences in characteristics (i.e., gender, age, pre-existing muscle mass, physical activity levels or dietary habits) indicating that there is not a simple explanation for the heterogeneity of the gains in lean mass [28], [32].

In recent years, we have focused on using the heterogeneous responses to exercise training and OMIC techniques to uncover molecular networks regulated by EET [33][29] or generate molecular predictors of trainability [34], directly in humans. The aim of the present work is to produce the first reproducible molecular signature of human muscle age, and examine how such a profile relates to new and established exercise adaptation gene networks. We generated new gene-chip profiles from muscle samples derived from two independent clinical cohorts, with a continuous range of ages (18–79 y) and which originate from distinct environments (UK and USA) and which were independently processed in the laboratory. We also generated a new data set of paired global RNA-responses to a supervised 20-wk RET program (N = 44), as well as utilizing various sets of published acute-RET and chronic-EET gene-chips (total N = 200) data sets. Finally, Ingenuity's new IPA up-stream analysis tool [35] was used to identify key features, within these novel age and exercise signatures, to provide independent and robust molecular insight into the heterogeneous nature of muscle hypertrophy, and human muscle age.

Results

The paired differential analysis, comparing expression in 38 from 44 subjects before and after RET yielded a dataset of ∼700 regulated genes (Dataset S1) and this gene-list related to a few upstream regulators in IPA (Dataset S2). For logical reasons we used only the 38 subjects that demonstrated a training effect [31]. This list included a 62-gene network (Figure 1A), representing the transcriptional action of all-trans-retinoic acid (Tretinoin) and this signature was highly ‘activated’ following 20-wk RET (Z-score = 4.5 for directional consistency; P-value for transcript overlap (p = 6×10⁻¹³)). In contrast, aryl hydrocarbon (AhR) and V-Ki-ras2 Kirsten rat sarcoma viral oncogene homolog (KRAS) signaling (Figure S1A) was inhibited, such that repression of connected genes was removed, and thus target mRNAs up-regulated.

We utilised the differential RNA response to 6-wks of EET [29] to address the question of specificity of this RET molecular profile to conditions of muscle growth. We have previously speculated that a core set of gene-networks will be common to all types of exercise and could represent basic determinants of tissue adaptability [22]. Indeed, we found that several significant upstream transcriptional regulators from the RET profile were also significant significantly regulated by EET (Dataset S1) and the majority of individual genes in the RET RNA signatures were activated in an identical manner between RET and EET (Figure 1B and Figure S1B). Thus while differential expression in the Derby Resistance Exercise Training (DRET) study clearly identifies a number of tissue remodeling related processes, these are not specific to exercise training modality.

We then utilized the dataset from the Trappe laboratory (GSE28422, ‘Study A’ from GEO). We removed the following samples (or the related paired sample) because the following samples appeared to have technical issues; GSM702424, GSM702439, GSM702431, GSM702421, GSM702414, GSM702436 and GSM702405 GSM702393, GSM702399, GSM702411. Unfortunately we were declined access to the phenotype data (individual muscle growth changes) for each subject (Dr S Trappe, personal communication) and thus we were unable to re-analyze their data to capture the growth related responses, as their data-set included a number of ‘non-responders’ for lean mass gains. Nonetheless, we used the acute RE response dataset found in the Trappe study (GSE28422) so that we could explore the overlap between the acute RNA responses and chronic RET adaptation. What we noted was that <2% of the RNA changes that occur with chronic adaptation in the DRET study, were regulated in response to acute resistance exercise.

Molecular networks associated with muscle hypertrophy

Of greater physiological importance was our effort to identify genes which link to the magnitude of muscle hypertrophy in humans. Quantitative SAM analysis [3], [36] identified 642 probe-sets and interestingly the majority of genes identified were negatively correlated with gains in lean mass (Dataset S2). The probe-set list was imported into IPA (filtered at a 5% FDR) and 384 genes could be mapped to the IPA database for pathway analysis. The enrichment score generated for the EIF2 canonical pathway enrichment was extremely significant (P-value = 1×10⁻⁶⁴) and the combination of these observations indicates that the gene-list was both predictable (valid with extensive literature) and statistically very robust (Figure S2). We identified a number of regulators that could be responsible for regulating the transcriptional signature that correlated to gain in lean mass (Dataset S2).

Figure 2A presents the most striking finding, an active rapamycin signature, equating to inhibition of mTOR signaling [37]. This signature was comprised of genes that almost entirely negatively correlated with lean-mass gain (Z-score = 2.8 for directional consistency; P-value for transcript overlap p = 1.4×10⁻³⁰). In short, subjects that demonstrated the largest gains in lean tissue mass following 20-wks of RET have suppressed mTOR signaling over the training period, a novel and controversial observation (all raw data was manually checked for consistency of direction). A second major transcriptional regulator was MYC, and the gene-list driving its inclusion (Z-score −5.8 for directional consistency; p-value for transcript overlap p = 4×10⁻⁴³) overlapped with the rapamycin network, equating to the inhibition of MYC activity. MYC can be an upstream activator of mTOR signaling in cell culture systems [38]. Thus these two robust observations are consistent, and notably the signature evidence is based on entirely independent data.

Close examination of the rapamycin associated gene network (Figure 2A) reveals ∼30 ribosomal RNA genes (a total of 70 genes were in the lean mass gain associated list but not all featured in the data-base of rapamycin regulated genes). The remaining genes (Figure 2A) belonged to metabolic processes and other facets of protein metabolism or signaling. To more easily appreciate the characteristics of those subjects that were found to have the greatest increase in lean tissue mass combined with a reduction in ribosomal gene abundance, we plotted the quartile response in lean mass (Figure 2B). Baseline lean mass could not explain our observation and in fact the four groups had the same mean age, physiological characteristics, while the highest and lowest quartiles for lean mass gain had exactly the same proportion of males and females (Table 1). Figure 2C specifically illustrates a subset of the ribosomal RNA genes, however all other rRNA genes were consistent with this plot. There was an almost universal pattern of down-regulation observed in high responders (shown in purple) while subjects that were unable to increase muscle mass (shown in blue) substantially up-regulated these genes, as would be expected from the qSAM statistical analysis (FDR<5%) in Figure 1A.

Thus, our analysis strategy enabled the discovery of an entirely novel in vivo feature of the mTOR growth pathway, while standard differential RNA expression analysis (pre vs. post sample) could not. We also plotted the relationship between physiological characteristics, protein-phosphorylation during acute net anabolic situations (resistance exercise coupled with feeding) with lean mass gain in these subjects using principal component analysis (PCA). Selected variables were scaled and principal components 1 and 2 are presented. In Figure 3A, it is abundantly clear that none of the metabolic or physiological characteristics shared variance with the main component capturing lean mass gain variation following RET. Likewise, while protein kinase abundance or protein kinase activation status varied in a manner consistent with the literature, none of these acute net anabolic responses were correlated to gains in RET induced lean mass or shared variance with gains in lean mass when studied prior to 20-wk RET (Figure 3B). In short, only the unbiased transcriptomics method was able to identify a biological profile distinguishing between high and low responders for lean mass gain.

Age muscle gene networks

Identification of the determinants of skeletal muscle mass has obvious implications for treating age-related sarcopenia. There is no longitudinal molecular analysis of ageing muscle in humans. However using cross-sectional gene-chip data-sets, attempts have been made to identify age-related gene expression changes. For example, Melov et al., reported that differences in gene expression between a cohort of young and old subjects can be removed through a RET program [23]. RET removed some aspect of the ‘inactivity’ related component of the difference between young and old subjects but when we contrasted their age-signature with other publically available muscle age datasets [23], [24], [26], [27] no overlap was apparent.

To investigate this issue further, we utilised the Melov et al., data and the data from Trappe and our lab. We used SAM analysis to compare young with old subjects in each study. For the Trappe study we had 13 young (20–30 y) and 11 old subjects (>80 y). We used baseline samples from DRET, 10 young (20–29 y) versus 16 old subjects (64–75 y), and finally we re-analysed the Melov et al., data [23] which had 26 younger (18–28 y) and 25 older subjects (65–85 y). As can be clearly observed in Figure 4 there was no overlap between the three studies indicating that a reproducible set of ‘ageing’ genes cannot be generated with this statistical or experimental approach. Our re-analysis of the GO analysis (using DAVID) of earlier studies [26], [39], using the appropriate back-ground files [40], also confirmed that mitochondrial RNA changes in ageing cannot be claimed as being a reproducible hallmark of ageing, despite the presumed association with inactivity. Re-analysis of the Melov et al., data did identify a mitochondrial gene expression signature (less significant than the original analysis due to comparison with a more appropriate ontological background) but in that study the older subjects were substantially de-trained and aerobic function was not presented. This gene-set is also known to reflect physical activity [29] and inactivity in humans [41] and thus it shouldn't be attributed to age per se anyway. To consolidate the conclusion that there was no common finding across these three studies, we entered the individual gene-lists in a gene ontology analysis to evaluate if some common pathways were enriched in each list, even if the member genes differed. Only 1 ontological group was common to 2 from the 3-gene lists and it related to mitochondrial processes indicating that even when the older subjects have a lower physical capacity a decline in mitochondrial genes is not always a prominent feature of age-related changes.

Therefore, an alternative approach to identify age-related gene expression profiles in human muscle was required. To achieve this we utilised QSAM, which we have previously applied and validated to some extent for human studies [3]. We applied QSAM to generate a list of muscle transcript levels, which correlated with subject age (20–75 yr), with correction for multiple testing. This allowed us to identify genes that either negatively or positively associated with the subject age. This has not been attempted before because previous studies did not have a sufficiently wide and continuous range of ages to generate such data [23], [24], [26], [39]. However, such an analysis would be of limited value if some of the observations could not be independently reproduced, using a distinct set of clinical samples. To this end, we generated 52 new U133+2 profiles (17–63 yr) from muscle biopsy samples from the HERITAGE Family Study [42] and were able to identify a set of 580 genes or transcripts (Dataset 3) which were correlated with age in both clinical studies and which was enriched in post-transcriptional and chronic disease traits (Figure S3) but not mitochondrial related gene-sets.

We found in IPA that the age-related dataset was consistent with the activation of the PGR (z-score = 2.6 and P-value = 0.001) and RXR (z-score = 2.0 and p-value = 0.0001) proteins and 5-fluorouracil agonism (Z-score = 2.2 and P-value = 0.0005, Figure 5A). Each mediator orchestrated a set of either positive (yellow) or negatively (green) age-correlated genes such that both overlap and direction were similar to the literature-constructed networks. Critically, these networks were not significantly related to the lean-mass associated gene-list (Figure 5B) or differentially regulated by either acute RE (Figure S4A) or chronic endurance exercise (Figure S4B). Thus it is unlikely that these new age-related observations reflect simple confounding factors, such as exercise training or being physically active.

There was also inhibition of certain protein mediators with age (Figure 6A) including c-MYC (z-score = −2.8 and p-value<0.0001) and CLDN7 (z-score = −2.6 and p-value = 0.05). Again, no clear relationship with acute exercise or endurance training was apparent (Figure S5), while a closer association with genes related to gains in lean mass was noted (Figure 6B) with some key exceptions. Notably inhibition of MYC was predicted in both the lean-mass and age-related gene lists (with gene-correlations in the same direction) which we would not expect if muscle ageing was simply the ‘opposite’ of muscle growth or lack of response to exercise training. Furthermore, large differences in gene expression still existed when comparing the age groups and the pre and post-training samples in the Trappe dataset (data not shown). The age-related expression signature was also related to RNA signatures in the Broad Connectivity database, including multiple serotonin antagonists and appears opposite to DNA topoisomerase inhibition (Dataset S4). Finally, we examined whether the age-related genes were over represented at genomic loci using Positional enrichment analysis [43]. Both Chromosome 1 (q12) and 13 (q21.33) had significant hits and the genes associated with those locations can be found in Figure 7A and 7B and the remaining analysis in Dataset S5.

Discussion

We have carried out unbiased molecular analyses on both new and pre-existing human muscle data-sets (acute and chronic exercise, RET or EET and ageing) from which we have been able to define: (i) a ‘generic’ set of molecular networks that are activated irrespective of exercise-training mode in humans, i.e., RET or EET, (ii) the differential effects of acute exercise versus chronic exercise training on molecular networks, (iii) the molecular networks that are specifically modulated in relation to the extent of human muscle hypertrophy, and (iv) the first reproducible set of age-related transcriptional changes further supporting our view that large human sample sizes, using unbiased molecular profiling techniques [4], [29], [44] is an important strategy for translational medical science.

Exercise mode and the specificity of the molecular response

Muscle hypertrophy is the most recognized adaptation to RET. However, there are numerous other adaptations that occur, to support the biochemical, physical and metabolic requirements of a growing muscle. For example, hypertrophy is associated with activation of muscle satellite cells to support growth [45] while RET stimulates angiogenesis proportional to muscle fiber growth (rather than increasing capillary-to-muscle fiber area as EET [46]). Furthermore, along with accretion of contractile proteins, the extracellular matrix (ECM) undergoes substantial remodeling after RET [47]. Clearly then, successful hypertrophy is the summation of complex intra/extra-muscular cross talk to co-ordinate hypertrophy. As such, we believe that teasing out the molecular networks regulating adaptation in vivo requires charting the relationship between the modulation of molecular factors with that of the physiological outcome(s) [29].

Our initial analysis revealed that activation of a Tretionin (all-trans retinoic acid (ATRA)) and inactivation of aryl hydrocarbon receptor (AhR) are common molecular responses to training, irrespective of exercise mode i.e., RET or EET. ATRA is the active form of vitamin A, which serves as a ligand for two families of widely expressed nuclear receptors; the retinoic acid receptors (RAR) that bind to ATRA and the retinoid X receptors (RXR) that bind to its stereoisomer, alitretionin (9-cis-RA). Although little work exists on ATRA signalling in skeletal muscle, Halevy et al. showed exogenous ATRA promoted myogenic cell differentiation [48] which is allied to the function of ATRA amongst various cell types [49]. Given the post-mitotic properties of myonuclei, this may point to a novel link between ATRA-like signalling and aspects of in vivo satellite cell activities induced by exercise training. In support of this idea, we also observed IGF-1 and IGF-2 up-regulation in both RET and EET, and this is thought to play a role in satellite cell activation and differentiation [50], [51][52]. ATRA-like signalling has also been shown to modulate endothelial cell maturation and angiogenesis in tube formation assays [53] suggesting a role for activation of this network in exercise-induced angiogenesis, while angiogenic factors are also associated with satellite cell activation and differentiation [54]. Indeed, we identified up-regulation of homeobox (HOX) genes (e.g. HOXB3/7), which have roles in vascular remodelling and angiogenesis [55]. Thus this collection of genes is likely contributing to vascular remodelling-induced both by RET [56] and EET [57]. In addition, the turnover of ECM components is activated when smooth muscle cells are exposed to exogenous ATRA, thereby suggesting this pathway is involved in ECM remodelling [58]. Indeed, the (activated) ATRA gene list was highly enriched in collagen genes. Finally, while one cannot rule out that the degree of modulation of ECM (which differed between RET and EET) may influence aerobic adaptation [29], [31] it would appear to be a constituent feature of muscle growth and remodelling per se.

The AhR is a ligand-activated transcription factor known to mediate the negative effects of environmental xenobiotic contaminants such as dioxin (i.e., TCDD; 2,3,7,8-tetrachlorodibenzo-[p]-dioxin). This receptor belongs to the basic-helix-loop-helix (bHLH)/PAS (Period [Per]-AhR nuclear translocator [Arnt]-Single minded [Sim]) family of heterodimeric transcriptional regulators [59]. There have been a number of studies in which physiological clues have been gathered as to the function of AhR. For example, AhR has shown tumour suppressor effects i.e., when receptor levels are down-regulated by siRNA [60] the G1/S transition of the cell cycle and cell proliferation is promoted, suggesting a growth inhibitory role of the receptor. There also exists links between AhR and angiogenesis. Vascular endothelial growth factor (VEGF), a major growth factor that regulates angiogenesis is transcriptionally regulated by hypoxia-inducible factor-1 alpha (HIF-1α) in response to tissue hypoxia and muscle contraction [61]. Hypoxia stabilizes HIF-1α, which forms heterodimers with HIF-1β, also known as the AhR nuclear translocator (ARNT). ARNT can heterodimerize with AhR, minimizing the ability of ARNT to interact with stabilized HIF-1 to induce VEGF production. Finally, matrix remodelling is negatively affected during AhR activation. For example, AhR activation blocks regenerative processes during zebra fish caudal fin regeneration while impairing expression of genes involved in the structure and remodelling of the ECM [62]. Thus, inactivation of AhR-signalling may facilitate ECM and vascular remodelling occurring in response to both RET and EET (Figure 1).

The notion that both ATRA and AhR pathways are potent regulators of cell growth, differentiation, ECM remodelling, vascularization, organogenesis and embryogenesis [63] underlines their key roles in tissue development and homeostasis. Moreover, since adaptation to RET and EET both involve aspects of cellular remodelling, muscle satellite cell activities [64], ECM and vascular remodelling, we suggest that ATRA and AhR molecular programs are playing a previously undefined but central role in regulating “generic” features of exercise adaptation [22]. Finally, since ATRA and AhR gene networks that were regulated during long-term exercise training (Figure 1), were not reflective of those modulated in the hours after a single bout of exercise [24], [33], this casts doubt over ascribing formative purpose of acute exercise gene networks, which more likely represent stress pathways instigated by unfamiliar activities or simply the acute energy crisis in exercised muscle (in agreement with the lack of a striking ontology profile). This may explain why acute mRNA changes do not overlap with the chronic exercise patterns or, in our hands, relate to the networks that associate with the degree of gain in lean mass (see below).

Networks activated in a manner directly related to muscle growth

Following the identification of what might be called generic ‘adaptability’ [22] molecular networks, we were interested to see if we could identify any molecular networks that were regulated in proportion to the degree of muscle hypertrophy in individual subjects. The justification for this approach was based on the marked heterogeneity in capacity for muscle growth in humans, with gains ranging from 0% to 22% [24], [28], [65]. In the DRET study, we found a similar range of changes in muscle size (−3% to +28%) and the analyses of the gain-related gene networks yielded striking results. We first discovered that there existed a correlation between capacity for human muscle lean mass gains and activity of both c-MYC and mTORc1 sensitive genes, including a large group of ribosomal RNA (rRNA) genes (∼70/560 total rRNAs). We further demonstrated that the nature of the association was not as one would be expecting from pre-clinical research, but rather there was a reduction in rRNA gene expression when greater muscle hypertrophy was observed (Figure 2B).

We speculate that human high-responders to hypertrophy potentially show superior efficiency in protein synthesis (i.e., protein yielded per RNA) and/or a reduction in proteolytic responses to RET. Nonetheless, it should also be noted that while rRNA expression negatively correlated to gain in lean mass, the highest responders tended to have higher levels of rRNA pre-training for individual rRNAs and hence the absolute abundance was similar post-training. This pattern of response clearly demonstrates that ‘more’ mTOR activation (e.g. more rRNA production) in humans is neither a hallmark nor a necessity for gains in lean mass in vivo. Our observations in humans also contrast with the molecular responses observed during acute synergist ablation induced hypertrophy [66], raising further doubts over the relevance of such pre-clinical models to inform about physiological muscle growth in humans. Regardless, our data demonstrate that high-responders for muscle hypertrophy evoke an “anti-growth” transcriptional response during a period of successful muscle growth with more studies being required in high or low lean-mass responders to investigate this phenomenon unambiguously.

While the unbiased transcriptional profiling yielded novel insights into human muscle growth responses, we also attempted to link training responsiveness to subject baseline physiology and the acute response of phospho-protein (AKT-mTORc1) signaling in response to an anabolic stimulus (RET combined with optimal nutrition [67]). We used correlation analysis to examine the relationship between key factors that have been speculated to be important for human muscle growth (e.g. body composition, fiber-type, metabolic and signaling molecule status). We found no significant shared variance between any of these variables and gains in lean mass, as presented in Figure 3 using PCA. We utilized PCA for visualizing the integration of physiological and molecular data because it enables an over-view in a single plot of how physiological and molecular aspects may vary within the major ‘units’ of variance of a particular dataset rather than plotting multiple individual scatter plots. By doing this, it becomes easier to visualize when paired relationships are behaving as expected, such as total protein and phospho-protein positions within each principal component. Unusual patterns could be used to identify problem areas such as when detection methods (e.g. antibodies) are poorly functioning.

As can be seen in Figure 3A, while principal component 1 was dominated by the variation in lean mass gains, basal lean mass, fat mass or fasting glucose did not vary. These data are in agreement with previous studies in which baseline physiological status was unable to select out high-responder status to hypertrophy [28]. We integrated the acute changes in phospho-protein signaling under the conditions of RET plus nutrition as an index of growth signaling potential in each subject. This seemed a valid approach as muscle hypertrophy is the product of nutrition and exercise-induced muscle protein synthesis [68]. While pre-training acute “anabolic” signals did not share variance with gains in muscle lean mass (Figure 3A), some weak relationships appeared between acute “anabolic” signaling elements following exposure to RET (Figure 3B) and gains in lean mass. Collectively, this underlines that while acute changes in phosphorylation may associate with those of acute remodeling processes (i.e., muscle protein synthesis responses) they are a poor indicator of future leans mass gains. Previously we have found a degree of dissociation between AKT-mTOR signals [69], [70] and muscle protein synthesis, while others have shown that acute synthesis does not, but some signaling molecules do, relate to future gains in lean mass [71]. We contend that using these individual signals as acute proxies for RET muscle growth is not going to be the most sensitive strategy.

Age-associated and exercise-associated molecular networks

There is a long established relationship between canonical pathways related to muscle growth and mechanisms that are associated with extended life-span in model organisms [72]. For example, inhibitors of mTORc1 and PI3K (a gene also known as “age-1”) activity can extend life-span in Caenorhabditis elegans, drosophila and mice [72]. As the number of people living beyond their eighth decade rises, it is expected that skeletal muscle atrophy and dysfunction (sarcopenia and dynapenia respectively) will become an increasing public health challenge [73]. Activation of mTORc1 and muscle IGF-1 signalling is associated with muscle cell growth experimentally while chronic inhibition of mTOR has been predicted to induce muscle frailty in humans [14], [74]. Likewise, reduced S6K1 (ribosomal protein S6 kinase) activity of the multifaceted regulator of cell growth would be thought to impair the retention of muscle mass in humans [75]. Thus there is a clear molecular basis to believe that processes under-pinning ageing, longevity, sarcopenia and muscle growth will be strongly inter-connected.

We found a list of ∼500 genes which track with age in human muscle across two independent cohorts, when our analyses utilised a continuum of ages. This contrasts with the irreproducibility of previous RNA versus muscle age datasets. Using this data we then evaluated which age-related gene networks link to muscle growth signalling or a variety of exercise scenarios [24], [29], [56], [76]. The motivation for this is that these factors cannot be independently controlled for in human studies and thus post hoc considerations are essential. The first clear observation was that the reproducible age-related gene-list and the lean-mass gain related gene-list both had inhibition of MYC as a key transcriptional feature immediately indicating that age and muscle growth are not exact opposites. Further, the key upstream regulators of the age-related gene list (e.g. PGR, RXR, Claudin7) contained a set of genes which were largely unrelated to the RET and EET transcriptomes, or the acute exercise responses [24], [76] for that matter. Claudin-7 activity was inhibited with age (Figure 6) and appears to relate to developmental differentiation and is strongly regulated by androgen signaling and HGF in vitro [77]. According to the IPA publication database the RXR ligand, CD437, induces a transcriptional signature that is consistent with age-related gene correlations, involving a network that is for 2/3^rd negatively correlated with increasing age and for 1/3^rd positively associated with age (i.e. with age KLF5 and IRS2 expression increases). Note that the ATRA signature, characteristic of increased physical activity in our hands, should operate through the RAR pathway. However, there is some overlap in the down-stream gene activation and thus both age-related changes and exercise activate some common features related to Vitamin A biology. In short, it is abundantly clear that the age-related changes in gene expression are not simply the ‘opposite’ [23] of the profiles seen with exercise or exercise training.

Positional enrichment analysis [43] was used to map the ∼500 age-related genes to chromosomal localisation, as attempts to link DNA variants with ageing have so far been only partly successful suggesting that an alternative approach may provide useful insights. We found several loci that yielded a significant enrichment score with ∼1q12 (e.g. PDE4DIP) and ∼13q21 (e.g. LMO7 (FBXO20)) yielding the most significant scores (Figure 7). PDE4DIP is a binding partner of phosphodiesterase 4D (PDE4D) which partners with Rheb to be a cAMP-specific negative regulator of mTORc1 [78]. When PDE4D binds Rheb it inhibits the ability of Rheb to activate mTORC1 and hence it is plausible that PDE4DIP (also called myomegalin) impacts on this relationship, as an increased interaction between Rheb and mTOR should promote growth. On chromosome 13, another growth and differentiation gene, Lmo7, was identified. Lmo7 impacts on myoblast differentiation, being required to induce Pax3 and MyoD expression [79]. Lmo7 is positively associated with age in muscle, while we have previously identified PAX3 as an up-stream transcriptional regulator of the EET induced transcriptome [29], highlighting further why the age-related transcriptome is not simply the opposite of exercise training. i.e., muscle ageing is not simply inactivity and thus is unlikely to be reversed by only activity (this does not at all contest that muscle function can be substantially retained by physical training in many but not all subjects).

Perspectives and limitations

Over the past 2 decades, therapeutic advances for complex chronic diseases have failed to generate all the progress predicted by the emergence of genetic technologies [80]. Part of the reason for this is that the investment in forward genetic pre-clinical models [81] has not yielded the expected validation of drug targets and it is now widely accepted that the clinical approach to chronic-disease management will have to reflect on numerous interactions between environmental and molecular factors [82]. Thus, an alternative approach, whereby identification of disease networks directly observed in clinical populations may have merit and lead to a more rapid translation of basic science [3], [4], [34], [40], [83]. Validation of such observations is dependent on access to sets of independent clinical data sufficiently large to have robust statistical power and diverse enough to be able to generalize the conclusions. It is safe to say that our approach is not at all universally favoured, as there is still a great reliance on so-called ‘validation’ studies involving forward and reverse genetic strategies in mice and cells. We suspect that such validation studies will, in the end prove to be context dependent and no easier to interpret than our in vivo molecular studies.

The singular advantage of our approach is that our data is generated under the precise conditions that would ultimately require therapeutic intervention. However, in the present analysis we focus on components of dynamic function, namely muscle mass, as the study of muscle performance will require determination of a wider range of muscle functional parameters (power, torque, velocity) and larger clinical studies, studies which do not exist at present. We also appreciate that our current observations would benefit greatly from follow-on genetic association studies in humans and by pharmacological or nutritional intervention studies. For example we have shown that a simple relationship between more mTOR activation and muscle growth does not exist in vivo in humans, albeit we are relying in part on the pharmacological activity of rapamycin to support this observation (off target effects may be present) and do not yet have kinetic data to understand the dynamic nature of this new relationship. We also failed to find a link between inter-subject variation in acute phospho-protein anabolic signalling and gains in lean mass. This may reflect the choice of time-point that we profiled the protein responses under. However, it must also be recognised that quantification of protein abundance changes suffers from numerous complications, including compression or exaggeration of dynamic range and challenges with specificity. Time will tell if our systems-biology translational medicine approach [84] exceeds the performance of traditional approaches taken to yield new therapeutic advances for human health.

Materials and Methods

This manuscript relies on several microarray studies to produce and independently validate molecular associations. We present a new 89-chip analysis (44 paired samples and 1 baseline-only sample) combined with direct comparisons to literature-based microarray data [23], [24]. The Trappe Lab chip data (GSE28422) was processed using identical methods to our new clinical data. Our QC analysis identified a number of chips that should be removed prior to statistical analysis. Inclusion of such data can both increase or decrease detection of differential gene expression (derived from technical and not biological variation). The analysis as a whole represents the most comprehensive examination of human molecular exercise-ageing pathways, and the application of independent chip resources provides a gold-standard level of validation.

Derby RET (DRET) subject characteristics

Subjects were recruited from an age range of 18 to 75 y. Before beginning the study all subjects were screened using a medical questionnaire, physical examination and resting ECG with exclusions for overt muscle wasting (>2 SD below age norms) [85], metabolic, respiratory/cardiovascular disorders or other major contraindications to a healthy status. All subjects had normal blood chemistry and were normotensive (BP<140/90). All subjects performed routine activities of daily living and recreation but did not participate in moderate to high intensity aerobic exercise and none had participated in RET in the last 24 months. Body composition was measured at screening and following RET by dual energy X-ray absorptiometry (DEXA) (Lunar Prodigy II, GE Medical Systems). Subject positioning on the DEXA bed was optimized to allow the region of interest (ROI) body compartments to be analyzed separately. The upper leg ROI was selected as the area inferior to the lowest visible point of the coccyx to the mid-point of the patella. All subjects gave their written, informed consent to participate after all procedures and risks were explained. This study was approved by the University of Nottingham Medical School Ethics Committee and complied with the Declaration of Helsinki. The clinical data from these subjects were first reported in 2012 [56]. For the purposes of this article a total of 45 from the original 51 subjects were utilised as this represents the total number of gene-chip profiles that passed the appropriate quality control processes (N = 89 U133+2 Affymetrix chips).

DRET protocol

The 20-wk fully supervised RET programme was designed to achieve skeletal muscle hypertrophy. Subjects trained three times per week, with each session lasting approximately 60 min. During four weeks of induction training (to ensure adoption and adherence to correct technique) intensity was increased from 40% to 60% 1-RM. For the remaining 16-wk of training intensity was set at 70% 1-RM with multiple sets of 12 repetitions, with two min rest between sets. 1-repetition maximum (1-RM) assessments were made every four weeks to ensure that the intensity of loading was constant. Subjects were excluded from the study for non-compliance, defined as: non-attendance for >6 consecutive sessions, less than 75% attendance, or failure to complete the set exercise regime on >15% attendance. Muscle biopsies (∼150 mg) were taken under fasted-non-exercised (“basal”) and optimal growth conditions (acute exercise-fed conditions, 2.5 h after a single bout of exercise) both before and after chronic-RET from the vastus lateralis muscle under local anaesthesia (2% lidocaine, with the use of a conchotome biopsy forceps, as previously described [56]). Blood was collected in pre-chilled tubes containing Lithium Heparin, plasma was separated by centrifugation and was then stored at −80°C until analyses. Plasma glucose concentration was measured on a clinical chemistry glucose analyser (ILAB 300 Plus). For insulin, blood was collected in pre-chilled tubes containing EDTA, plasma was separated by centrifugation within 30 min of collection and was then stored at −80°C until final analyses. Plasma insulin concentration was determined using high sensitivity insulin ELISA systems (Immunodiagnostic systems limited).

Immuno-blotting

Muscle biopsies (∼20 mg) were homogenized in ice-cold extraction buffer (10 µL.mg⁻¹) containing 50 mM Tris-HCl (pH 7.4), 0.1% Triton X-100, 1 mM EDTA, 1 mM EGTA, 50 mM NaF, 0.5 mM activated sodium orthovanadate (Sigma Aldrich, Poole, UK) and a complete protease inhibitor cocktail tablet (Roche, West Sussex, UK). Homogenates were centrifuged at 10,000 g for 10 min at 4°C. Bradford assays were used to determine supernatant protein concentrations after which samples were standardized to 1 mg.mL⁻¹ in Laemmli loading buffer. Samples were heated at 95°C for 5 min before 15 µg of protein/lane was loaded on to Criterion XT Bis-Tris 12% SDS-PAGE gels (Bio-Rad, Hemel Hempstead, UK) for electrophoresis at 200 V for ∼60 min. Gels were equilibrated in transfer buffer (25 mM Tris, 192 mM glycine, 10% methanol) for 30 min before proteins were electro-blotted on to 0.2 µm PVDF membranes (Bio-Rad) at 100 V for 30 min. After blocking with 5% low-fat milk in TBS-T (Tris-Buffered Saline and 0.1% Tween-20; both Sigma-Aldrich, Poole, UK) for 1 h, membranes were rotated overnight with primary antibody (all AbCam, Cambridge, UK) against our target proteins (AKT, mTOR, p70S6K1, 4EBP1, eEF2) at a concentration of 1∶2000 at 4°C. Membranes were washed (3×5 min) with TBS-T and incubated for 1 h at room temperature with HRP-conjugated anti-rabbit secondary antibody (New England Biolabs, UK), before further washing (3×5 min) with TBS-T and incubation for 5 min with ECL (Immunstar; Bio-Rad). Blots were imaged and quantified by assessing peak density, after ensuring bands were within the linear range of detection using the Chemidoc XRS system (Bio-Rad, Hemel Hempstead, UK). Protein phosphorylation was corrected for loading to actin loading control before the protein signals were subject to PCA to explore relationships between ‘anabolic signals’, specifically in terms of muscle hypertrophy responsiveness.

RNA extraction and chip processing

Total RNA was isolated from muscle biopsies taken before and after (72 h after the final training session) RET by chloroform-phenol based extraction. In brief, paired tissue samples (obtained before and after RET) of ∼20 mg each were processed simultaneously in 1 mL TRIzol (Invitrogen) using a Mini-Beadbeater-8 (Biospec Inc.) for 15 sec on the “homogenize” setting. After 5 min of incubation at room temperature, 200 µL of chloroform (Sigma-Aldrich) was added and samples shaken vigorously by hand. Samples were briefly incubated on ice prior to centrifugation at 12,000 g for 15 min. The supernatant was removed and mixed with isopropanol (Sigma-Aldrich) and spun once more at 12,000 g for 10 min, after 10 min of incubation. After a single washing step with 75% EtOH RNA pellets were re-suspended in 40 µL DEPC-treated water (Ambion) and quantified using a NanoDrop Spectrophotometer (NanoDrop Technologies). RNA purity was assessed using the A₂₆₀/A₂₈₀, A₂₆₀/A₂₃₀ ratios and stored at −80°C. Samples were put through this process in pairs (pre-post samples) while the order of subject processing was carried out to distribute ‘non-responders’ equally.

Reverse transcription of RNA was carried out using the Affymetrix 3′ IVT express kit. 100 ng of total RNA was reverse transcribed as per manufacturer's protocol, and quantified using a Nanodrop ND-1000 instrument. aRNA was fragmented and labeled as per manufacturers protocol and hybridized to Affymetrix U133+2 arrays (Affymetrix, USA). Arrays were washed, stained and scanned following Affymetrix standard procedures, using an Affymetrix 3000 7G scanner and Affymetrix 450 wash station. A visual inspection of each array was carried out.

Microarray analysis and bioinformatic procedures

Low-level processing of all arrays was undertaken using Bioconductor in R. The Affy package was used to carry out MAS5 based normalization and generate present, marginal and absent (PA) scores. NUSE plots were generated and combined with PCA, outlier samples were identified where both the NUSE plot and PCA was supportive of its exclusion (∼2% of arrays). For baseline correlation analysis, all samples that passed QC were utilized (N = 45). This procedure was applied to the data-set originating from the Trappe laboratory (GSE28422, [24]) and outliers removed from the dataset that failed the QC process, leaving n = 96 for analysis. Pre-exercise training muscle biopsy samples from the HERITAGE family study (N = 50) were also analyzed to yield a second independent data set with a continuous span of age-ranges (see below). The Trappe and HERITAGE datasets therefore represent independent datasets which we utilized, where possible, to validate the pathway analysis of our study. Such confirmation benchmarks results using thousands of data-points and is more desirable that targeted real-time qPCR confirmation (where the gene selection is biased and the sample size inappropriate to make statistical conclusions). Annotation of all CEL files used ‘hgu133plus2cdf_2.9.1.tgz’ while annotation of probe-set lists was then updated using the Ingenuity Pathway Analysis database, as of August 2012.

Our first objective was to identify the gene-networks regulated in proportion to gains in upper leg muscle mass (hypertrophy), the same location as our biopsy sample. Such analysis relies on the established principal that adaptation responses (for the majority of phenotypes) to exercise training in outbred populations is highly variable, typically reflecting genetic and epigenetic variation and in genomic variation. We utilized quantitative SAM analysis [3], [36] to generate a list of genes which vary in a positive and negative manner with changes in DEXA assessed upper leg lean mass. This was applied to PA filtered data and the statistical parameter generated is a q-value (false discovery rate). This provided for the first time a candidate list of gene-changes that may exhibit primary or secondary influence over muscle growth in humans. The gene-list was then subject to IPA based pathway analysis and in particular the Upstream Analysis tool in IPA was utilized. This analysis has similarities to the Molecular Connectivity Database [86] where pre-existing collections of RNA signatures are compared with our lean-mass related gene list, and significant overlaps identified. An overlap P-value is generated based on the degree of overlap between the gene-set within the IPA database (which reflects the RNA molecules changed in response to a ‘mediator’ such as a transcription factor or a drug) and our data set, adjusting for data set sizes using the Fischer's Exact Test. We accepted a stringent P-value of p<0.001 as being significant. A second parameter is the activation “z-score” where the directional change in RNA is compared between the IPA mediator data-set and our lean-mass gain gene list. The z-score informs on whether the drug/protein mediator is likely to be ‘active’ or ‘inhibited’ during gains in lean mass. Thus, if we discovered that an antagonist is ‘inhibited’ in our analysis, this indicates that the drug target is activated. However, in the present study the data-input refers to genes, which positively or negatively correlate with lean mass gains e.g. if we find a “Statin” signature was inhibited, it is interpreted that HMG-CoA reductase regulated genes are negatively correlated with lean mass gain.

The two-step process presented above generates a focused gene-list with a high statistical rigor for true positive associations. This type of analysis also utilizes the full range of physiological response observed, however it assumes that expression of important genes will relate in a linear manner to lean mass gain and thus can not discover all appropriate associations. We then contextualize the statistical findings both in terms of subject characteristics and through comparison of the response of these significant networks with independent gene-array data (e.g. [24]). At this stage we utilized descriptive statistics, plotting the significant network genes as simple expression values relative to the quartile distribution of lean-mass gains to allow for clear discussion of the results. As these plots are based on the z-scores and P-values as above, no further statistical analysis is presented.

Following identification of our primary objectives we then carried out a classic differential expression analysis using SAM. Given that we have established that chronic differential expression patterns, following exercise training, are dependent on the presence of physiological adaptation we removed 6 subjects that demonstrated no gain in lean mass. This yielded a list of differentially expressed genes that could then be compared with the RET gene-list generated from the Trappe laboratory data and our published exercise studies [24]. Secondary analysis, where subject age or baseline lean-mass was related to baseline gene-expression was carried out using quantitative SAM analysis as described above [3], [36]. This allowed us to present comparisons of the RET gene-list with other modes of exercise, such as endurance exercise training or disease [3], [29] and age-related analysis [23], [24]. PCA was utilized to visualize the association between selected physiological and protein expression parameters and training induced changes in muscle lean mass. PCA was implemented in R, using prcomp() command, which calculated a singular value decomposition and plots the selected principal components using the plot command in R. All data was individually transformed to a median value within that data set so that all variables were within a consistent data range. In each case the majority (∼65%) of the total variance was captured by the first two principal components.

Finally, positional gene enrichment analysis (PGE) was used to identify whether the classification genes (or the classifier network genes) were significantly enriched within given chromosomal regions [43]. This analysis is based on the following rules: Rule 1: it contains at least two genes of interest, Rule 2: there is no smaller region containing the same genes of interest, Rule 3: there is no bigger region with more genes of interest and the same genes not of interest, Rule 4: there is no larger encompassing region with a higher percentage of genes of interest, Rule 5: there is no smaller encompassed region with a better P-value, Rule 6: it does not contain any region having less than expected genes of interest. The approach of PGE exhaustively evaluates the over-representation at all chromosomal resolution levels simultaneously.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10