Power training and postmenopausal hormone therapy affect transcriptional control of specific co-regulated gene clusters in skeletal muscle.

At the moment, there is no clear molecular explanation for the steeper decline in muscle performance after menopause or the mechanisms of counteractive treatments. The goal of this genome-wide study was to identify the genes and gene clusters through which power training (PT) comprising jumping activities or estrogen containing hormone replacement therapy (HRT) may affect skeletal muscle properties after menopause. We used musculus vastus lateralis samples from early stage postmenopausal (50-57 years old) women participating in a yearlong randomized double-blind placebo-controlled trial with PT and HRT interventions. Using microarray platform with over 24,000 probes, we identified 665 differentially expressed genes. The hierarchical clustering method was used to assort the genes. Additionally, enrichment analysis of gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was carried out to clarify whether assorted gene clusters are enriched with particular functional categories. The analysis revealed transcriptional regulation of 49 GO/KEGG categories. PT upregulated transcription in "response to contraction"-category revealing novel candidate genes for contraction-related regulation of muscle function while HRT upregulated gene expression related to functionality of mitochondria. Moreover, several functional categories tightly related to muscle energy metabolism, development, and function were affected regardless of the treatment. Our results emphasize that during the early stages of the postmenopause, muscle properties are under transcriptional modulation, which both PT and HRT partially counteract leading to preservation of muscle power and potentially reducing the risk for aging-related muscle weakness. More specifically, PT and HRT may function through improving energy metabolism, response to contraction as well as by preserving functionality of the mitochondria.

RCT Entities:   Population Interventions Outcomes

Introduction

Menopause, the transition from the reproductive to non-reproductive stage, is a complex physiological process, often accompanied by additional effects of aging (Nelson 2008). In addition to reproductive viability, deterioration in musculoskeletal properties begins to accumulate (Kallman et al. 1990; Phillips et al. 1993; Samson et al. 2000), which may lead to severe consequences for the quality of life and the ability to recover from physical traumas ultimately reducing the healthspan of women. Therefore, early preventative actions to maintain adequate muscle performance in later life are needed.

So far the best prevention strategy against aging-related weakness is physical exercise (Greenlund and Nair 2003; Taaffe 2006). Strength, endurance, and power training have been found to exhibit beneficial albeit differential effects on aging musculature (e.g., Hoppeler et al. 1985; Häkkinen et al. 2000; Roth et al. 2001; Timmons et al. 2005; Porter 2006; Hazell et al. 2007; Henwood et al. 2008; Orsatti et al. 2008). Progressive power training such as undertaken in the current study, is associated with enhanced quality of life, i.e., physical functioning and well-being (Katula et al. 2008), neuromuscular functions involving rapid force production (Sipilä et al. 2001; Taaffe et al. 2005), and bone formation (Cheng et al. 2002). Moreover, muscle power is essential for many daily tasks such as climbing stairs, rising from chair, and preventing a fall after a slip leading a number of recent studies (Hazell et al. 2007; Paterson et al. 2007; Katula et al. 2008; Orr et al. 2008) to conclude that strategies improving muscle power and composition are essential in order to prevent old age disability. Therefore, power training is promising training mode for postmenopausal women whose risk for osteoporosis and neuromuscular dysfunction might otherwise be increased. However, the genome-wide effects of this specific type of training on skeletal muscle have not, until now, been investigated.

In addition to physical exercise, the restoration of female sex hormones to near pre-menopausal levels may represent one way to attenuate functional decline and muscle loss after menopause. Indeed, several randomized controlled trials (RCT) have shown that hormone replacement therapy (HRT) increases muscle mass and improves neuromuscular performance when used in the beginning of the postmenopausal period (Skelton et al. 1999; Sipilä et al. 2001; Dobs et al. 2002; Taaffe et al. 2005). In addition, a co-twin study, in which monozygotic HRT users were compared to their non-using sisters, showed that long-term HRT treatment is associated with greater muscle power, higher walking speed, and decreased total body and thigh fat content (Ronkainen et al. 2009). However, not all trials with HRT have reported beneficial effects on muscle and physical function (Ribom et al. 2002; Tanko et al. 2002; Kenny et al. 2005) indicating the need for additional research to clarify this discrepancy.

In our yearlong RCT, the exercise and hormone replacement therapy-study (Ex/HRT), the effects of plyometric power training and HRT on bone as well as muscle structure and function in postmenopausal women were studied (Sipilä et al. 2001; Cheng et al. 2002; Taaffe et al. 2005). We found a significant 6% increase in muscle power with power training and 7% with HRT, while a 5% reduction was observed in control participants without any treatments (Sipilä et al. 2001). Running speed increased by 4% both after HRT and power training and reduced by 2% in the controls (Taaffe et al. 2005). Moreover, knee extensor muscle density increased significantly by both power training and HRT (Taaffe et al. 2005) reflecting beneficial changes in muscle composition in favor of muscle tissue itself on the expense of adipose tissue (Goodpaster et al. 2000). Knee extensor muscle cross-sectional area increased by 6% with HRT (Sipilä et al. 2001). Our earlier microarray study investigating the effects of 1-year HRT use showed that the postmenopausal women without treatment have notable changes in their muscle transcriptome and that these changes were largely counteracted by the use of HRT (Pöllänen et al. 2007). To the best of our knowledge, other microarray studies on the effects of HRT or power training on skeletal muscle of postmenopausal women have not previously been undertaken.

The purpose of this ancillary genome-wide analysis of the Ex/HRT study was to identify gene clusters that are affected by progressive power training, use of HRT or recent menopause, and thus may induce the adaptive changes in skeletal muscle characteristics. To accomplish our goal, we combined previously obtained microarray data with completely new data from eight power trainers and four control women (totally ten HRT users, eight trainers, and nine non-treated women). This allowed us to identify both unique and overlapping effects of the two treatments on muscle transcriptome.

Materials and methods

Study design and interventions

The detailed description of the study design has been reported by Sipilä et al. 2001. The trial is registered in the Current Controlled Trials with ISCTN number ISRCTN49902272 (http://www.controlled-trials.com/ISRCTN49902272). Briefly, we enrolled 1,333 women at age 50–57 years living in the area of the municipality of Jyväskylä in Finland. From these women, 118 went through an extensive medical and physical examination including determination of their menopausal status as well as eligibility to the study. To be eligible for the trial, participants had to have no serious medical conditions; no current or previous use of medication including estrogen, fluoride, calcitonin, biophosphonates, or steroids; last menstruation at least 0.5 years but not more than 5 years ago; FSH level above 30 IU/L; and no contraindications for exercise and HRT. Finally, we randomly assigned 80 women fulfilling the inclusion criteria, i.e., being at the very early stage following menopause, to one of four study groups: power training (PT, n = 20), HRT (n = 20), PT+HRT (n = 20), and control (CO, n = 20).

For the current microarray analysis, baseline and 12-month muscle samples were available from eight PT, ten HRT, and nine CO women. Seventeen participants were lost during the intervention, and muscle samples from both baseline and post-intervention (12-months) were not available from 14 subjects. Moreover, one participant did not consent for the microarray study. The analysis does not include the PT+HRT group due to the small number of tissue samples available (n = 6). The participants available for the microarray study did not differ significantly from the whole study sample in any of the tested muscle phenotype variables at baseline or in changes during the intervention.

Muscle biopsies obtained from the mid-part of musculus vastus lateralis (midpoint between the greater trochanter and the lateral joint line of the knee) were taken at baseline and after completion of the study as previously described (Pöllänen et al. 2007). The biopsy protocol was standardized within the study in order to avoid sampling-induced variation. All muscle biopsies were taken from the side of the dominant hand approximately at the same time of day by the same experienced physician. The 12-month biopsy was taken 1 cm above the previous biopsy cicatrix within 3 to 5 days after completion of the intervention. The histological evaluation of the samples did not reveal any signs of damage, such as central nuclei, in any of the samples.

The detailed description of the interventions is reported in Sipilä et al. 2001. In brief, PT participants underwent a progressive plyometric training program for the lower limbs, which comprised two supervised sessions per week and four home-based unsupervised sessions per week. The training was performed in a circuit format including bounding, drop jumping, hopping, and skipping performed at high velocity in order to improve muscle power production and to produce high-impact loading for bones. The training program progressed in the number of rotations performed, volume of work undertaken as well as height of obstacles for bounding and height for drop jumping. Each supervised session included three to four resistance training exercises for the upper body and commenced with a warm-up period and concluded with a cool-down period of stretching activities. HRT and CO subjects were advised to maintain their daily routines without altering their physical activity patterns. The HRT treatment was conducted double-blind. All study participants used either continuous, combined HRT preparation containing estradiol (2 mg) and or noretisterone acetate (1 mg) preparation (Kliogest, Novo Nordisk, Copenhagen, Denmark) or placebo (composed of lactose monohydrate, cornstarch, gelatin, talc, and magnesium stearate, which were auxiliary substances in the Kliogest tablet) one tablet every day.

The study was carried out in accordance with the Declaration of Helsinki (48th WMA General assembly; 1996 1989) of the World Medical Association (www.wma.net) and was approved by the ethics committee of the Central Finland Health Care District. Informed consent was given by all subjects.

RNA preparation, cRNA generation, and array hybridization

The RNA preparation, cRNA generation, and array hybridization procedures have been previously described (Pöllänen et al. 2007). Briefly, Trizol-reagent (Invitrogen, Carlsbad, CA, USA) was used to isolate RNA from biopsy samples homogenized on FastPrep FP120 apparatus (MP Biomedicals, Illkrich, France). The Experion (Bio-Rad Laboratories, Hercules, CA, USA) was used to inspect RNA concentration and quality. Only pure, good-quality RNA was used in further analysis (260/280 ratio >1.8). An Illumina RNA amplification kit (Ambion, Austin, TX, USA) was used according to the manufacturer's instructions to obtain biotin-labeled cRNA from 500 ng of total RNA. Experion was used to perform quality control after amplification. Hybridization to the HumanRef-8 v1.0 or HumanWG-6 v1.0 BeadChips (Illumina Inc., San Diego, CA) as well as washing and scanning was performed according to the Illumina BeadStation 500x manual (revision C). Both samples (baseline and follow-up) from each study subject were always hybridized onto the same chip. The slides were scanned by confocal laser scanning system (Illumina BeadReader Rev. C, Illumina Inc., San Diego, CA, USA). The data were acquired by the BeadStudio Direct Hybridization V.1.5.0.34. The Turku Centre for Biotechnology (Turku, Finland) carried out the cRNA generation, array hybridizations, and quality control of the raw data. The data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus (GEO; Edgar et al. 2002) and are accessible through GEO Series accession number GSE16907 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE16907). The MIAME guidelines were followed during array data generation, preprocessing, and analysis.

Array data preprocessing and validation

Raw gene expression data from eight PT and four CO participants were obtained by using Human WG-6 BeadChips (Illumina), while HumanRef-8 BeadChips (Illumina) were utilized for ten HRT and five CO participants. Raw data from both BeadChips were combined, and only the probes, which are identical in both platforms, were included in further analysis. Consequently, over 24,000 probes for approximately 21,000 NCBI RefSeq-transcripts were included. Signal for each transcript comes from approximately 30 independent technical replicates (beads) as described previously (Pöllänen et al. 2007). Quality control and visualization of the raw, combined expression data were performed with aid of box plots, hierarchical clustering, correlation matrix, and principal component analysis. Microarray data from the three study groups, i.e., PT, HRT, and CO, were normalized separately using the quantile normalization method implemented in the Affy-package (Gautier et al. 2004) of the R/Bioconductor analysis software (www.r-project.org, www.bioconductor.org). Samples (n = 6) from three CO subjects were hybridized onto both Illumina platforms in order to compare the performance of previous and current arrays with the Pearson correlation coefficient test. Even though the correlation of the data produced by the same samples on different Illumina platforms was high (r = 0.88–0.94), the clear batch effect caused by two different platforms was visible. This was corrected in statistical analyses by using array type as a covariate.

Data analysis

After preprocessing the microarray data, the Limma package (Smyth 2004) of R was used to identify differentially expressed genes within each condition. The fold changes (FC) and p values calculated by Limma methods were used as filtering criteria in order to detect genes up- or downregulated within each study group. The repeated measurements design and the batch effect caused by the usage of two different array types were taken into account in statistical testing. At this point, we chose to use quite low stringency thresholds (p value <0.05 and |FC| >1.2) in order to list a reasonable number of the genes most likely to be differentially expressed. Decreasing or increasing the thresholds resulted in poorer cluster structures in the clustering of the filtered genes. Genes determined to be differentially expressed within at least one of the study groups were used in further analyses, which allow identifying co-regulated gene clusters, i.e., genes, which show consistent behavior being either up- or downregulated within each study group. We applied hierarchical clustering with Pearson's metrics to form clusters of most similarly behaving genes. The Pearson's correlation thresholds were set to 0.95. Additionally, extensive analysis of functional categories, i.e., gene ontology (GO) terms (Ashburner et al. 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al. 2006), was carried out to clarify whether certain gene clusters are enriched with particular functional categories. Hypergeometric test and false positive discovery rate calculation were applied to identify the most significant categories. GenoSyst Ltd carried out the hierarchical clustering with enrichment analysis. The results obtained were further analyzed and processed in our lab. PatternViewer 2.39 (GenoSyst Ltd, Turku, Finland) was used to visualize the results.

Results

Co-regulated gene clusters

Differentially expressed genes were defined as genes that differed between follow-up and baseline samples within a study group with statistical threshold p value <0.05 and |FC| >1.2. These thresholds led to the discovery of 665 genes. In the PT group, 328 genes were upregulated and 182 downregulated, whereas in the HRT group, 34 genes were upregulated and 20 downregulated. The corresponding numbers for the CO group were 91 and 70, respectively. A non-redundant list of differentially expressed genes was used in the analysis of co-regulated gene clusters. Genes were sorted into eight significant and consistent clusters (Fig. 1, Table 1), which showed either study group-specific regulation or regulation into the same direction within all study groups. As seen in Fig. 1, genes within each cluster had a highly consistent expression pattern forming congruent cluster structure. The complete list of genes belonging to each cluster is presented in Table S1 as supporting information. For clarity, the official or commonly used acronym for gene names are used in the text and complete gene names with proper annotation are provided as supporting information in Table S2.

Fig. 1

Graphical view of the cluster structure within each cluster (CL #1–CL #8). Each line represents the mean fold change (FC) of a single gene

Table 1

Differentially expressed genes formed eight distinctive clusters

Clusters	Classification of the clusters	Number of genes in cluster	Number of genes with GO annotation	Number of genes with KEGG annotation	Number of unclassified genes
CL #1	Uniform upregulation	224	199	61	23
CL #2	PT-specific upregulation	82	69	20	13
CL #3	Uniform downregulation	139	119	46	20
CL #4	Upregulation in the PT and HRT groups	81	71	22	9
CL #5	HRT-specific upregulation	45	41	21	4
CL #6	Downregulation in the PT and HRT groups	43	19	3	24
CL #7	HRT-specific downregulation	31	24	12	7
CL #8	PT-specific downregulation	20	12	5	8

Enrichment analysis

The identified gene clusters (CL #1–8) were diligently investigated in order to clarify whether they were enriched with particular GO terms or KEGG pathways. The algorithms used in the analysis allow determination of the percentage of false positive discovery (FP-%), which enables monitoring the expected proportion of false positives among the discovered patterns. The FP-% is obtained by comparing the permutated (at least 1,000 simulations) cluster-annotation pairs to the actually observed cluster-annotation pairs, and it allows efficient and visual selection of the significant cut-offs with an acceptable level of false discovery (Fig. 2). Due to the hypothesis-generating rather than hypothesis-proofing nature of the microarray studies, we selected FP-% <15%. For the GO terms, a p value <0.0005 was used as the cut-off value. With this cut-off, there are 26 significantly enriched GO terms from which theoretically three might be false positives (FP-% ∼12%). For the KEGG pathways, a p value <0.005 is sufficient to tolerate similar proportion of false positives (FP-% ∼10%). This cut-off value identifies 23 significant KEGG pathways. The results of the enrichment analyses are presented in Tables 2 and 3 and thoroughly examined below with a summary of the main results presented in Fig. 3.

Fig. 2

Distribution of false positive rate (FP-%) at different p values. The horizontal black line represents the selected threshold FP-%. The p values <0.0005 and <0.005 are sufficient to provide FP-% <15% confidence level for enrichment analysis of GO and KEGG, respectively, functional categories

Table 2

Enrichment of GO terms into the co-regulated gene clusters. Note, numbers of genes are presented as they occur in the GO terms, meaning that they may be included in several different terms. Corrected p value is obtained by multiplying by the number of hypothesis tested

GO ID	Term description	Number of genes in cluster/array	p value	Corrected p value
CL #1: uniform upregulation, 199 GO-annotated genes
GO:0006936	Muscle contraction	7/84	3.29E−5	0.019
GO:0006941	Striated muscle contraction	5/36	3.86E−5	0.023
GO:0005977	Glycogen metabolic process	4/21	6.63E−5	0.039
GO:0007517	Muscle development	6/85	3.04E−4	0.179
GO:0007519	Skeletal muscle development	4/31	3.21E−4	0.189
GO:0004368	Glycerol-3-phosphate dehydrogenase activity	2/3	3.40E−4	0.201
GO:0030731	Guanidinoacetate N-methyltransferase activity	2/3	3.40E−4	0.201
GO:0016600	Flotillin complex	2/3	3.40E−4	0.201
CL #2: PT specific upregulation, 69 GO-annotated genes
GO:0003779	Actin binding	9/329	3.42E−6	0.001
GO:0005737	Cytoplasm	29/4142	1.91E−4	0.050
CL #3: uniform downregulation, 119 GO-annotated genes
GO:0047115	Trans-1,2-dihydrobenzene-1,2-diol dehydrogenase activity	3/4	1.02E−6	0.001
GO:0004033	Aldo-keto reductase activity	3/10	2.98E−5	0.013
GO:0047026	3-Alpha-hydroxysteroid dehydrogenase A-specific activity	2/2	4.07E−5	0.018
GO:0047042	3-Alpha-hydroxysteroid dehydrogenase B-specific activity	2/2	4.07E−5	0.018
GO:0006937	Regulation of muscle contraction	3/12	5.41E−5	0.024
GO:0005739	Mitochondrion	17/944	1.07E−4	0.048
GO:0045737	Positive regulation of cyclin-dependent protein kinase activity	2/3	1.22E−4	0.054
GO:0048306	Calcium-dependent protein binding	3/17	1.63E−4	0.073
O:0032052	Bile acid binding	2/4	2.42E−4	0.108
GO:0015721	Bile acid and bile salt transport	2/4	2.42E−4	0.108
GO:0001527	Microfibril	2/4	2.42E−4	0.108
GO:0016491	Oxidoreductase activity	11/481	2.67E−4	0.119
GO:0005578	Proteinaceous extracellular matrix	8/265	3.09E−4	0.138
GO:0031072	Heat shock protein binding	4/57	4.87E−4	0.217
CL #5: HRT-specific upregulation, 41 GO-annotated genes
GO:0005739	Mitochondrion	12/944	5.56E−7	0.001
CL #7: HRT-specific downregulation, 24 GO-annotated genes
GO:0001558	Regulation of cell growth	3/114	4.15E−4	0.052

Table 3

Enrichment of KEGG pathways into the co-regulated gene clusters. Note, the number of genes is presented as they occur in the KEGG pathways, meaning that they may be included in several different pathways. Corrected p value is obtained by multiplying by the number of hypothesis tested

KEGG ID	Term name	Number of genes in cluster/array	p value	Corrected p value
CL #1: uniform upregulation, 61 KEGG-annotated genes
04910	Insulin signaling pathway	9/163	1.06E−4	0.010
00010	Glycolysis/gluconeogenesis	6/78	2.81E−4	0.026
00030	Pentose phosphate pathway	4/29	3.36E−4	0.031
05120	Epithelial cell signaling in Helicobacter pylori infection	5/79	2.25E−3	0.205
04370	VEGF signaling pathway	5/91	4.17E−3	0.379
CL #2: PT specific upregulation, 20 KEGG-annotated genes
04910	Insulin signaling pathway	5/163	3.11E−4	0.015
04920	Adipocytokine signaling pathway	3/83	3.75E−3	0.184
CL #3: uniform downregulation, 46 KEGG-annotated genes
00980	Metabolism of xenobiotics by cytochrome P450	7/65	1.31E−6	9.83E−5
00190	Oxidative phosphorylation	6/114	4.62E−4	0.035
05216	Thyroid cancer	3/37	4.16E−3	0.312
05010	Alzheimer's disease	3/38	4.49E−3	0.337
03320	PPAR signaling pathway	4/77	4.64E−3	0.348
00072	Synthesis and degradation of ketone bodies	2/12	4.84E−3	0.363
CL #4: upregulation in the PT and HRT, 22 KEGG-annotated genes
00710	Carbon fixation	3/29	2.28E−4	0.008
04720	Long-term potentiation	4/85	3.99E−4	0.015
05214	Glioma	4/88	4.56E−4	0.017
05040	Huntington's disease	3/48	1.03E−3	0.038
04020	Calcium signaling pathway	5/212	1.65E−3	0.061
04912	GnRH signaling pathway	4/132	2.09E−3	0.077
00010	Glycolysis/gluconeogenesis	3/78	4.15E−3	0.154
04916	Melanogenesis	6/127	1.06E−5	3.93E−4
04910	Insulin signaling pathway	6/163	4.43E−5	0.002
CL #6: downregulation in the PT and HRT, 3 KEGG annotated genes
04360	Axon guidance	2/174	3.31E−3	0.007

Fig. 3

Main functional categories found to be affected in postmenopausal women. The number of genes is presented in non-redundant manner, i.e., “cytoplasm” includes a total of 29 genes from which nine are also included in the “response to contraction” and four in “energy metabolism” and therefore excluded from the “other cytoplasm-related” genes. The same approach was used in presenting the number of genes in all functional categories to avoid overpresenting the numbers. ↑ upregulation, ↓ downregulation

Power training-specific functional categories

The PT-specific gene clusters were CL #2 and CL #8, which consisted of 82 upregulated and 20 downregulated genes, respectively (Fig. 1b, h, Table 1). CL #2, which was specifically upregulated in the PT group and downregulated in the HRT and CO groups (Fig. 1b), including 13 genes without GO or KEGG annotations, 29 genes with annotation “cytoplasm,” and 9 genes with annotation “actin binding” (Table 2). GO term “cytoplasm” constitutes 42% of the whole cluster and refers to all cellular actions taking place within the cytoplasmic entity of cells excluding the cell membrane and the nucleus. One of the most highly upregulated genes in CL #2 was the STARS gene, which belongs to both “cytoplasm” and “actin binding” GO terms. STARS expression was 50% upregulated among the PT women, whereas the expression was 16% downregulated within the HRT and 46% downregulated within the CO group (Supporting Table S1). In addition, five genes encoding kinases, five encoding proteins related to apoptosis or proteolysis, four encoding signal transducers, four encoding proteins involved in regulation of transcription, and three encoding structural proteins were included in “cytoplasm” and/or “actin binding” GO terms. Moreover, CL #2 was significantly enriched with “insulin” and “adipocytokine signaling” KEGG pathways (Table 3). The genes contributing most to the significance of these upregulated pathways in the PT group were MKNK2, GLUT4, FASN, AKT2, PRKAR2A, and STK11. The enrichment analysis did not reveal any statistically significant functional categories from CL #8. This is probably due to the fact that 40% of the genes in this cluster were unclassified having neither GO nor KEGG annotation.

HRT-specific functional categories

The HRT-specific gene clusters were CL #5, which consisted of 45 upregulated genes, and CL #7, which included 31 downregulated genes (Fig. 1e, g, Table 1). Eleven genes did not have GO or KEGG annotation. According to enrichment analysis (Table 2), the use of HRT significantly upregulated transcription related to GO term “mitochondrion” and downregulated transcription of genes involved in “regulation of cell growth”. We did not find any significant HRT-specific enrichment of KEGG-annotated pathways. Genes from “mitochondrion” included four genes encoding mitochondrial ribosomal proteins (MRPS12, MRPS36, MRPL27, and MRPL33), one encoding mitochondrial chaperon (HSPE1), four genes directly involved with energy metabolism through oxidative phosphorylation (ATP5H and COX7A2), citric acid cycle (IDH3A) or fatty acid oxidation (PECI), and three genes involved with other mitochondrial functions (MTP18, PTPMT1, and TST). CL #7 had one significant category, the “regulation of cell growth”, which included three genes (ING1L, CHAD, and UBE2E3).

Functional categories affected by both PT and HRT

Both interventions, that is, one year of either PT or HRT, resulted in parallel upregulation of 81 genes (CL #4) and downregulation of 43 genes (CL #6) in contrast to that observed in the CO group (Fig. 1d, f, Table 1). Nine KEGG pathways were significantly enriched among the upregulated genes (Table 3). These pathways included signaling related to carbohydrate metabolism (KEGG pathways 00710, 00010, and 04190 in Table 3) and calcium signaling (KEGG pathways 04020, 04720, 04912, 04916, 05040, and 05214 in Table 3). One of the most dramatically affected genes was ALDOA with expression over 70% upregulated in the PT, 20% upregulated in the HRT, and 10% downregulated in the CO group (Supporting Table S1). This enzyme catalyzes conversion of C6-carbohydrates to C3-carbohydrates being the initial step in the energy producing instead of energy consuming phase of glycolysis. Other upregulated genes contributing to carbohydrate metabolism through “glycolysis/glyconeogenesis” or “insulin signaling” were PKM2, which catalyzes the last step of glycolysis, TPI1, PHKG1, FOXO1A, MKNK2, CALM1, and CALM3. According to KEGG annotations CALM1, CALM3, and PHKG1 are also involved in the “calcium signaling pathway”. In addition, DCTN1, CaMK2A, CaMK2B, FZD7, and DVL1 genes, which encode proteins also involved in “calcium signaling”, were upregulated in the PT and HRT women.

The only statistically significant downregulated pathway among both PT and HRT women was “axon guidance” (Table 3), including two moderators of G-protein signaling, the RGS3 and SRGAP2. Only three of the 43 genes belonging to CL #6 had KEGG annotation and 19 out of 43 had GO annotation, leaving 24 genes of CL #6 without any annotations at all.

Functional categories, which were changed into same direction in all study groups

The largest and most heterogeneous clusters were CL #1 with 224 genes and CL #3 with 139 genes (Fig. 1a, c, Table 1). Transcription of genes in these clusters was either upregulated or downregulated within all three study groups, i.e., neither PT nor HRT was able to reverse the regulation of these genes. CL #1 with uniform upregulation is mainly formed of GO terms (Table 2) related to muscle development (including GO terms GO:0007517 and GO:0007519), muscle contraction (including GO:0006936, GO:0006941, and GO:0030731), and muscle energy metabolism (including GO:0005977, GO:0004368, and GO:0016600) and KEGG pathways (Table 3) related to energy (carbohydrate) metabolism (including pathways 04910, 00010, and 00030) as well as VEGF and epithelial cell signaling. It is notable that enrichment related to muscle carbohydrate utilization and especially glycolysis is present in several clusters. The PT specifically upregulated five “insulin signaling” genes, both PT and HRT upregulated six “glycolysis” or “insulin signaling” genes, and finally 16 genes that showed uniform upregulation within all study groups had “glycolysis”, “glycogen metabolism” or “insulin signaling” annotation.

Other significantly and uniformly upregulated genes, which also might have direct effects on muscle functionality, were classified as “muscle development” (eight genes) or “muscle contraction” (12 genes). Among these were two genes encoding myosin heavy chain polypeptides (MYH1 and MYH8). To avoid confusion, the MYH1 gene, which encodes for fast myosin heavy chain IIx polypeptide, will be referred to as MYH IIx in this report. Transcription of MYH IIx was substantially upregulated in all study groups (FC: 2.40, 1.85 and 1.74 in the PT, HRT, and CO, respectively). The relative proportion of MyHC IIx protein was also increased in all study groups according to gel electrophoresis analysis (FC: 1.9, 1.8 and 2.2 in the PT, HRT, and CO, respectively; data not shown). Transcription of MYH8, which encodes perinatal myosin heavy chain polypeptide, was hardly detectable as upregulated in the HRT group but substantially more so in the PT and CO groups (FC: 1.84, 1.02, and 1.79 in the PT, HRT, and CO, respectively). The rest of the genes in “muscle development”—related terms were myostatin, FHL3, MAPK12, MYOD1, MEF2D, MYF6, STK23, and MYLPF. These genes were upregulated in the PT group by 25–54% and to a lesser extent in the HRT and CO groups (3–21%), except MYOD1, which had 30% upregulation in the CO group, 14% upregulation in the HRT group, and only 4% in the PT group (Supporting Table S1).

Of the genes downregulated in all groups (CL #3), the largest common terms were related to mitochondrial energy metabolism with 18 genes contributing to the significance of three functional categories (terms GO:0005739, 00190, and 00072 in Tables 2 and 3) and oxidation–reduction reactions with 14 genes belonging to one or more of the following categories: GO:0047115, GO:0016491, GO:0004033, GO:0047026, GO:0047042, GO:0032052, GO:0015721, and 00980 listed in Tables 2 and 3. The next largest significant functional category was “proteinaceous extracellular matrix” with nine downregulated genes. Also GO/KEGG categories: “regulation of muscle contraction”, “signaling through protein binding”, “lipid metabolism,” and “regulation of cell cycle” were affected.

Discussion

This study investigated global transcriptional changes following power type of training and estrogen containing hormonal replacement in skeletal muscle of postmenopausal women. We found 665 differentially expressed genes, which were assorted into eight non-redundant clusters. PT appeared to have broader effects on muscle metabolism, as seen in the number of affected genes (510), while effects of HRT were narrowed down to 54 changed genes. Moreover, the “without treatment”—condition led to the discovery of 161 genes with a significantly changed expression pattern. Using enrichment analysis, we were able to identify transcript groups, which were specific to PT, specific to HRT, shared with PT and HRT, and finally common to all studied, 50–57-year-old postmenopausal women (Fig. 3).

Specific transcription responding to muscle contraction

PT specifically affected transcription of nine genes known to respond to muscle contraction. Among these was STARS, which is considered to be an important mediator linking changes in actin dynamics and sarcomere structure with muscle gene expression (Kuwahara et al. 2005). Recently, Lamon et al. (2009) showed that resistance training stimulates the expression of STARS, which was associated with skeletal muscle hypertrophy. Here we show, that also plyometric power training, which improves muscle power and composition (Sipilä et al. 2001; Taaffe et al. 2005) stimulates the expression of STARS. Our finding indicates that STARS expression is induced by contraction per se and that the signaling through STARS can be directed to different pathways either to directly promote hypertrophy (Lamon et al. 2009), or to improve muscle composition and neuromuscular function (Taaffe et al. 2005). The components of such pathways are most likely co-regulated and therefore correlate closely with expression of STARS. In our data, the change in the mRNA expression of STARS did not correlate with change in any of the downstream genes investigated by Lamon et al. However, we find high correlation with other genes. Therefore, the ten most promising novel candidate genes (r = 0.54–0.68, p < 0.00001) for STARS-related regulation of muscle composition and function could be STIM, MRPL30, PPM1A, AFAP1L1, SYNE2, LONRF1, ETFA, AGL, and FHOD1.

PT and HRT have overlapping effects on the carbohydrate metabolism-related transcription

Some mediators of the insulin signaling pathway were upregulated solely by PT (five genes) and some by both PT and HRT (six genes), suggesting overlapping effects of these two interventions. PT specifically upregulated the mRNA expression of AKT2 and several other kinases as well as the expression of GLUT4. AKT protein kinases have three isoforms from which AKT2 is most important in modulating glucose homeostasis by regulating GLUT4 translocation into the plasma membrane in skeletal muscle (Gonzalez and McGraw 2009). Our result may indicate that the import of glucose into muscle cells is increased in the PT women assuming that the enhanced mRNA expression will actually lead to a larger amount of GLUT4 proteins. Ten additional genes related to the insulin signaling pathway were upregulated among all study participants. These include OGT, which has been recently identified as a suppressor of insulin signaling and its hepatic overexpression causes insulin resistance and dyslipidemia (Yang et al. 2008). The upregulation of the OGT gene was most prominent within the CO women (22%), and less obvious among the PT (12%) and HRT (6%) women (Supporting Table S1). OGT catalyzes the attachment of O-GlcNAc to proteins providing post-translational modification comparable and often competitive to phosphorylation. The rate of O-GlcNAcetylation is highly dependent on the availability of intracellular glucose, from which 2–5% is converted into UDP-GlcNAc, the donor substrate for GlcNAcetylation (Marshall et al. 1991). Therefore, OGT functions as an energy sensor, capable of regulating signal transduction, transcription, and proteosomal degradation (Wells et al. 2003; Copeland et al. 2008). In addition, many genes directly related to glycolysis were upregulated either in both PT and HRT groups (three genes) or in all three study groups (seven genes). The mRNA expressions of all major enzymes (PGM1, PFKM, FBP2, ALDOA, and TPI1) needed for conversions of glucose to C3-carbohydrates were upregulated at the mRNA level. Also the mRNA expression of PKM2, which is a glycolytic enzyme that catalyzes the transfer of a phosphoryl group from phosphoenolpyruvate to ADP thus generating ATP (Ikeda and Noguchi 1998), was upregulated in the PT and HRT groups. Taken together, it seems that the tendency to increase transcription related to carbohydrate metabolism after menopause is further enhanced by PT and HRT and finally manifests as better muscle performance. It is also possible that after menopause, the excess of sugar metabolites starts to accumulate, which produces a need to enhance both insulin signaling and glycolysis. Further studies are needed to verify if this strain observed here at mRNA level is due to an increment in energy need per se, due to truncation in energy production or due to inefficiency in energy utilization after menopause and whether the increments at mRNA level will lead to actual increments of corresponding proteins/enzymes.

Parallel activation of calcium signaling pathway by PT and HRT

Yet another candidate pathway for improving muscle performance is the calcium signaling pathway. Myoplasmic calcium (Ca2+) is a well-known secondary messenger involved in excitation–contraction coupling (Melzer et al. 1995) and is believed to be involved in an array of cellular functions including carbohydrate metabolism, gene transcription, protein synthesis, and mitochondrial biogenesis, further increasing the diversity of the pathway (Hook and Means 2001; Wu et al. 2002; Wright 2007; Rose and Richter 2008; Illario et al. 2009; Rose et al. 2009). In our study, five members of the calcium signaling pathway, namely two calmodulins: CALM1 and CALM3, and three calmodulin and/or Ca2+-activated kinases: PHKG1, CaMK2A, and CaMK2B, showed consistent transcriptional activation. Training increased expression of these genes by 20–30%, and the use of HRT by 0–15%, while the same genes were slightly downregulated in the CO group (2–13%, Supporting Table S1). It is well established that, contraction-induced elevation in Ca2+ concentration leads to rapid activation of CALMs, which in turn activate multifunctional CaMK kinases among other downstream targets (Walsh 1983). Our results suggest that the transcription of CaMK2, the most abundant CaMK in skeletal muscle, is regulated by physical exercise and that long-term PT leads to elevated basal mRNA levels. This is in line with protein level studies in which the activity of CaMK2 has been shown to be rapidly induced by exercise in an intensity dependent manner (Rose et al. 2005). It is also widely acknowledged that physical exercise increases the mRNA and protein expression of GLUT4. The mechanism for an exercise-induced increment is not entirely clear, but it has been suggested that it may go through an increase in intracellular Ca2+ content, which activates CaMK and calcineurin, which in turn activates specific transcription factors including MEF2 (Holmes and Dohm 2004). Our results support this hypothesis and suggest that PT also directly affects the transcription of CaMK, but not that of calcineurin, whose transcription was not induced in our study. The mechanism by which HRT is able to resist decrements in CALM1, CALM3, PHKG1, CaMK2A, and CaMK2B transcription, which was seen among the CO women, is not known.

HRT aids in maintaining muscle mitochondrial function via transcriptional regulation

One important and interesting aspect of aging-related phenomena discovered in previous microarray studies is the enhanced negative transcriptional regulation of oxidative phosphorylation and other mitochondrial functions with advancing age (Welle et al. 2003; Welle et al. 2004; Zahn et al. 2006; Melov et al. 2007; Welle et al. 2008). This downregulation was also observed in our study with early postmenopausal, middle-aged women. In our data, 29 genes with GO annotation “mitochondrion” or related to “oxidoreductase activity” were downregulated in all study groups within the 1-year interval. However, downregulation of “mitochondrion” genes was only minimal among the HRT women being 0–13% for most of the genes and 16–23% for nine genes related to oxidation events (Supporting Table S1). Much more severe downregulation was observed in the PT (5–55%) and CO (3–46%) groups. Moreover, there was specific upregulation of another set of 12 “mitochondrion”—genes within the HRT women (Fig. 3). Even though the upregulation of these genes was modest (0–13%) among the HRT women, the same genes were downregulated substantially more among the PT (12–53%) and CO (2–33%) women (Supporting Table S1). This data indicates that HRT may be beneficial for maintaining proper functionality of mitochondria, which seems to start accumulating negative transcriptional regulation already in the early postmenopausal years. Furthermore, contrary to strength-based resistance (Melov et al. 2007) and endurance training (Mahoney et al. 2005), the plyometric power training used in the current study was not effective in resisting negative transcriptional regulation of mitochondria. We also had some muscle samples (n = 46) taken 6 months after the initiation of the Ex/HRT intervention, from which we determined the amount of mitochondrial DNA relative to genomic DNA in order to count the number of mitochondria. However, we did not observe differences in the number of mitochondria between the samples of PT, HRT, and CO women (data not shown), which indicates that the regulation is on functionality rather than the number of mitochondria.

Transcriptional regulation of myogenesis and muscle contraction-related genes in all study groups

We found that several functional categories directly related to myogenesis, muscle structure, and function were affected during the first years postmenopause (Fig. 3). Eight “muscle development”—genes were upregulated within all study groups. For example, from the four members of MYOD1-family of myogenic regulatory factors (MYOD1, myogenin, MYF5, and MYF6), both MYOD1 and MYF6 were found to be upregulated within all study groups. These transcription factors have both specific and redundant functions during muscle development, and their activation is required for satellite cell differentiation during myogenesis and muscle regeneration (Yablonka-Reuveni and Rivera 1994; Yablonka-Reuveni et al. 1999). Also, the MEF2 transcription factor family includes four genes, which have sequential expression patterns during muscle development. The MEF2 gene, which we found to be activated among all study participants, was MEF2D, whose binding to the promoter of MYOD1 is required for increase and maintenance of MYOD1 expression during differentiation into myotubes (L’honore et al. 2007). Taken together, it seems that transcriptional regulation of myogenesis and muscle regeneration is activated in the early postmenopausal women.

Furthermore, functional categories “regulation of contraction,” “muscle contraction,” and “extracellular matrix” were enriched among all study groups (including 11 upregulated and 14 downregulated genes). In addition, PT specifically affected another group of genes that are known to respond to contraction (nine upregulated genes). In studies by Welle and co-workers (Welle et al. 2003; Welle et al. 2004), MYH8 encoding the perinatal myosin heavy chain, was found to be upregulated with aging. Similarly, we found that MYH8 was upregulated within all study groups albeit the upregulation was barely visible within the HRT group (2%) but much more pronounced in the PT (84%) and CO (79%) groups (Supporting Table S1). Expression of MYH8 dominates the early skeletal muscle development, and it is re-expressed in regenerating muscles (d’Albis et al. 1988; Weiss et al. 1999). Regeneration requires the activation of quiescent satellite cells to proliferate and differentiate into myoblasts, which eventually fuse with existing muscle fibers to repair damage or to help fibers to cope with physiological demands. Even though the expression level of perinatal MYH8 was low compared to the adult MYH genes, MYH I, MYH IIa, and MYH IIx, the observed substantial increment in the expression in the PT and CO groups may indicate increased demand for myofiber regeneration. Furthermore, this demand appeared to be absent among women using HRT possibly because they were the only group able to maintain and even increase their muscle mass.

Strengths and caveats of the current study

Most of the expression array studies investigating the effects of aging (Welle et al. 2003; Welle et al. 2004; Giresi et al. 2005; Zahn et al. 2006) are cross-sectional studies and compare old (65–80 years old) to younger (19–30 years old) participants neglecting the middle-aged population. Furthermore, only one study has included women (Welle et al. 2004). Therefore, our study with early postmenopausal 50–57-year-old women extends the work undertaken in previous studies. We studied the transcriptional changes using the age group, which according to earlier studies (Kallman et al. 1990; Phillips et al. 1993; Samson et al. 2000) is just beginning to accumulate deteriorations in their muscle function. This age range is optimal for detecting the first alterations in transcriptome even before they are fully translated to the phenotype. The relatively small sample size may, however, affect the results. On the other hand, the long duration of our study enables us to examine true adaptation to treatments both at the phenotype and transcriptome level. For skeletal muscle tissue, 1 year is a rather long follow-up since muscle is highly adaptive to changes in, e.g., loading conditions. Furthermore, the RCT design with repeated measures analysis and using a single HRT preparation strengthens our study by eliminating physiological and genetic differences in responses to training or HRT and, therefore, enables the detection of true intra-individual responses.

To date, many studies have used quantitative PCR (qPCR) to confirm microarray results. For example, in 16 microarray studies using human muscle tissue as a source for RNA (Roth et al. 2002; Rome et al. 2003; Welle et al. 2003; Zambon et al. 2003; Welle et al. 2004; Wittwer et al. 2004; Mahoney et al. 2005; Teran-Garcia et al. 2005; Schmutz et al. 2006; Chen et al. 2007; Melov et al. 2007; Pöllänen et al. 2007; Sjogren et al. 2007; Thalacker-Mercer et al. 2007; Welle et al. 2008; Stepto et al. 2009) and indicating the use of qPCR as validation method, 2–9 genes have been re-analyzed using either the same or different RNA pool. In ten cases, the microarray results were totally confirmed, and in six cases, there were 1–2 failures (out of nine, eight, seven, six, and five genes studied). Importantly, all failures have been on genes with the low level of expression. This also holds true for the array platform used in the current paper as demonstrated by Melov et al (2007) and Pöllänen et al (2007). Building on this literature, the quality and reproducibility of microarray experiments are nowadays well established, and additional validation with qPCR seems to be unnecessary, at least as far as the lowest expression levels are not considered. A more valuable extension would come from analyzing, e.g., protein samples, which we unfortunately were unable to do for other than the MyHC proteins, due to lack of samples.

In addition to the molecular level data, we have conducted comprehensive muscle phenotype analyses including tests for muscle power, strength, performance, composition, and mass (Sipilä et al. 2001; Taaffe et al. 2005). The results obtained showed the effectiveness of the PT and HRT on these phenotypes, which we believe is due to the observed changes in muscle transcriptome.

Conclusions

This study showed that several functional gene categories, which may directly affect the preservation of neuromuscular performance after menopause, are already altered within the first years of the postmenopausal period and that some of these are modulated by 12-month power training and/or hormonal replacement. The newly observed transcriptional changes partially explain phenotypic improvements after PT or HRT (Sipilä et al. 2001; Taaffe et al. 2005) and provide novel insights for further studies. For example, contraction-induced gene regulation via STARS, transcriptional activation of ALDOA and other glycolysis enzymes, regulation of signal transduction via OGT and insulin, and maintaining functionality of mitochondria as well as regulation of adult myogenesis and myofiber regeneration are all important but not thoroughly studied mechanisms involved in the regulation of muscle power, performance and quality after menopause. Taken together, our main findings indicate that PT and/or HRT may function through improving energy metabolism and response to contraction as well as by preserving functionality of mitochondria in skeletal muscle of postmenopausal women.

Below is the link to the electronic supplementary material.

(XLS 8988 kb)

(XLS 29 kb)